amrex-astro / castro Goto Github PK

We changed the way we customize a problem setup. We need to write up all the customization stuff, including:

• problem_post_init
• tagging
• new derived fields

we do not trace source terms for ppm_type = 0 currently

once we merge the radiation code into Castro (being done on the unified_code branch), then we will rely on the helmeos in Microphysics.git instead of the previous helmeos_without_rad. The Microphysics helmeos will disable photon pressure when the code is compiled with USE_RAD=TRUE.
We need to document this.

Note that this also means that if you want radiation (say for neutrinos) and also want the radiation pressure to be included in helmeos, this is currently not supported.

This needs to be clearly documented

need to document runtime diagnostics

The function (Castro::getTempDiffusionTerm) that calculates the thermal diffusion term, div. sigma grad T, seems to always operate with S_old

So it seems that we are never time-centering this term in our source term predictor-corrector update strategy.

Is the solution as simple as changing

   MultiFab& S_old = get_old_data(State_Type);                

to

   MultiFab& S_new = get_new_data(State_Type);                

for the corrector?

Note that once we compute this term, it looks like we incorporate it into new_sources in the correct fashion for the correction.

In advance_hydro in Castro_advance.cpp:

Suppose DIFFUSION is defined and add_ext_src is 0.

Then we don't call get getSource routine but we do still add diffusion to ext_src_old and add ext_src_old to stateout

We also add the diffusion to ext_src_new but then we don't call time_center_source_terms

The driver ca_umdrv and ca_umdrv_rad should be merged, and probably made dimension-agnostic.

The 1-d radiation solver has separate arrays for ergdnv and lamgdnv in umdrv. But the 2- and 3-d use the Godunov indices to pack these into q1, ...

The HLLC solver approximates the flux on the interfaces directly. For Cartesian coords, the pressure term is part of the momentum flux, so there is no issue.

For axisymmetric coords, this is not the case. The solution is to either include it in the flux and then have an explicit (non-conservative) geometric source term or ignore it in the flux. In either case, we still need a valid edge-centered, time-centered pressure.

we need to document the options for angular momentum conservation. Some text for this is in Max's thesis

The radiation routines (RAD_?D.F) use a common block defined in StateLayout.H. These could instead get the same information from meth_params_module

There are an abundance of new quantities stored in the primitive variable vector, q, that are not implemented in CastroRadation.

QDPDR
QDPDE
QC
QGAMC
QCMSL

we probably don't want these in the primitive varaible vector (with the exception of QGAMC), since these are all reconstructed and traced, and I don't think we need that (alternately, we could create a QVAR_TRACE that defines the variable count up to which we do the reconstruction and tracing)

none of the problems in the regression test suite exercise the neutrino portions of the code

The new scheme for protecting against density undershoots, and castro.limit_hydro_fluxes_on_small_dens needs to be documented. Best to document all of the density reset stuff together

When we do a gravity sync after a reflux for Poisson gravity, we use a sync source to the energy that is non-conservative. Since the default in Castro is now to use the conservative energy scheme for the main gravity source term, we should do the same for the sync source.

The sync source is described in the first Castro code paper, it's the equations just after Equation 39. The momentum source term is straightforward and consistent with our current approach -- you just compute the 'corrected' source term accounting for the original gravitational force and the reflux, and then subtract off the term that was already added during the main advance. The energy source term for the sync is essentially computed by using the idea that Sr_E = v . Sr_mom, an idea that we use elsewhere when we're not concerned about strict conservation.

If we want this to be conservative in the same manner as the main source term, what we should do is to construct the full flux-based energy source term, using the final gravitational field and corrected fluxes, and then subtract off the energy flux that was contributed during the main advance.

I am not sure yet exactly what the sync source should look like or whether we even have enough information in post_timestep to do it; will follow up on that.

To make it easier to discover, document, and add runtime parameters, we should modify Gravity.cpp, Radiation.cpp, and Diffusion.cpp to use the _cpp_parameters functionality.

eta1, ... are not considered in the radhydro version of ctoprim

ppm_temp_fix is not implemented in 1-d, nor in any of the radiation tracing

It also needs to be documented.

In the C prototypes for the Fortran functions we call directly from Castro, we often declare the data pointer to the FAB as:

const Real * state

