amrex-astro/maestro

A low Mach number stellar hydrodynamics code

maestro's Issues

wdconvect stopped early

I'm running wdconvect problem with OMP activated only. The input file is inputs_3d_C.128.
But the simulation seems very unstable and stops at some point occasionally. It gives me the following message: amrex::Abort::0::ERROR: ncenter invalid in Diag() !!!

moving initial model stuff to its own repo

The Util/initial_models/ code should not live in Maestro, since it is used by Castro and now Maestro-Ex too. I've created a new repo, https://github.com/AMReX-Astro/initial_models to hold this code. It only requires FBoxLib and Microphysics to build. It should take over all the functionality. In a bit I'll do a PR to remove the initial model code from Maestro.

add rotation

This is a long-term project (probably requiring a postdoc).

Maestro uses a 1-d base state to represent HSE and evolves the departures from this. For a rapidly rotating star, this will breakdown -- there is not a 1-d HSE model.

Instead we will need to do a 2-d base state, r, theta, that describes the structure. Generating the initial model itself will be challenging, but a relaxation technique and perhaps based on a self-consistent field method as done with binary pairs (that is partially implemented in Castro) would work.

An idea is that we initialize a spherical star, evolve a step, reaverage the Cartesian state to r, theta, resetting the base state, and then evolve again. Once equilibrium is reached, we reset the calculation with this initial calculation and turn on the reactions, etc. to play the simulation forward.

The change to the code would be rather invasive -- we will need to change all references to the base state to take a 2-d array. It is not clear how to extend the current averaging method to this case (although it may actually be simpler), and it is not clear how to include a w_0 term since there will no longer be a simple elliptic equation that drives the expansion.

The second complication is that the constraint itself will not be able to be written as a projectionable operation. Since there will be no equivalent to beta_0 that can hide the density / pressure stratification inside the divergence constraint. This is similar to the issue discussed in Klein & Paulius 2012 for the general EOS (that we get around with an average gamma1 with perhaps a lagging). SDC might be a path forward here, using a lagged U in the grad p0 term now. But this will require mathematical analysis.

(GPU) Invalid global read in burner_loop_3d_sph

Using: PGI 16.10, CUDA Toolkit 8.0
Problem: wdconvect_urca in SCIENCE
Network: ignition_reaclib/URCA-simple
Device: NVIDIA GTX 960, driver version 375.20
Branch: microphysics

cuda-memcheck complains about the following:

========= CUDA-MEMCHECK
========= Invalid __global__ read of size 8
=========     at 0x000001d8 in /home/dwillcox/codes/astro/MAESTRO/Exec/SCIENCE/wdconvect_urca/../../../Microphysics/networks/burner.f90:46:burner_module_burner_
=========     by thread (95,0,0) in block (31,0,0)
=========     Address 0x00000000 is out of bounds
=========     Device Frame:burner_loop_3d_sph_937_gpu (burner_loop_3d_sph_937_gpu : 0x1530)
=========     Saved host backtrace up to driver entry point at kernel launch time
=========     Host Frame:/lib64/libcuda.so (cuLaunchKernel + 0x2c5) [0x204235]
=========     Host Frame:/opt/pgi/linux86-64/2016/cuda/7.0/lib64/libcudart.so.7.0 [0x131b0]
=========     Host Frame:/opt/pgi/linux86-64/2016/cuda/7.0/lib64/libcudart.so.7.0 (cudaLaunch + 0x143) [0x2d653]
=========     Host Frame:/opt/pgi/linux86-64/16.10/lib/libaccnc.so (__pgi_uacc_cuda_launchk2 + 0x19ec) [0x331c]
=========     Host Frame:./main.Linux.PGI.acc.exe [0x3f0690]
=========     Host Frame:./main.Linux.PGI.acc.exe [0x3e9ec7]
=========     Host Frame:./main.Linux.PGI.acc.exe [0x3e47c5]
=========     Host Frame:./main.Linux.PGI.acc.exe [0xf3745]
=========     Host Frame:./main.Linux.PGI.acc.exe [0x43b2ea]
=========     Host Frame:./main.Linux.PGI.acc.exe [0x1d0059]
=========     Host Frame:./main.Linux.PGI.acc.exe [0x9bd4]
=========     Host Frame:/lib64/libc.so.6 (__libc_start_main + 0xf0) [0x20700]
=========     Host Frame:./main.Linux.PGI.acc.exe [0x9ab9]

