we need to standardize the naming for the impl namespace for each problem

Curerntly, the code to plot the mesh is part of the generation. we need to extract it and makeing it self-contained.
So that generating the mesh and visualizing it are separate steps/scripts

we currently select which order to use (1st, or 5th weno) based on the stencil size of the mesh provided. This is constraining a bit, since the 7 point stencil of the weno can support 1st order too.
So we need to add a way to select the order and then check that the selected order is compatible with the mesh provided.

1D

• linear advection
  • rhs (1st/weno3/weno5)
  • first-order Jacobian
  • OpenMP
  • Python bindings
• diffusion reaction
  • rhs
  • Jacobian
  • OpenMP
  • Python bindings
• 1d Euler
  • rhs (1st/weno3/weno5)
  • first-order Jacobian
  • OpenMP
  • Python bindings

2D

• 2d euler
  • rhs (1st/weno3/weno5)
  • first-order Jacobian
  • OpenMP
  • Python bindings
• 2d swe
  • rhs (1st/weno3/weno5)
  • first-order Jacobian
  • OpenMP
  • Python bindings
• 2d diffusion-reaction
  • rhs
  • Jacobian
  • OpenMP
  • Python bindings
• 2d compre nav-sto
  • rhs (1st/weno3/weno5)
  • first-order Jacobian
  • OpenMP
  • Python bindings

3D

• 3d euler
  • rhs (1st/weno3/weno5)
  • first-order Jacobian
  • OpenMP
  • Python bindings
• 3d compre nav-sto
  • rhs (1st/weno3/weno5)
  • first-order Jacobian
  • OpenMP
  • Python bindings

currently assuming first and last cells are on boundary, but this does not work for sample mesh

• min/max for mesh coordinates
• boolean to check if graph has periodic BC built in
  others?

i get a segfault from the test

• the simple meshing file needs to be extented to support a step inside

• the mesh for the domain wth the step should still be uniform

Test problems below are partially taken from here
A list of specific references for each is shown below, and will be updated.

1D

• linear advection

  • use domain:[-1, 1]
  • IC either provided by user or default one is \phi = sin(\pi x)
  • BC: periodic

• Euler Smooth

  • use domain: [-1, 1]
  • default IC is \rho =1 + 0.2*sin(\pi x), u= 1, p = 1
  • BC: periodic
  • analytical density at time t : \rho =1 + 0.2*sin(\pi (x-t))

• Sod problem

  • use domain: [-0.5, 0.5]
  • IC:
    • x<=0: \rho =1, u= 1, p = 1
    • x>0 : \rho =0.125, u= 1, p = 0.1
  • BC: homogeneous Neumann

• Lax problem

  • use domain: [-5, 5]
  • IC:
    • x<=0: \rho = 0.445, u = 0.698, p = 3.528
    • x>0 : \rho = 0.5, u = 0., p = 0.571
  • BC: homogeneous Neumann

2D

• Euler Smooth

  • use domain: [-1, 1]^2
  • default IC is \rho =1 + 0.2*sin(\pi (x+y)), u = 1, v = 1, p = 1
  • BC: periodic
  • analytical density at time t : \rho =1. + 0.2*sin(\pi (x+y- 2 t))

• Sedov problem (full version, no symmetry use)

  • reference
  • use domain: [-0.6, 0.6]^2
  • BC: homogeneous Neumann on all boundaries
  • IC is a high pressure concentrated small spherical region of radius R = 2. * std::min(dx, dy)
    • r<=R: \rho =1, u = 0, v = 0, p = (\gamma-1)/(pi*R*R)
    • r>R : \rho =1, u = 0, v = 0, p = 5e-5
  • dx, dy are the cell widths and the factor of 2 is used to spread the source over 2 cells for numerical reasons, and r = sqrt(x^2+y^2) is the distance from the origin.

• Sedov problem (exploit symmetry to simulate only one quadrant)

  • use domain: [0, 1.2]^2
  • BC: reflective on y=0 and x=0 and homogeneous Neumann for x=1.2 and y=1.2
  • IC involves a small spherical region of radius R = 3. * std::min(dx, dy) with high pressure
    • r<=R: \rho =1, u = 0, v = 0, p = ((\gamma-1)* 0.851072)/(pi*R*R)
    • r>R : \rho =1, u = 0, v = 0, p = 2.5e-5
  • dx, dy are the cell widths, and r = sqrt(x^2+y^2) is the distance from the origin.

