Merge branch 'development' into modify_profile

Exec/reacting_tests/reacting_convergence/ci-benchmarks/react_converge_128_true_sdc.out

@@ -21,16 +21,16 @@
  rho_e_sdc_source            -6.1984065959e+22       6.3082607681e+21
  Temp_sdc_source                             0                      0
  rho_He4_sdc_source              -543231.29383           78510.973919
- rho_C12_sdc_source              -1694.4890575           10.073635142
+ rho_C12_sdc_source              -1694.4890575           10.073635144
  rho_O16_sdc_source           -0.0098093776861       7.1506350196e-06
  rho_Fe56_sdc_source         -5.4492579275e-05       7.8513247889e-06
  pressure                     1.4155320805e+22       4.2608130858e+22
- kineng                        5.228354476e-13       1.6647276226e+18
+ kineng                       3.6057617076e-14       1.6647276226e+18
  soundspeed                       212069503.63           257404342.21
  Gamma_1                          1.5601126452           1.5885713572
- MachNumber                   6.8192135042e-18        0.0086982771596
+ MachNumber                   1.7908132003e-18        0.0086982771596
  entropy                           348439780.9           349209883.57
- magvort                      6.3108872418e-30       0.00018541051835
+ magvort                                     0       0.00018541051835
  divu                            -0.1163387912          0.55816306007
  eint_E                       4.5262357143e+16       6.9937847678e+16
  eint_e                       4.5262357143e+16       6.9937843728e+16
@@ -49,11 +49,11 @@
  z_velocity                                  0                      0
  t_sound_t_enuc               0.00023710316663           0.0195732693
  enuc                         2.9131490847e+15       4.5102586513e+17
- magvel                        1.446147223e-09           2067446.1363
- radvel                      -0.00067837194793           2067446.1363
+ magvel                       3.7977686648e-10           2067446.1363
+ radvel                      -0.00067839201193           2067446.1363
  circvel                                     0           11.820144682
- magmom                       0.00072307361144       1.6422547233e+12
- angular_momentum_x                         -0                     -0
+ magmom                       0.00018988843323       1.6422547233e+12
+ angular_momentum_x                          0                      0
  angular_momentum_y                          0                      0
- angular_momentum_z          -1.2410862734e+14       1.2410862747e+14
+ angular_momentum_z           -1.241086281e+14       1.2410862748e+14
 

Source/driver/Castro.cpp

@@ -216,6 +216,7 @@ Castro::read_params ()
     initialize_cpp_runparams();
 
     ParmParse pp("castro");
+    ParmParse ppa("amr");
 
     using namespace castro;
 
@@ -441,6 +442,13 @@ Castro::read_params ()
         }
     }
 
+    // Post-timestep regrids only make sense if we're subcycling.
+    std::string subcycling_mode;
+    ppa.query("subcycling_mode", subcycling_mode);
+    if (use_post_step_regrid == 1 && subcycling_mode == "None") {
+        amrex::Error("castro.use_post_step_regrid == 1 is not consistent with amr.subcycling_mode = None.");
+    }
+
 #ifdef AMREX_PARTICLES
     read_particle_params();
 #endif
@@ -545,7 +553,6 @@ Castro::read_params ()
 
    }
 
-   ParmParse ppa("amr");
    ppa.query("probin_file",probin_file);
 
     Vector<int> tilesize(AMREX_SPACEDIM);

Source/radiation/Make.package

@@ -31,8 +31,6 @@ CEXE_headers += filt_prim.H
 
 FEXE_headers += RAD_F.H
 
-ca_F90EXE_sources += RAD_$(DIM)D.F90 
-
 CEXE_sources += trace_ppm_rad.cpp
 
 ca_F90EXE_sources += rad_params.F90

Source/radiation/RAD_F.H

@@ -114,10 +114,6 @@ extern "C" {
                BL_FORT_FAB_ARG_3D(state),
                BL_FORT_FAB_ARG_3D(kappar));
 
-  void rfface(BL_FORT_FAB_ARG(fine), 
-              BL_FORT_FAB_ARG(crse), 
-              const int& idim, const int* irat);
-
   void FORT_RADBNDRY2(amrex::Real* bf, ARLIM_P(blo), ARLIM_P(bhi), 
                       int* tfab, ARLIM_P(dlo), ARLIM_P(dhi),
                       const amrex::Real* dx, const amrex::Real* xlo, const amrex::Real& time);

Source/radiation/Radiation.cpp

@@ -2100,17 +2100,17 @@ void Radiation::deferred_sync(int level, MultiFab& rhs, int indx)
           for (FabSetIter fsi(ref_sync_flux[lo_face]);
                fsi.isValid(); ++fsi) {
 
-            rfface(BL_TO_FORTRAN(ref_sync_flux[lo_face][fsi]),
-                   BL_TO_FORTRAN_N(crse_sync_flux[lo_face][fsi], indx),
-                   dir, ref_rat.getVect());
+            rfface(ref_sync_flux[lo_face][fsi].array(),
+                   crse_sync_flux[lo_face][fsi].array(indx),
+                   dir, ref_rat);
           }
 
           for (FabSetIter fsi(ref_sync_flux[hi_face]);
                fsi.isValid(); ++fsi) {
 
-            rfface(BL_TO_FORTRAN(ref_sync_flux[hi_face][fsi]),
-                   BL_TO_FORTRAN_N(crse_sync_flux[hi_face][fsi], indx),
-                   dir, ref_rat.getVect());
+            rfface(ref_sync_flux[hi_face][fsi].array(),
+                   crse_sync_flux[hi_face][fsi].array(indx),
+                   dir, ref_rat);
           }
         }
 