The relevant lines of the burner module are:

    do n = 1, nspec
       state_in % xn(n) = x_in(n)
    enddo

The error message indicates the GPU is attempting to read from global memory for x_in, which contradicts the -Minfo compiler message for the calling subroutine which states:

937, Local memory used for rhowdot,x_in,x_out

From Mat on the PGI user forum, however:

By default in OpenACC, scalars are private and arrays are shared. In cases such as your code when arrays are reused, having them shared across all threads would cause a race condition. Hence, the user must privatize the array. The array is private to the loop so a "private" clause on a "vector" loop means every thread gets it's own copy of the array. A "private" on a "worker" loop means that it will be private to each worker but shared between the vectors within that worker. A "private" on a "gang" loop means each gang will get it's own copy which are shared amongst all the vectors within that gang.

Since privatized arrays are placed in global device memory and may lead to poor memory access patterns, it's advisable to avoid using "private" if the algorithm can be rewritten.

So why does the compiler report using local memory for rhowdot, x_in, x_out if privatized arrays should be in global device memory?

merge the microphysics branch

We want Maestro to use Microphysics, as does Castro. Most of the networks are working as expected (accounting for changes to the integratoin). The xrb_simple problem has some issues, probably related to T integration.

Use CUDA EOS and reactions

I'm working on swapping out the burner calls in react_state to use the CUDA-accelerated integration in the cudadevice branch of Microphysics.

Because the reaction integrator calls the eos, it must be a device subroutine, so I'm also replacing all the eos calls in MAESTRO to call the CUDA interface to the eos.

Maestro fails to compile with latest AMReX