• Riemann problem

  • use domain: [0.0, 1.0]^2
  • BC: homogeneous Neumann on all boundaries
  • various initial conditions are supported:
    • IC 1:
      • x>=0.5, y>=0.5: \rho = 0.5313, u = 0, v = 0, p = 0.4
      • x<0.5, y>=0.5 : \rho = 1, u = 0.7276., v = 0, p = 1
      • x<0.5, y<0.5. : \rho = 0.8, u = 0, v = 0., p = 1
      • x>0.5, y<0.5. : \rho = 1., u = 0, v = 0.7276, p = 1
    • IC 2:
      • x>=0.8, y>=0.8 : \rho = 1.5, u = 0, v = 0, p = 1.5
      • x<0.8, y>=0.8 : \rho = 0.5323, u = 1.206, v = 0, p = 0.3
      • x<0.8, y<0.8 : \rho = 0.138, u = 1.206, v = 1.206, p = 0.029
      • x>0.8, y<0.8 : \rho = 0.5323, u = 0, v = 1.206, p = 0.3

• Double Mach reflection

  • in progress

• Wind tunnel with step

  • in progress

3D

• Euler Smooth

  • use domain: [-1, 1]^3
  • default IC is \rho = 1 + 0.2*sin(\pi (x+y+z)), u = 1, v = 1, w = 1, p = 1
  • BC: periodic
  • analytical density at time t : \rho =1 + 0.2*sin(\pi (x+y+z-3 t))

• Sedov problem (exploit symmetry to simulate only 1/8)

  • use domain: [0, 1.2]^3
  • BC: reflective on planes next to source, homogeneous Neumann for the others
  • IC involves a small spherical region of radius R = 3. * min(dx, dy, dz) with high pressure
    • r<=R: \rho =1, u = 0, v = 0, p = (3 (\gamma-1)* 0.851072)/(4 pi R^3)
    • r>R : \rho =1, u = 0, v = 0, p = 2.5e-5
  • dx, dy, dz are the cell widths, and r = sqrt(x^2+y^2+z^2) is the distance from the origin.

TODO:

• there is a need of the name change for a Python script which builds dynamic library with CMake
• name build confuses Python interpreter and python -m build is not using build package, to actually build the , but is running the script
• changes needs to be done in CI as well

need to add a python script to visualize results

• 1D cases
• 2D cases

Goal: self-explanatory

Steps:

• make a new branch from develop
• look inside pressio-demoapps/docs/src
• doc is written in rst
• fill math formulation for each problem
• once done, open PRs for individual problems, or for example groups of problems (e.g. one PR for 1d, one for 2d ,etc)

Building the documentaiton:

cd pressio-demoapps/docs
make html

then open index.html in a browser. If you rebuild, you can refresh the page.
To build, you need sphinx and some other things.

we need to ensure that # of nonzeros after initialization is the same as after first evaluation of jacobian.
this is because for Eigen, one can add elements, so we need to make sure we don't add new elements.