Source/radiation/rad_util.H

@@ -295,4 +295,73 @@ amrex::Real kavg(Real a, Real b, Real d, int opt)
     return k;
 }
 
+AMREX_INLINE
+void rfface (Array4<Real> const fine,
+             Array4<Real const> const crse,
+             int idim, const IntVect& irat)
+{
+    const Dim3 flo = amrex::lbound(fine);
+    const Dim3 fhi = amrex::ubound(fine);
+
+    const Dim3 clo = amrex::lbound(crse);
+    const Dim3 chi = amrex::ubound(crse);
+
+    Real ifac = 1.0_rt;
+    Real jfac = 1.0_rt;
+    Real kfac = 1.0_rt;
+    Real rfac = 1.0_rt;
+
+    if (idim == 0) {
+#if AMREX_SPACEDIM >= 2
+        jfac = irat[1];
+#endif
+#if AMREX_SPACEDIM == 3
+        kfac = irat[2];
+#endif
+
+        Real rfac = jfac * kfac;
+
+        int i = flo.x;
+        for (int k = flo.z; k <= fhi.z; ++k) {
+            for (int j = flo.y; j <= fhi.y; ++j) {
+                fine(i,j,k) = crse(clo.x, j / jfac, k / kfac) / rfac;
+            }
+        }
+    }
+    else if (idim == 1) {
+#if AMREX_SPACEDIM >= 2
+        ifac = irat[0];
+#endif
+#if AMREX_SPACEDIM == 3
+        kfac = irat[2];
+#endif
+
+        rfac = ifac * kfac;
+
+        int j = flo.y;
+        for (int k = flo.z; k <= fhi.z; ++k) {
+            for (int i = flo.x; i <= fhi.x; ++i) {
+                fine(i,j,k) = crse(i / ifac, clo.y, k / kfac) / rfac;
+            }
+        }
+    }
+    else {
+#if AMREX_SPACEDIM >= 2
+        ifac = irat[0];
+#endif
+#if AMREX_SPACEDIM == 3
+        jfac = irat[1];
+#endif
+
+        rfac = ifac * jfac;
+
+        int k = flo.z;
+        for (int j = flo.y; j <= fhi.y; ++j) {
+            for (int i = flo.x; i <= fhi.x; ++i) {
+                fine(i,j,k) = crse(i / ifac, j / jfac, clo.z) / rfac;
+            }
+        }
+    }
+}
+
 #endif

Source/sdc/Castro_sdc.cpp

@@ -208,8 +208,6 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt)
             {
                 normalize_species_sdc(i, j, k, U_center_arr);
             });
-            // ca_normalize_species(AMREX_INT_ANYD(bx1.loVect()), AMREX_INT_ANYD(bx1.hiVect()),
-            //                      BL_TO_FORTRAN_ANYD(U_center));
 
             // convert the C source to cell-centers
             C_center.resize(bx1, NUM_STATE);
@@ -234,6 +232,13 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt)
                 sdc_update_centers_o4(i, j, k, U_center_arr, U_new_center_arr, C_center_arr, dt_m, sdc_iteration);
             });
 
+            // enforce that the species sum to one after the reaction solve
+            amrex::ParallelFor(bx1,
+            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                normalize_species_sdc(i, j, k, U_new_center_arr);
+            });
+
             // compute R_i and in 1 ghost cell and then convert to <R> in
             // place (only for the interior)
             R_new.resize(bx1, NUM_STATE);
@@ -350,6 +355,13 @@ Castro::construct_old_react_source(MultiFab& U_state,
 
             make_cell_center(obx, U_state.array(mfi), U_center_arr, domain_lo, domain_hi);
 
+            // sometimes the Laplacian can make the species go negative near discontinuities
+            amrex::ParallelFor(obx,
+            [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                normalize_species_sdc(i, j, k, U_center_arr);
+            });
+
             // burn, including one ghost cell
             R_center.resize(obx, NUM_STATE);
             Elixir elix_r_center = R_center.elixir();

Source/sdc/Castro_sdc_util.H

@@ -52,21 +52,14 @@ sdc_solve(const Real dt_m,
     int ierr;
     Real err_out;
 
-    // for debugging
-    GpuArray<Real, NUM_STATE> U_orig;
-
-    for (int n = 0; n < NUM_STATE; ++n) {
-        U_orig[n] = U_old[n];
-    }
-
     if (sdc_solver == NEWTON_SOLVE) {
         // We are going to assume we already have a good guess
         // for the solve in U_new and just pass the solve onto
         // the main Newton solve
         sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_out, ierr);
 
         // failing?
-        if (ierr != NEWTON_SUCCESS) {
+        if (ierr != newton::NEWTON_SUCCESS) {
             Abort("Newton subcycling failed in sdc_solve");
         }
     } else if (sdc_solver == VODE_SOLVE) {
@@ -85,7 +78,7 @@ sdc_solve(const Real dt_m,
         sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_out, ierr);
 
         // Failing?
-        if (ierr != NEWTON_SUCCESS) {
+        if (ierr != newton::NEWTON_SUCCESS) {
             Abort("Newton failure in sdc_solve");
         }
     }