We get the following errors when we try to compile Maestro (note, this also affects compiling the tools like fcompare:

mpicxx  -std=c++11 -Wall -g -O2 -ftree-vectorize -DBL_Linux -DFORTRAN_BOXLIB -DBL_FORT_USE_UNDERSCORE -I/home/zingale/development/AMReX//Src/Base  -c -o t/Linux.gfortran.mpi/o/AMReX_CArena.o /home/zingale/development/AMReX//Src/Base/AMReX_CArena.cpp
In file included from /home/zingale/development/AMReX//Src/Base/AMReX_Array.H:6:0,
                 from /home/zingale/development/AMReX//Src/Base/AMReX_CArena.H:9,
                 from /home/zingale/development/AMReX//Src/Base/AMReX_CArena.cpp:5:
/home/zingale/development/AMReX//Src/Base/AMReX_SPACE.H:27:2: error: #error AMREX_SPACEDIM must be defined
 #error AMREX_SPACEDIM must be defined
  ^~~~~
/home/zingale/development/AMReX//Src/Base/AMReX_SPACE.H:31:2: error: #error AMREX_SPACEDIM must be either 1, 2, or 3
 #error AMREX_SPACEDIM must be either 1, 2, or 3
  ^~~~~
/home/zingale/development/AMReX//Src/Base/AMReX_SPACE.H:47:26: error: ‘AMREX_SPACEDIM’ was not declared in this scope
     const int SpaceDim = AMREX_SPACEDIM;
                          ^~~~~~~~~~~~~~
/home/zingale/development/AMReX//Src/Base/AMReX_SPACE.H:47:26: note: suggested alternative: ‘AMREX_SPACE_H’
     const int SpaceDim = AMREX_SPACEDIM;
                          ^~~~~~~~~~~~~~
                          AMREX_SPACE_H
In file included from /home/zingale/development/AMReX//Src/Base/AMReX_CArena.H:9:0,
                 from /home/zingale/development/AMReX//Src/Base/AMReX_CArena.cpp:5:
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:51:19: error: ‘AMREX_SPACEDIM’ was not declared in this scope
     std::array<T*,AMREX_SPACEDIM> GetArrOfPtrs (std::array<T,AMREX_SPACEDIM>& a)
                   ^~~~~~~~~~~~~~
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:51:19: note: suggested alternative: ‘AMREX_SPACE_H’
     std::array<T*,AMREX_SPACEDIM> GetArrOfPtrs (std::array<T,AMREX_SPACEDIM>& a)
                   ^~~~~~~~~~~~~~
                   AMREX_SPACE_H
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:51:33: error: template argument 2 is invalid
     std::array<T*,AMREX_SPACEDIM> GetArrOfPtrs (std::array<T,AMREX_SPACEDIM>& a)
                                 ^
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:51:62: error: ‘AMREX_SPACEDIM’ was not declared in this scope
     std::array<T*,AMREX_SPACEDIM> GetArrOfPtrs (std::array<T,AMREX_SPACEDIM>& a)
                                                              ^~~~~~~~~~~~~~
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:51:62: note: suggested alternative: ‘AMREX_SPACE_H’
     std::array<T*,AMREX_SPACEDIM> GetArrOfPtrs (std::array<T,AMREX_SPACEDIM>& a)
                                                              ^~~~~~~~~~~~~~
                                                              AMREX_SPACE_H
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:51:76: error: template argument 2 is invalid
     std::array<T*,AMREX_SPACEDIM> GetArrOfPtrs (std::array<T,AMREX_SPACEDIM>& a)
                                                                            ^
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H: In function ‘int amrex::GetArrOfPtrs(int&)’:
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:53:34: error: invalid types ‘int[int]’ for array subscript
         return {AMREX_D_DECL(&a[0], &a[1], &a[2])};
                                  ^
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:53:41: error: invalid types ‘int[int]’ for array subscript
         return {AMREX_D_DECL(&a[0], &a[1], &a[2])};
                                         ^
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:53:48: error: invalid types ‘int[int]’ for array subscript
         return {AMREX_D_DECL(&a[0], &a[1], &a[2])};
                                                ^
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:53:17: error: there are no arguments to ‘AMREX_D_DECL’ that depend on a template parameter, so a declaration of ‘AMREX_D_DECL’ must be available [-fpermissive]
         return {AMREX_D_DECL(&a[0], &a[1], &a[2])};
                 ^~~~~~~~~~~~
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:53:17: note: (if you use ‘-fpermissive’, G++ will accept your code, but allowing the use of an undeclared name is deprecated)
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H: At global scope:
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:57:25: error: ‘AMREX_SPACEDIM’ was not declared in this scope
     std::array<T const*,AMREX_SPACEDIM> GetArrOfConstPtrs (const std::array<T,AMREX_SPACEDIM>& a)
                         ^~~~~~~~~~~~~~
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:57:25: note: suggested alternative: ‘AMREX_SPACE_H’
     std::array<T const*,AMREX_SPACEDIM> GetArrOfConstPtrs (const std::array<T,AMREX_SPACEDIM>& a)
                         ^~~~~~~~~~~~~~
                         AMREX_SPACE_H
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:57:39: error: template argument 2 is invalid
     std::array<T const*,AMREX_SPACEDIM> GetArrOfConstPtrs (const std::array<T,AMREX_SPACEDIM>& a)
                                       ^
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:57:79: error: ‘AMREX_SPACEDIM’ was not declared in this scope
     std::array<T const*,AMREX_SPACEDIM> GetArrOfConstPtrs (const std::array<T,AMREX_SPACEDIM>& a)
                                                                               ^~~~~~~~~~~~~~
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:57:79: note: suggested alternative: ‘AMREX_SPACE_H’
     std::array<T const*,AMREX_SPACEDIM> GetArrOfConstPtrs (const std::array<T,AMREX_SPACEDIM>& a)
                                                                               ^~~~~~~~~~~~~~
                                                                               AMREX_SPACE_H
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:57:93: error: template argument 2 is invalid
     std::array<T const*,AMREX_SPACEDIM> GetArrOfConstPtrs (const std::array<T,AMREX_SPACEDIM>& a)
                                                                                             ^
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H: In function ‘int amrex::GetArrOfConstPtrs(const int&)’:
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:59:34: error: invalid types ‘const int[int]’ for array subscript
         return {AMREX_D_DECL(&a[0], &a[1], &a[2])};
                                  ^
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:59:41: error: invalid types ‘const int[int]’ for array subscript
         return {AMREX_D_DECL(&a[0], &a[1], &a[2])};
                                         ^
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:59:48: error: invalid types ‘const int[int]’ for array subscript
         return {AMREX_D_DECL(&a[0], &a[1], &a[2])};
                                                ^
/home/zingale/development/AMReX//Src/Base/AMReX_Array.H:59:17: error: there are no arguments to ‘AMREX_D_DECL’ that depend on a template parameter, so a declaration of ‘AMREX_D_DECL’ must be available [-fpermissive]
         return {AMREX_D_DECL(&a[0], &a[1], &a[2])};
                 ^~~~~~~~~~~~
make: *** [/home/zingale/development/AMReX//Tools/F_mk/GMakerules.mak:112: t/Linux.gfortran.mpi/o/AMReX_CArena.o] Error 1

flame problem does not compile

The flame setup does not compile (but flame_1d does).

Update SDC Source to work with StarKiller Microphysics

The interface from Source_SDC in Maestro to the burner should be updated to work with the StarKiller Microphysics repository.

The SCIENCE/flame problem is a good test, compile with SDC=t

At present the interface is broken so compilation fails.

tiling

We should implement tiling into Maestro routines. This is probably feasible for anything that doesn't write to the base state at this time. VARDEN shows how to do the tiling in F90 BoxLib

resolution of g-modes

With the test_stability problem we found out, that maestro is not doing very well around peaks in the brunt-vaisala frequency (N). We attributed this to unresolved (reflected, interacting ....) g-modes in these areas.

Fabian Miczek proposed a much simpler test setup in his PhD thesis where one follows a single g-mode in a stably stratified ideal gas with constant N. Such a wave should move at a constant amplitude with constant speed, which gives us the chance to compare the code to this analytic solution.

I set up this problem (g_wave_package) in my fork in the commit 9ee6406 . (including a notebook for a quick first analysis)

An amplitude decrease is somewhat expected due to numerical dissipation, but the one found here seems a bit large. On the other hand, it looks like we can resolve the speed rather well (at least for small dt), but for some reason the wave splits and moves in both directions (which might explain the amplitude decrease).

separate the initial_models stuff into its own repo

The initial_models code is used by Maestro and Castro. It should be separated.

Test Coverage

I have been working on determining the level of coverage the current test suite has on all of the parameters that are available in MAESTRO. Currently, I have produced two lists of parameters: one that includes problem specific parameters and one that excludes problem specific parameters. These are lists of all of the parameters that aren't changed from their default values in at least one of the tests in the suite. I have attached the two lists here, these reports were run with the most recent execution of the MAESTRO test suite on 2018-04-16

Out of the parameters that aren't being covered, which are important enough that we should develop tests for them?
Coverage-NoSpecific.txt
Coverage.txt

SYSTEM_BLAS is not present anymore?

We used to allow setting SYSTEM_BLAS := t to link in the system BLAS routines instead of compiling them ourselves. This was documented in CHANGES, but it doesn't seem to exist anymore.

grid imprint in velocity field for some simulations

Some of @dwillcox 's Urca simulations should the grid structure imprint in the velocity field. It is possible that somewhere we are missing a ghost cell fill (e.g., after the projection, before plotfile, ...). It is most pronounced in the radial velocity field.

write a script to get all dependencies

We need a shell script that can clone all the dependencies for Maestro

Corner Coupling Bug

Look at IAMR development branch commit a3f9f6d096992dd63cf4733632b8afafba5083e2. There is a bug in the 3d corner coupling terms for conservative differencing. This manifested itself in IAMR test cases as a constant density field not remaining constant under divergence-free flow. This needs to be fixed in MAESTRO and MAESTROeX.

more to dimension agnostic routines

We should eliminate the different dimension-specific Cartesian routines (*_1d, *_2d, *_3d), in favor of a single _cart routine. An example of this is done in sanity.f90.

multilevel restart for wdconvect broken

inputs_3d.64.5dr.eq.dx_2levels and inputs_3d.128.5dr.eq.dx_2levels don't restart correctly, seeing differences in velocity of 10^-6 after 1 time step. single-level wdconvect restarts correctly, as does multilevel reacting_bubble