• advection velocity for 1d, but need to change impl since it assumes unitary velocity
• gas constant for euler
• shock angle and initial shock Mach for double mach
• all coefficients for shallow water (depends on #12)

in 2d,3d we first handle the boundary cells, and then the inner ones.
We should do similarly for 1d.
That allows us to use the auxiliary stencil vector for boundary cells but just use the state (without copying data) for the inner cells.

• 1st order local Jacobians for Rusanov Flux
• Convert tests to gtest format
• Set up complex step derivative in tests
• Global Jacobian contruction
  • Dirichlet BCs
  • Neumann BCs
  • Reflection BCs
  • Periodic BCs
• Weno reconstruction contributions
• 1st order local Jacobians for Roe Flux

• update the Eigen-based swe problem class
• update the corresponding public include file
• add test for it

Currently only support evaluation of the RHS, need to add FD Jacobians to all classes

We should add diffusion such that it is consistent with the FV formulation.

Need to change:

• meshing script: should be pretty easy to change mesh scripts to do this
• mesh class in the code (now periodic means fully periodic along all axes)

need to better constraint the templates in weno.hpp

seems we have everything right except for top boundary

Using same BCs as first and WENO3, still gets unstable right away near wedge location.
Currently, the code throws if user selects WENO5 for double mach.

CMakeLists links main1d.cc to both the 2d and 3d tests; making an issue so we don't forget.

Modify script to enable non periodic BC

Add Roe flux capability.

Following chapter 11 of "Riemann Solvers and Numerical Methods" by Toro, the algorithm should

1. Compute Roe averaged velocities, total enthalpy, and speed of sound
2. Compute Roe averaged eigenvalues using an analytical expression, potentially with an "entropy fix".
3. Compute the Roe averaged right eigenvectors using analytical expressions for each entry
4. Project the jump in each conserved quantity onto the right eigenvectors. Again, there are analytical expressions for this in Toro's book or Roe's original paper.
5. Compute fluxes using equations 11.27-11.29 in Toro.

All of this should be possible following the same API used for the Rusanov flux in pressiodemoapps/impl_apps/eulerCommon/fluxes.hpp

Important note: there is a typo in Roe's original paper! The right hand side of Equation 22b should be \delta u_4 - w * \delta u_1.

Since we have specific problems, we need to add inside the impl classes checks such that when the problems are constucted , the mesh object passed in satisfies the physical bounds for that specific problem

e.g., for the Euler2d:

explicit EigenApp2d(const mesh_t & meshObj, int icId = 1)  
: meshObj_(meshObj), icId_(icId)  {    
 if (probid == 1){ 
   // throw if domain not compatible with 2d Sedov, else throw
 }
else if (probid==2){
   // throw if domain not compatible with 2d Riemann 
}
 ///
}

The operator () for pybind::array_t in line 906 of this file is not efficient. I added this myself originally since natively array_t did not support it.
One reason behind the perf issue is that looking at the impl (traces back to the parent array class, see that same file), it still performs a bound check.

Problem is that inside pressio-demoapps, we are using this to access the graph entries. We should change things to use the mutable and mutable_unchecked that seem to perform well since these proxy classes do not perform bounds check. However, we need to figure out how to handle these proxies properly.

• related to #6
• function should accept prim vars on one side and return the prim on the other
• a test should be added to only test this function is correct

need to add action for testing.
since we have C++ and Python code, the CI should cover both.
Check the readme for how to build/run the C++ tests and bindings.

Need to add for:

• explicit time stepping
• implicit

References:

• https://homepages.dias.ie/jmackey/jmac/node10.html

• others?

• problem setup: domain, BC, etc

• add initial condition for the problem here: https://github.com/Pressio/pressio-demoapps/blob/main/include/impl/euler2d/initial_condition.hpp

• add Eigen implementation:

  • file to use : https://github.com/Pressio/pressio-demoapps/blob/main/include/impl/euler2d/eigen_app.hpp
  • use a proper value for probid (see the other cases) to identify this particular problem

• add the public alias for this problem here: https://github.com/Pressio/pressio-demoapps/blob/main/include/euler.hpp

• add a test for it, use the Sedov test as template: https://github.com/Pressio/pressio-demoapps/tree/main/tests/sedov2d

  • to create the mesh files, see python script here
    If we use the domain in the reference above, the command is something:

  python create_full_mesh.py --outDir <where-youwant-the-mesh-files> \
  -s 7 --bounds 0.0 0.325 0.0 1.0 --periodic false -n 260 80

  the -s 7 is to create a mesh with a 7 point stencil.

should not do anymore python setup.py install but sometihing else.

in fact I get:

/usr/local/lib/python3.9/site-packages/setuptools/command/install.py:34: SetuptoolsDeprecationWarning: setup.py install is deprecated. Use build and pip and other standards-based tools.
  warnings.warn(

• make clear how dofs are stored, the difference between primitive and convervative, etc

need to change how constructors are exposed now. We have a mix of things: inside includable headers we have some that pretty much as just variadic and forward everything to problem classes. We need to expose public constructors ONLY in the includable headers.

seems like cmd line option passed to install command does not work

missing for 1d

Add the shallow water equations along with a case setup to pressio demo apps. To do:

• Add code for SWE equations
• test for the velocity function in the case of full mesh
• test for the velocity function in the case of sample mesh
• test for the full time integration
• validate against existing SWE solver

Obj: to create problem classes in python that use the c++ implementation

how:

• follow how we did bindings for pressio-tools or pressio4py
• should be relatively easy to do , maybe some issues can come up related to sparse jacobians

Priorites:

• Sedov problem for John to use