diff --git a/.codespell-ignore-words b/.codespell-ignore-words index fdcff634ec..5a235e1f40 100644 --- a/.codespell-ignore-words +++ b/.codespell-ignore-words @@ -11,3 +11,4 @@ thi nd ue bion +aas diff --git a/.github/workflows/c-linter.yml b/.github/workflows/c-linter.yml index 6327a1f53a..602e92f48f 100644 --- a/.github/workflows/c-linter.yml +++ b/.github/workflows/c-linter.yml @@ -24,6 +24,16 @@ jobs: sudo apt-get update -y -qq sudo apt-get -qq -y install curl clang-tidy cmake jq clang cppcheck clang-format bear g++>=9.3.0 gfortran>=9.3.0 + - name: Install hypre + run: | + wget -q https://github.com/hypre-space/hypre/archive/refs/tags/v2.28.0.tar.gz + tar xfz v2.28.0.tar.gz + cd hypre-2.28.0/src + ./configure --with-cxxstandard=17 --without-MPI + make -j 2 + make install + cd ../../ + - name: Get cpp linter repo run: | cd external @@ -38,6 +48,7 @@ jobs: - name: Run cpp linter run: | + export AMREX_HYPRE_HOME=${PWD}/hypre-2.28.0/src/hypre python3 external/cpp-linter-action/run_on_changed_files.py ${{ github.event.pull_request.base.sha }} ${{ github.event.pull_request.head.sha }} \ -ignore-files="amrex|Microphysics" \ -config-file="${GITHUB_WORKSPACE}/.clang-tidy" \ diff --git a/.github/workflows/gpu_action.yml b/.github/workflows/gpu_action.yml index 1c4d0eb531..10843ef256 100644 --- a/.github/workflows/gpu_action.yml +++ b/.github/workflows/gpu_action.yml @@ -31,7 +31,17 @@ jobs: sudo apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/3bf863cc.pub sudo add-apt-repository "deb https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/ /" sudo apt-get update - sudo apt-get -y install cuda-toolkit + sudo apt-get -y install cuda-toolkit-11-8 + + - name: Install hypre + run: | + wget -q https://github.com/hypre-space/hypre/archive/refs/tags/v2.28.0.tar.gz + tar xfz v2.28.0.tar.gz + cd hypre-2.28.0/src + CUDA_HOME=/usr/local/cuda HYPRE_CUDA_SM=60 ./configure --with-cxxstandard=17 --with-cuda --enable-unified-memory --without-MPI + make -j 2 + make install + cd ../../ - name: Get cpp linter repo run: | @@ -48,4 +58,5 @@ jobs: - name: Compile problems for GPU run: | export PATH=$PATH:/usr/local/cuda/bin + export AMREX_HYPRE_HOME=${PWD}/hypre-2.28.0/src/hypre python3 external/cpp-linter-action/run_on_changed_files.py ${{ github.event.pull_request.base.sha }} ${{ github.event.pull_request.head.sha }} -header-filter=Castro -ignore-files="amrex|Microphysics" -gpu diff --git a/.github/workflows/mhd-compare.yml b/.github/workflows/mhd-compare.yml new file mode 100644 index 0000000000..8df15f2574 --- /dev/null +++ b/.github/workflows/mhd-compare.yml @@ -0,0 +1,45 @@ +name: OrszagTang + +on: [pull_request] +jobs: + OrszagTang-3d: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + with: + fetch-depth: 0 + + - name: Get submodules + run: | + git submodule update --init + cd external/Microphysics + git fetch; git checkout development + cd ../amrex + git fetch; git checkout development + cd ../.. + + - name: Install dependencies + run: | + sudo apt-get update -y -qq + sudo apt-get -qq -y install curl cmake jq clang g++>=9.3.0 libopenmpi-dev openmpi-bin + + - name: Compile OrszagTang + run: | + cd Exec/mhd_tests/OrszagTang + make USE_MPI=TRUE -j 2 + + - name: Run OrszagTang-3d + run: | + cd Exec/mhd_tests/OrszagTang + mpirun -np 2 ./Castro3d.gnu.MPI.ex inputs.test amr.plot_files_output=1 + + - name: Build the fextrema tool + run: | + cd external/amrex/Tools/Plotfile + make programs=fextrema -j 2 + + - name: Check the extrema + run: | + cd Exec/mhd_tests/OrszagTang + ../../../external/amrex/Tools/Plotfile/fextrema.gnu.ex plt00010 > fextrema.out + diff fextrema.out ci-benchmarks/OrszagTang-3d.out diff --git a/CHANGES.md b/CHANGES.md index 6186a3f336..7db160e12f 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,5 +1,12 @@ # 23.07 + * The parameter castro.state_nghost, which allowed State_Type to have ghost + zones, has been removed. (#2502) + + * The additional ghost zone in State_Type, used when radiation is enabled, + has been removed. The checkpoint version number has been updated to avoid + restarting from a checkpoint with the wrong number of ghost zones. (#2495) + * The parameter gravity.no_composite was removed (#2483) * The parameter spherical_star was removed (#2482) diff --git a/Docs/source/rotation.rst b/Docs/source/rotation.rst index 0368ee9ff1..2bb60ea4a4 100644 --- a/Docs/source/rotation.rst +++ b/Docs/source/rotation.rst @@ -224,7 +224,7 @@ a rotating frame The kinetic energy equation can be obtained from :eq:`eq:v-rot` by -mulitplying by :math:`\rho\vbt`: +multiplying by :math:`\rho\vbt`: .. math:: \begin{align} diff --git a/Docs/source/software.rst b/Docs/source/software.rst index 0828e497a7..34d71d51a5 100644 --- a/Docs/source/software.rst +++ b/Docs/source/software.rst @@ -256,10 +256,7 @@ The current ``StateData`` names Castro carries are: simply be advected, but we will allow rotation (in particular, the Coriolis force) to affect them. - ``State_Type`` ``MultiFab`` s have no ghost cells by default for - pure hydro and a single ghost cell by default when ``RADIATION`` - is enabled. There is an option to force them to have ghost cells by - setting the parameter ``castro.state_nghost`` at runtime. + ``State_Type`` ``MultiFab`` s have no ghost cells. Note that the prediction of the hydrodynamic state to the interface will require 4 ghost cells. This accommodated by creating a separate diff --git a/Exec/mhd_tests/OrszagTang/ci-benchmarks/OrszagTang-3d.out b/Exec/mhd_tests/OrszagTang/ci-benchmarks/OrszagTang-3d.out new file mode 100644 index 0000000000..d936b2e915 --- /dev/null +++ b/Exec/mhd_tests/OrszagTang/ci-benchmarks/OrszagTang-3d.out @@ -0,0 +1,19 @@ + plotfile = plt00010 + time = 0.00017983829222995601 + variables minimum value maximum value + density 0.22099460191 0.22100305489 + xmom -0.2210626652 0.2210626652 + ymom -0.2211744105 0.2211744105 + zmom -1.3394504115e-20 1.2796480966e-20 + rho_E 0.19805685456 0.46260143162 + rho_e 0.19794567403 0.19795835472 + Temp 7.2178022365e-09 7.2179887773e-09 + rho_X 0.22099460191 0.22100305489 + pressure 0.1326236016 0.13263209766 + x_velocity -1.0002839593 1.0002839593 + y_velocity -1.0007881112 1.0007881112 + z_velocity -6.0608661064e-20 5.7902777063e-20 + B_x -0.28206060038 0.28206060038 + B_y -0.28220293547 0.28220293547 + B_z -1.5570588641e-18 1.630216695e-18 + diff --git a/Exec/reacting_tests/reacting_bubble/inputs_2d_test b/Exec/reacting_tests/reacting_bubble/inputs_2d_test index 3005709b0f..09dcdee812 100644 --- a/Exec/reacting_tests/reacting_bubble/inputs_2d_test +++ b/Exec/reacting_tests/reacting_bubble/inputs_2d_test @@ -23,10 +23,6 @@ castro.yr_ext_bc_type = 1 castro.hse_interp_temp = 1 - -# this is here just for test suite coverage -castro.state_nghost = 1 - # WHICH PHYSICS castro.do_hydro = 1 castro.do_react = 1 diff --git a/Exec/reacting_tests/reacting_bubble/inputs_3d_test b/Exec/reacting_tests/reacting_bubble/inputs_3d_test index abce4000f6..d3d3c40020 100644 --- a/Exec/reacting_tests/reacting_bubble/inputs_3d_test +++ b/Exec/reacting_tests/reacting_bubble/inputs_3d_test @@ -27,9 +27,6 @@ castro.fill_ambient_bc = 1 castro.ambient_fill_dir = 2 castro.ambient_outflow_vel = 1 -# this is here just for test suite coverage -castro.state_nghost = 1 - # WHICH PHYSICS castro.do_hydro = 1 castro.do_react = 1 diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H index 184dd46842..1f862bc8c1 100644 --- a/Source/driver/Castro.H +++ b/Source/driver/Castro.H @@ -429,16 +429,8 @@ public: /// @param time the current simulation time /// @param dt the timestep to advance (e.g., go from time to /// time + dt) -/// @param amr_iteration where we are in the current AMR subcycle. Each -/// level will take a number of steps to reach the -/// final time of the coarser level below it. This -/// counter starts at 1 -/// @param amr_ncycle the number of subcycles at this level /// - advance_status do_advance_ctu(amrex::Real time, - amrex::Real dt, - int amr_iteration, - int amr_ncycle); + advance_status do_advance_ctu (amrex::Real time, amrex::Real dt); #ifndef AMREX_USE_GPU diff --git a/Source/driver/Castro_advance_ctu.cpp b/Source/driver/Castro_advance_ctu.cpp index 9e98e9b18c..e6eeb2fb34 100644 --- a/Source/driver/Castro_advance_ctu.cpp +++ b/Source/driver/Castro_advance_ctu.cpp @@ -13,15 +13,8 @@ using namespace amrex; advance_status -Castro::do_advance_ctu(Real time, - Real dt, - int amr_iteration, - int amr_ncycle) +Castro::do_advance_ctu (Real time, Real dt) { - - amrex::ignore_unused(amr_iteration); - amrex::ignore_unused(amr_ncycle); - // this routine will advance the old state data (called Sborder here) // to the new time, for a single level. The new data is called // S_new here. The update includes reactions (if we are not doing @@ -36,21 +29,6 @@ Castro::do_advance_ctu(Real time, const Real prev_time = state[State_Type].prevTime(); const Real cur_time = state[State_Type].curTime(); - MultiFab& S_new = get_new_data(State_Type); - - MultiFab& old_source = get_old_data(Source_Type); - MultiFab& new_source = get_new_data(Source_Type); - -#ifdef MHD - MultiFab& Bx_old = get_old_data(Mag_Type_x); - MultiFab& By_old = get_old_data(Mag_Type_y); - MultiFab& Bz_old = get_old_data(Mag_Type_z); - - MultiFab& Bx_new = get_new_data(Mag_Type_x); - MultiFab& By_new = get_new_data(Mag_Type_y); - MultiFab& Bz_new = get_new_data(Mag_Type_z); -#endif - // Perform initialization steps. status = initialize_do_advance(time, dt); @@ -59,98 +37,74 @@ Castro::do_advance_ctu(Real time, return status; } - // If we are Strang splitting the reactions, do the old-time contribution now + // Perform all pre-advance operations and then initialize + // the new-time state with the output of those operators. -#ifdef REACTIONS - status = do_old_reactions(time, dt); + status = pre_advance_operators(prev_time, dt); if (status.success == false) { return status; } -#endif - - // Initialize the new-time data. This copy needs to come after the - // reactions. - - MultiFab::Copy(S_new, Sborder, 0, 0, NUM_STATE, S_new.nGrow()); - // Construct the old-time sources from Sborder. This will already - // be applied to S_new (with full dt weighting), to be correctly - // later. Note -- this does not affect the prediction of the - // interface state, an explicit source will be traced there as + // Construct the old-time sources from Sborder. This will already + // be applied to S_new (with full dt weighting), to be corrected + // later. Note that this does not affect the prediction of the + // interface state; an explicit source will be traced there as // needed. -#ifdef GRAVITY - construct_old_gravity(amr_iteration, amr_ncycle, prev_time); -#endif + status = do_old_sources(prev_time, dt); - do_old_sources( -#ifdef MHD - Bx_old, By_old, Bz_old, -#endif - old_source, Sborder, S_new, prev_time, dt); + if (status.success == false) { + return status; + } -#ifdef SIMPLIFIED_SDC -#ifdef REACTIONS - // the SDC reactive source ghost cells on coarse levels might not - // be in sync due to any average down done, so fill them here + // Perform any operations that occur after the sources but before the hydro. - MultiFab& react_src = get_new_data(Simplified_SDC_React_Type); + status = pre_hydro_operators(prev_time, dt); - AmrLevel::FillPatch(*this, react_src, react_src.nGrow(), cur_time, Simplified_SDC_React_Type, 0, react_src.nComp()); -#endif -#endif + if (status.success == false) { + return status; + } - // Do the hydro update. We build directly off of Sborder, which - // is the state that has already seen the burn + // Do the hydro update. We build directly off of Sborder, which + // is the state that has already seen the burn. - if (do_hydro) - { #ifndef MHD - status = construct_ctu_hydro_source(time, dt); + status = construct_ctu_hydro_source(prev_time, dt); #else - status = construct_ctu_mhd_source(time, dt); + status = construct_ctu_mhd_source(prev_time, dt); #endif - if (status.success == false) { - return status; - } + if (status.success == false) { + return status; } - // Construct and apply new-time source terms. + // Perform any operations that occur after the hydro but before + // the corrector sources. -#ifdef GRAVITY - construct_new_gravity(amr_iteration, amr_ncycle, cur_time); -#endif + status = post_hydro_operators(cur_time, dt); - do_new_sources( -#ifdef MHD - Bx_new, By_new, Bz_new, -#endif - new_source, Sborder, S_new, cur_time, dt); + if (status.success == false) { + return status; + } - // If the state has ghost zones, sync them up now - // since the hydro source only works on the valid zones. + // Construct and apply new-time source terms. - if (S_new.nGrow() > 0) { - clean_state( -#ifdef MHD - Bx_new, By_new, Bz_new, -#endif - S_new, cur_time, 0); + status = do_new_sources(cur_time, dt); - expand_state(S_new, cur_time, S_new.nGrow()); + if (status.success == false) { + return status; } // Do the second half of the reactions for Strang, or the full burn for simplified SDC. -#ifdef REACTIONS - status = do_new_reactions(cur_time, dt); + status = post_advance_operators(cur_time, dt); if (status.success == false) { return status; } -#endif // REACTIONS + + // Perform finalization steps. status = finalize_do_advance(cur_time, dt); @@ -252,7 +206,6 @@ Castro::retry_advance_ctu(Real dt, const advance_status& status) } return do_retry; - } @@ -404,7 +357,7 @@ Castro::subcycle_advance_ctu(const Real time, const Real dt, int amr_iteration, // We do the hydro advance here, and record whether we completed it. - status = do_advance_ctu(subcycle_time, dt_subcycle, amr_iteration, amr_ncycle); + status = do_advance_ctu(subcycle_time, dt_subcycle); if (in_retry) { in_retry = false; diff --git a/Source/driver/Castro_advance_sdc.cpp b/Source/driver/Castro_advance_sdc.cpp index 2e5a501cc3..54ad501e54 100644 --- a/Source/driver/Castro_advance_sdc.cpp +++ b/Source/driver/Castro_advance_sdc.cpp @@ -91,7 +91,7 @@ Castro::do_advance_sdc (Real time, // TODO: this is not using the density at the current stage #ifdef GRAVITY - construct_old_gravity(amr_iteration, amr_ncycle, prev_time); + construct_old_gravity(prev_time); #endif if (apply_sources()) { diff --git a/Source/driver/Castro_io.cpp b/Source/driver/Castro_io.cpp index 39bb65019b..f6cad5cb5c 100644 --- a/Source/driver/Castro_io.cpp +++ b/Source/driver/Castro_io.cpp @@ -55,11 +55,12 @@ using namespace amrex; // 9: Rotation_Type was removed from Castro // 10: Reactions_Type was removed from checkpoints // 11: PhiRot_Type was removed from Castro +// 12: State_Type's additional ghost zone, used when radiation is enabled, has been removed namespace { int input_version = -1; - int current_version = 11; + int current_version = 12; } // I/O routines for Castro diff --git a/Source/driver/Castro_setup.cpp b/Source/driver/Castro_setup.cpp index 3a8f2f5d07..c17bd8937b 100644 --- a/Source/driver/Castro_setup.cpp +++ b/Source/driver/Castro_setup.cpp @@ -348,14 +348,7 @@ Castro::variableSetUp () bool state_data_extrap = false; bool store_in_checkpoint; -#if defined(RADIATION) - // Radiation should always have at least one ghost zone. - int ngrow_state = std::max(1, state_nghost); -#else - int ngrow_state = state_nghost; -#endif - - BL_ASSERT(ngrow_state >= 0); + int ngrow_state = 0; store_in_checkpoint = true; desc_lst.addDescriptor(State_Type,IndexType::TheCellType(), diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 03ffc689a2..e6e06797ec 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -19,11 +19,6 @@ state_interp_order int 1 # 2: preserve linear combinations and prevent new extrema lin_limit_state_interp int 0 -# Number of ghost zones for state data to have. Note that -# if you are using radiation, choosing this to be zero will -# be overridden since radiation needs at least one ghost zone. -state_nghost int 0 - # do we do the hyperbolic reflux at coarse-fine interfaces? do_reflux int 1 diff --git a/Source/gravity/Castro_gravity.H b/Source/gravity/Castro_gravity.H index 3d16552ea0..b9caa2a1bd 100644 --- a/Source/gravity/Castro_gravity.H +++ b/Source/gravity/Castro_gravity.H @@ -1,23 +1,17 @@ /// /// Construct gravitational field at old timestep /// -/// @param amr_iteration where we are in the current AMR subcycle -/// @param amr_ncycle the number of subcycles at this level /// @param time current time /// - void construct_old_gravity(int amr_iteration, int amr_ncycle, - amrex::Real time); + void construct_old_gravity (amrex::Real time); /// /// Construct gravitational field at new timestep /// -/// @param amr_iteration where we are in the current AMR subcycle -/// @param amr_ncycle the number of subcycles at this level /// @param time current time /// - void construct_new_gravity(int amr_iteration, int amr_ncycle, - amrex::Real time); + void construct_new_gravity (amrex::Real time); /// /// Construct gravitational source terms at old timestep diff --git a/Source/gravity/Castro_gravity.cpp b/Source/gravity/Castro_gravity.cpp index 8f0149f50b..f0cecf9f44 100644 --- a/Source/gravity/Castro_gravity.cpp +++ b/Source/gravity/Castro_gravity.cpp @@ -11,14 +11,12 @@ using namespace amrex; void -Castro::construct_old_gravity(int amr_iteration, int amr_ncycle, Real time) +Castro::construct_old_gravity (Real time) { - - amrex::ignore_unused(amr_iteration); - amrex::ignore_unused(amr_ncycle); - BL_PROFILE("Castro::construct_old_gravity()"); + const Real strt_time = ParallelDescriptor::second(); + MultiFab& grav_old = get_old_data(Gravity_Type); MultiFab& phi_old = get_old_data(PhiGrav_Type); @@ -114,17 +112,30 @@ Castro::construct_old_gravity(int amr_iteration, int amr_ncycle, Real time) gravity->get_old_grav_vector(level, grav_old, time); + if (verbose > 0) + { + const int IOProc = ParallelDescriptor::IOProcessorNumber(); + Real run_time = ParallelDescriptor::second() - strt_time; + +#ifdef BL_LAZY + Lazy::QueueReduction( [=] () mutable { +#endif + ParallelDescriptor::ReduceRealMax(run_time,IOProc); + + amrex::Print() << "Castro::construct_old_gravity() time = " << run_time << "\n" << "\n"; +#ifdef BL_LAZY + }); +#endif + } } void -Castro::construct_new_gravity(int amr_iteration, int amr_ncycle, Real time) +Castro::construct_new_gravity (Real time) { - - amrex::ignore_unused(amr_iteration); - amrex::ignore_unused(amr_ncycle); - BL_PROFILE("Castro::construct_new_gravity()"); + const Real strt_time = ParallelDescriptor::second(); + MultiFab& grav_new = get_new_data(Gravity_Type); MultiFab& phi_new = get_new_data(PhiGrav_Type); @@ -241,6 +252,21 @@ Castro::construct_new_gravity(int amr_iteration, int amr_ncycle, Real time) } + if (verbose > 0) + { + const int IOProc = ParallelDescriptor::IOProcessorNumber(); + Real run_time = ParallelDescriptor::second() - strt_time; + +#ifdef BL_LAZY + Lazy::QueueReduction( [=] () mutable { +#endif + ParallelDescriptor::ReduceRealMax(run_time,IOProc); + + amrex::Print() << "Castro::construct_new_gravity() time = " << run_time << "\n" << "\n"; +#ifdef BL_LAZY + }); +#endif + } } void Castro::construct_old_gravity_source(MultiFab& source, MultiFab& state_in, Real time, Real dt) @@ -393,7 +419,6 @@ void Castro::construct_old_gravity_source(MultiFab& source, MultiFab& state_in, }); #endif } - } void Castro::construct_new_gravity_source(MultiFab& source, MultiFab& state_old, MultiFab& state_new, diff --git a/Source/gravity/Gravity.cpp b/Source/gravity/Gravity.cpp index 4f7b5ddd7c..cb7d1c1c95 100644 --- a/Source/gravity/Gravity.cpp +++ b/Source/gravity/Gravity.cpp @@ -705,6 +705,8 @@ Gravity::multilevel_solve_for_new_phi (int level, int finest_level_in) amrex::Print() << "... multilevel solve for new phi at base level " << level << " to finest level " << finest_level_in << std::endl; } + const Real strt = ParallelDescriptor::second(); + for (int lev = level; lev <= finest_level_in; lev++) { BL_ASSERT(grad_phi_curr[lev].size()==AMREX_SPACEDIM); for (int n=0; nget_state_data(Rad_Type).prevTime(); MultiFab& S_new = castro->get_new_data(State_Type); - AmrLevel::FillPatch(*castro,S_new,ngrow,time,State_Type,0,S_new.nComp(),0); + FillPatchIterator fpi_new(*castro, S_new, ngrow, time, State_Type, 0, S_new.nComp()); + MultiFab& S_new_border = fpi_new.get_mf(); Array lambda; if (radiation::limiter > 0) { @@ -54,11 +55,9 @@ void Radiation::MGFLD_implicit_update(int level, int iteration, int ncycle) Er_lag.setBndry(-1.0); Er_lag.FillBoundary(parent->Geom(level).periodicity()); - MultiFab& S_lag = castro->get_old_data(State_Type); - for (FillPatchIterator fpi(*castro,S_lag,ngrow,oldtime,State_Type, - 0,S_lag.nComp()); fpi.isValid(); ++fpi) { - S_lag[fpi].copy(fpi()); - } + MultiFab& S_old = castro->get_old_data(State_Type); + FillPatchIterator fpi_old(*castro, S_old, ngrow, oldtime, State_Type, 0, S_old.nComp()); + MultiFab& S_lag = fpi_old.get_mf(); MultiFab kpr_lag(grids,dmap,nGroups,1); MGFLD_compute_rosseland(kpr_lag, S_lag); @@ -111,18 +110,16 @@ void Radiation::MGFLD_implicit_update(int level, int iteration, int ncycle) #ifdef _OPENMP #pragma omp parallel #endif - for (MFIter mfi(S_new,true); mfi.isValid(); ++mfi) { - + for (MFIter mfi(S_new_border, true); mfi.isValid(); ++mfi) { const Box &gbx = mfi.growntilebox(1); const Box &bx = mfi.tilebox(); - rho[mfi].copy(S_new[mfi],gbx, URHO, gbx,0,1); + rho[mfi].copy(S_new_border[mfi], gbx, URHO, gbx, 0, 1); - rhoe_new[mfi].copy(S_new[mfi], bx, UEINT, bx,0,1); + rhoe_new[mfi].copy(S_new_border[mfi], bx, UEINT, bx, 0, 1); rhoe_old[mfi].copy(rhoe_new[mfi], bx); - temp_new[mfi].copy(S_new[mfi],gbx, UTEMP, gbx,0,1); - + temp_new[mfi].copy(S_new_border[mfi], gbx, UTEMP, gbx, 0, 1); } // Planck mean and Rosseland @@ -206,7 +203,7 @@ void Radiation::MGFLD_implicit_update(int level, int iteration, int ncycle) it++; if (it == 1) { - eos_opacity_emissivity(S_new, temp_new, + eos_opacity_emissivity(S_new_border, temp_new, temp_star, // input kappa_p, kappa_r, jg, djdT, dkdT, dedT, // output @@ -404,9 +401,9 @@ void Radiation::MGFLD_implicit_update(int level, int iteration, int ncycle) etaT, etaTz, eta1, coupT, kappa_p, jg, mugT, - S_new, level, delta_t, ptc_tau, it, conservative_update); + S_new_border, level, delta_t, ptc_tau, it, conservative_update); - eos_opacity_emissivity(S_new, temp_new, + eos_opacity_emissivity(S_new_border, temp_new, temp_star, // input kappa_p, kappa_r, jg, djdT, dkdT, dedT, // output @@ -474,9 +471,9 @@ void Radiation::MGFLD_implicit_update(int level, int iteration, int ncycle) if (!converged && it > n_bisect) { bisect_matter(rhoe_new, temp_new, rhoe_star, temp_star, - S_new, grids, level); + S_new_border, grids, level); - eos_opacity_emissivity(S_new, temp_new, + eos_opacity_emissivity(S_new_border, temp_new, temp_star, // input kappa_p, kappa_r, jg, djdT, dkdT, dedT, // output diff --git a/Source/sdc/Castro_sdc_util.H b/Source/sdc/Castro_sdc_util.H index 435554a411..7f84100c57 100644 --- a/Source/sdc/Castro_sdc_util.H +++ b/Source/sdc/Castro_sdc_util.H @@ -541,7 +541,7 @@ sdc_vode_solve(const Real dt_m, } C_react[NumSpec+1] = C[UEINT]; - dvode_t dvode_state; + dvode_t dvode_state; burn_t burn_state; diff --git a/Source/sdc/vode_rhs_true_sdc.H b/Source/sdc/vode_rhs_true_sdc.H index c491b46e33..5a1c18e86c 100644 --- a/Source/sdc/vode_rhs_true_sdc.H +++ b/Source/sdc/vode_rhs_true_sdc.H @@ -17,8 +17,9 @@ // This is a generic interface that calls the specific RHS routine in the // network you're actually using. +template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void rhs(const Real time, burn_t& burn_state, dvode_t& vode_state, RArray1D& dUdt, const bool in_jacobian=false) +void rhs(const Real time, burn_t& burn_state, DvodeT& vode_state, RArray1D& dUdt, const bool in_jacobian=false) { // note: dUdt is 1-based @@ -74,9 +75,9 @@ void rhs(const Real time, burn_t& burn_state, dvode_t& vode_state, RArray1D& dUd } -template +template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void jac (const Real time, burn_t& burn_state, dvode_t& vode_state, MatrixType& pd) +void jac (const Real time, burn_t& burn_state, DvodeT& vode_state, MatrixType& pd) { // Call the specific network routine to get the Jacobian. diff --git a/Source/sources/Castro_sources.H b/Source/sources/Castro_sources.H index b4ac46288c..8dc34db4b3 100644 --- a/Source/sources/Castro_sources.H +++ b/Source/sources/Castro_sources.H @@ -47,6 +47,16 @@ amrex::Real time, amrex::Real dt, bool apply_to_state = true); +/// +/// Construct source terms at old time (convenience wrapper for CTU advance) +/// +/// @param time the current simulation time +/// @param dt the timestep to advance (e.g., go from time to +/// time + dt) +/// + advance_status do_old_sources (amrex::Real time, amrex::Real dt); + + /// /// Construct the old-time source /// @@ -83,6 +93,16 @@ amrex::Real time, amrex::Real dt, bool apply_to_state = true); +/// +/// Construct new time sources (convenience wrapper for CTU advance) +/// +/// @param time the current simulation time +/// @param dt the timestep to advance (e.g., go from time to +/// time + dt) +/// + advance_status do_new_sources (amrex::Real time, amrex::Real dt); + + /// /// Construct the new-time sources. /// @@ -316,3 +336,45 @@ /// void fill_geom_source(amrex::Real time, amrex::Real dt, amrex::MultiFab& cons_state, amrex::MultiFab& geom_src); + + +/// +/// Perform all operations that occur prior to computing the predictor sources +/// and the hydro advance. +/// +/// @param time current time +/// @param dt timestep +/// + + advance_status pre_advance_operators (amrex::Real time, amrex::Real dt); + +/// +/// Perform all operations that occur after computing the predictor sources +/// but before the hydro advance. +/// +/// @param time current time +/// @param dt timestep +/// + + advance_status pre_hydro_operators (amrex::Real time, amrex::Real dt); + + +/// +/// Perform all operations that occur after the hydro source +/// but before the corrector sources. +/// +/// @param time current time +/// @param dt timestep +/// + + advance_status post_hydro_operators (amrex::Real time, amrex::Real dt); + + +/// +/// Perform all operations that occur after the corrector sources. +/// +/// @param time current time +/// @param dt timestep +/// + + advance_status post_advance_operators (amrex::Real time, amrex::Real dt); diff --git a/Source/sources/Castro_sources.cpp b/Source/sources/Castro_sources.cpp index 86bc343d1e..ab31fb5513 100644 --- a/Source/sources/Castro_sources.cpp +++ b/Source/sources/Castro_sources.cpp @@ -172,7 +172,30 @@ Castro::do_old_sources( }); #endif } +} + +advance_status +Castro::do_old_sources (Real time, Real dt) +{ + advance_status status {}; + + MultiFab& S_new = get_new_data(State_Type); + + MultiFab& old_source = get_old_data(Source_Type); + +#ifdef MHD + MultiFab& Bx_old = get_old_data(Mag_Type_x); + MultiFab& By_old = get_old_data(Mag_Type_y); + MultiFab& Bz_old = get_old_data(Mag_Type_z); +#endif + + do_old_sources( +#ifdef MHD + Bx_old, By_old, Bz_old, +#endif + old_source, Sborder, S_new, time, dt); + return status; } void @@ -235,6 +258,30 @@ Castro::do_new_sources( } +advance_status +Castro::do_new_sources (Real time, Real dt) +{ + advance_status status {}; + + MultiFab& S_new = get_new_data(State_Type); + + MultiFab& new_source = get_new_data(Source_Type); + +#ifdef MHD + MultiFab& Bx_new = get_new_data(Mag_Type_x); + MultiFab& By_new = get_new_data(Mag_Type_y); + MultiFab& Bz_new = get_new_data(Mag_Type_z); +#endif + + do_new_sources( +#ifdef MHD + Bx_new, By_new, Bz_new, +#endif + new_source, Sborder, S_new, time, dt); + + return status; +} + void Castro::construct_old_source(int src, MultiFab& source, MultiFab& state_in, Real time, Real dt) { @@ -471,3 +518,104 @@ Castro::print_all_source_changes(Real dt, bool is_new) evaluate_and_print_source_change(source, dt, source_name); } + +// Perform all operations that occur prior to computing the predictor sources +// and the hydro advance. + +advance_status +Castro::pre_advance_operators (Real time, Real dt) +{ + amrex::ignore_unused(time); + amrex::ignore_unused(dt); + + advance_status status {}; + + // If we are Strang splitting the reactions, do the old-time contribution now. + +#ifndef TRUE_SDC +#ifdef REACTIONS + status = do_old_reactions(time, dt); + + if (status.success == false) { + return status; + } +#endif +#endif + + // If we are using gravity, solve for the potential and gravatational field. + +#ifdef GRAVITY + construct_old_gravity(time); +#endif + + // Initialize the new-time data. This copy needs to come after all Strang-split operators. + + MultiFab& S_new = get_new_data(State_Type); + + MultiFab::Copy(S_new, Sborder, 0, 0, NUM_STATE, S_new.nGrow()); + + return status; +} + +// Perform all operations that occur after computing the predictor sources +// but before the hydro advance. + +advance_status +Castro::pre_hydro_operators (Real time, Real dt) +{ + amrex::ignore_unused(time); + amrex::ignore_unused(dt); + + advance_status status {}; + +#ifdef SIMPLIFIED_SDC +#ifdef REACTIONS + // The SDC reactive source ghost cells on coarse levels might not + // be in sync due to any average down done, so fill them here. + + MultiFab& react_src = get_new_data(Simplified_SDC_React_Type); + + AmrLevel::FillPatch(*this, react_src, react_src.nGrow(), time + dt, Simplified_SDC_React_Type, 0, react_src.nComp()); +#endif +#endif + + return status; +} + +// Perform all operations that occur after the hydro source +// but before the corrector sources. + +advance_status +Castro::post_hydro_operators (Real time, Real dt) +{ + amrex::ignore_unused(time); + amrex::ignore_unused(dt); + + advance_status status {}; + +#ifdef GRAVITY + construct_new_gravity(time); +#endif + + return status; +} + +// Perform all operations that occur after the corrector sources. + +advance_status +Castro::post_advance_operators (Real time, Real dt) +{ + advance_status status {}; + +#ifndef TRUE_SDC +#ifdef REACTIONS + status = do_new_reactions(time, dt); + + if (status.success == false) { + return status; + } +#endif +#endif + + return status; +}