But this means that state is a pointer to a Real that is constant, i.e., we should not modify the data that state points to (this follows the read-it-backwards technique advocated here:
https://isocpp.org/wiki/faq/const-correctness#ptr-to-const

We should only be using const for arrays that we mark intent(in) in Fortran.

Otherwise, we should not use const.

Note: it seems that compilers are not enforcing this currently.

The 1-d PLM hydro (trace_1d.f90) does not implement reference states. The 2- and 3-d version does.

We should not allow negative energy by default -- this is probably not what users intend.

The sponge is applied by subtracting dt * (rho * u)^(n+1) / sponge_timescale from the momentum. This is not time-centered; we should apply it in a predictor-corrector fashion like other source terms. As a consequence of this change, the source terms to the hydrodynamics would also "see" the sponge, making the hydro update more accurate.

Formally, boxlib supports double and single precision, and can be compiled to support either via the PRECISION variable in the makefiles.

However, all of the Fortran routines that we write assume that we pass data through as double *, so only double precision will work.

radiation should be a full-fledged member of Castro, and we should merge the code from that repo into here as as starting point of eliminating duplication.

we need to document how we limit based on burning / species change. Some text for this is in wdmerger paper II draft.

In 1D/2D, the call to normalize_species_fluxes in consup is done as the first step, while in 3D, the call is done after the artificial viscosity has been added. This is inconsistent, and also suboptimal for 1D and 2D because the artificial viscosity may cause the species fluxes to no longer be normalized.

In Castro_setup.cpp, for the hybrid momenta I used the scalar_bc for setting physical boundary conditions. But I think this may be incorrect. In particular, radial momentum and angular momentum are position dependent, so something like a piecewise constant extrapolation may not give a meaningful value.

We store both rho edot and edot in the Reactions_Type state data. We should not need both

Some of the runtime diagnostics are not weighted for axisymmetric coords. E.g., in sponge_nd.f90,

             state(i,j,k,UMX:UMZ) = state(i,j,k,UMX:UMZ) / (ONE + alpha * sponge_factor)                                        

             state(i,j,k,UEDEN) = state(i,j,k,UEDEN) + HALF * sum(state(i,j,k,UMX:UMZ)**2) * rhoInv - ke_old                    

Currently, we dimension the primitive variable vector as:

#ifdef RADIATION
   double precision :: q(QRADVAR)
#else 
   double precision :: q(QVAR)
#endif

To make life easier, we should introduce new parameters:

• NQVAR will always be the pure hydro primitive variables (the primitive version of the NVAR state variables)
• NQRADVAR will be the radiation components, namely energy density for the groups, QPTOT and QREITOT
• NQ will be used for dimensioning q, and will be NQ = NQVAR + NQRADVAR

The source term primitive array, srcQ will still be dimensioned just as NQVAR

The implementation of the slipwall and noslipwall boundaries in Castro_setup.cpp appear to be identical. noslipwall should force the transverse velocity to 0.

The following error message is generated when creating the actual_network.f90 file when compiling Castro with Python 3.4:

write_network.py: creating actual_network.f90
write_network.py: working on network file ../../Microphysics/networks/general_null/gammalaw.net...
Traceback (most recent call last):
File "../../Microphysics/networks/general_null/write_network.py", line 329, in
write_network(networkTemplate, netFile, outFile)
File "../../Microphysics/networks/general_null/write_network.py", line 229, in write_network
fout.write(string.replace(line,"@@nspec@@", str(len(speciesList))))
AttributeError: 'module' object has no attribute 'replace'

Castro/Source/Src_1d/trace_ppm_1d.f90 doesn't seem to do source term tracing

We have > 40 problem setups now. We should arrange them into subdirectories by logical grouping. Like:

hydro_tests/
rad_tests/
science/

The source term predictor adds (dt / 2) * (dS/dt) to the source terms seen by the hydrodynamics. dS/dt is estimated using a backward difference over the last timestep, (S^{n} - S^{n-1}) / (t^{n} - t^{n-1}). However, in the code, the new-time and old-time source terms are actually the corrector and predictor, respectively, not the pure old-time and new-time values of the sources. So the derivative we calculate is not accurate.

Multilevel refinement breaks HSE in 1D models. This became a noticeable issue after the 16.08 release.

Level solves on fine grids are currently done using a relative error tolerance chosen by the user; a typical tolerance on level 0 for might be 1e-10 or 1e-11. Practical experience shows that this same tolerance is hard to achieve if you use many levels of refinement, and typically one has to loosen the tolerances on higher levels to achieve convergence.

A related issue is that if the refinement on a grid is not where the mass is, then the multigrid solve can easily fail to converge for the level solve(s) on the fine grid(s) because in that case the RHS is effectively zero, so we're solving Laplace's equation, yet the norm we compare against for the initial residual is based on the Dirichlet boundary conditions provided by the coarse grid. In practice this seems to be too hard of a problem to solve, and the error bottoms out at a level which in absolute terms is approximately machine roundoff level compared to the norm on the coarse grid.

A solution being considered is to use the maximum density on the domain (over all levels) as an absolute error tolerance, and to not use relative tolerances at all. This has the effect that the fine grids with no mass should be much easier to converge on; it also means that if the finest grid has the highest density and this ends up being smoothed out on the coarse grid, that the coarse grid solution will be slightly different because it will be easier to converge. However, this is not necessarily a bad thing, as it is not really physically meaningful for the coarse grid to be using a tolerance threshold based only on that level's density.

Need to document enthalpy diffusion, species diffusion, and viscosity

-- are there timestep limiters for all?

ppm_type = 0 doesn't work with radiation

We need to put the hydro on the GPUs

Things to consider:

• we should not pass array slides as arguments. We do this a lot with qaux as, e.g., qaux(:,:,:,QGAMC)
• we should do it piece by piece by moving as much of the hydro routines to be interfaced through C++

before max leaves...

some problem setups use this -- should document in the software section

In ppm_2d.f90, we do:

                    alpham = alpham*D2LIM/max(abs(D2),1.d-10)
                    alphap = alphap*D2LIM/max(abs(D2),1.d-10)

in the ppm_type == 2 section. The constant 1.d-10 seems arbitrary and may not be general

In Castro/Source/Src_nd/advection_util_nd.F90:src_to_prim, we calculate the source terms for the primitive variables in terms of the source terms for the conserved variables. Right now, we calculate the source term for the primitive variable QREINT (== rho e) by taking the source term for the total energy and subtracting the kinetic energy contribution. This is sub-optimal because it means sources like gravity directly affect the (rho e) equation. (It also pollutes the pressure source term, since the pressure source term is written in terms of the (rho e) source term.) Furthermore, there is no need to do this; since (rho e) is a conserved variable, there is a source term for (rho e) itself, src(i,j,k,UEINT), and this can be used as the primitive source term.

The call to ca_initdata, which is the user-facing routine intended to be used for the user to fill in the initial state corresponding to their problem setup, is not done with OpenMP. This means that if the user wants to speed up the initialization they need to code in the OpenMP directives by hand. We should make ca_initdata consistent with the rest of the code by threading it at the MFIter level with OpenMP, and tiling it. This should speed up initialization significantly.

We will have to fix the few problem setups that do write in the OpenMP by hand, and we will have to notify users that their ca_initdata implementation needs to be thread-safe and tiling safe (so that, e.g., they should not do something like state(:,:,:,:) = 0.0, and instead only work on the tiling region from lo(:) to hi(:)).

There are many problem setups that have their own hydrostatic boundary condition implementation. We should make a centralized implementation for this that can be set through the inputs directly.

we do parabolic reconstruction on the sources for all primitive state variables but only implement the source tracing for the momentum (velocity) sources

We currently call computeTemp near the very end of Castro::post_timestep(), notably after the final average down and after the call to sum_integrated_quantities. This seems incorrect to me; we should do it before the average down (even though it should give the same answer) and before the diagnostics (in case they rely on a current temperature evaluation).

In Diffusion.cpp, we treat the upper boundary as:

            if (phys_bc->hi(dir) == Symmetry   ||                                              
                phys_bc->hi(dir) == SlipWall   ||                                              
                phys_bc->hi(dir) == NoSlipWall ||                                              
                phys_bc->hi(dir) ==  Inflow    ||                                              
                phys_bc->hi(dir) == Outflow)                                                   
            {                                                                                  
              mg_bc[2*dir + 1] = MGT_BC_NEU;                                                   
            }                                                                                  
            else if (phys_bc->hi(dir) ==  Inflow)                                              
            {                                                                                  
              mg_bc[2*dir + 0] = MGT_BC_DIR;                                                   
            } else {                                                                           
              BoxLib::Error("Failed to set hi mg_bc in Diffusion::make_mg_bc" );               
            }

Note that Inflow is specified in both the if and the else. I think it should only be in the else -- this would be consistent with the lo face version.

There are multiple instances of do_rad_sphere in the radiation solver. These should be removed and instead located in a problem-specific routine