fix 4th order diffusion convergence + diffusion_test fixes (#2592)

In the conversion from Fortran to C++, the order of and cell-centered T diffusive flux was swapped.
AMReX-Astro · Oct 5, 2023 · af29792 · af29792
1 parent 104a57d
commit af29792
Show file tree

Hide file tree

Showing 3 changed files with 32 additions and 26 deletions.
diff --git a/Exec/unit_tests/diffusion_test/Prob.cpp b/Exec/unit_tests/diffusion_test/Prob.cpp
@@ -21,8 +21,8 @@ void Castro::problem_post_simulation(Vector<std::unique_ptr<AmrLevel> >& amr_lev
     Castro& castro = dynamic_cast<Castro&>(*amr_level[n]);
     Real time = castro.get_state_data(State_Type).curTime();
 
-    const int* domain_lo = castro.geom.Domain().loVect();
-    const int* domain_hi = castro.geom.Domain().hiVect();
+    auto domain_lo = castro.geom.Domain().loVect3d();
+    auto domain_hi = castro.geom.Domain().hiVect3d();
 
     // the state data
     MultiFab& S = castro.get_new_data(State_Type);
@@ -33,12 +33,18 @@ void Castro::problem_post_simulation(Vector<std::unique_ptr<AmrLevel> >& amr_lev
 #ifdef TRUE_SDC
     // if we are fourth-order, we need to convert to averages
     if (sdc_order == 4) {
-      for (MFIter mfi(*analytic); mfi.isValid(); ++mfi) {
+        FArrayBox tmp;
 
-        const Box& gbx = mfi.growntilebox(1);
-        ca_make_fourth_in_place(AMREX_INT_ANYD(gbx.loVect()), AMREX_INT_ANYD(gbx.hiVect()),
-                                BL_TO_FORTRAN_FAB((*analytic)[mfi]),
-                                AMREX_INT_ANYD(domain_lo), AMREX_INT_ANYD(domain_hi));
+        for (MFIter mfi(*analytic); mfi.isValid(); ++mfi) {
+
+            const Box& bx = mfi.tilebox();
+            tmp.resize(bx, 1);
+            auto tmp_arr = tmp.array();
+
+            castro.make_fourth_in_place(bx,
+                                        analytic->array(mfi),
+                                        tmp_arr,
+                                        domain_lo, domain_hi);
 
       }
     }

diff --git a/Exec/unit_tests/diffusion_test/README.md b/Exec/unit_tests/diffusion_test/README.md
@@ -69,10 +69,10 @@ analytic solution, giving:
 A convergence test of the 4th-order SDC algorithm can be run as:
 
 ```
-./Castro1d.gnu.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=64
-./Castro1d.gnu.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=128
-./Castro1d.gnu.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=256
-./Castro1d.gnu.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=512
+./Castro1d.gnu.TRUESDC.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=64
+./Castro1d.gnu.TRUESDC.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=128
+./Castro1d.gnu.TRUESDC.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=256
+./Castro1d.gnu.TRUESDC.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=512
 ```
 
 Note: this is Cartesian, not spherical (since we don't have spherical

diff --git a/Source/hydro/fourth_order.cpp b/Source/hydro/fourth_order.cpp
@@ -883,40 +883,40 @@ Castro::fourth_add_diffusive_flux(const Box& bx,
     if (idir == 0) {
 
       if (is_avg) {
-        // we are working with the cell-center state
-        dTdx = (-q_arr(i+1,j,k,temp_comp) + 27*q_arr(i,j,k,temp_comp) -
-                27*q_arr(i-1,j,k,temp_comp) + q_arr(i-2,j,k,temp_comp)) / (24.0_rt * dx[0]);
-
-      } else {
         // we are working with the cell-average state
         dTdx = (-q_arr(i+1,j,k,temp_comp) + 15*q_arr(i,j,k,temp_comp) -
                 15*q_arr(i-1,j,k,temp_comp) + q_arr(i-2,j,k,temp_comp)) / (12.0_rt * dx[0]);
+
+      } else {
+        // we are working with the cell-center state
+        dTdx = (-q_arr(i+1,j,k,temp_comp) + 27*q_arr(i,j,k,temp_comp) -
+                27*q_arr(i-1,j,k,temp_comp) + q_arr(i-2,j,k,temp_comp)) / (24.0_rt * dx[0]);
       }
 
     } else if (idir == 1) {
 
       if (is_avg) {
-        // we are working with the cell-center state
-        dTdx = (-q_arr(i,j+1,k,temp_comp) + 27*q_arr(i,j,k,temp_comp) -
-                27*q_arr(i,j-1,k,temp_comp) + q_arr(i,j-2,k,temp_comp)) / (24.0_rt * dx[1]);
-
-      } else {
         // we are working with the cell-average state
         dTdx = (-q_arr(i,j+1,k,temp_comp) + 15*q_arr(i,j,k,temp_comp) -
                 15*q_arr(i,j-1,k,temp_comp) + q_arr(i,j-2,k,temp_comp)) / (12.0_rt * dx[1]);
+
+      } else {
+        // we are working with the cell-center state
+        dTdx = (-q_arr(i,j+1,k,temp_comp) + 27*q_arr(i,j,k,temp_comp) -
+                27*q_arr(i,j-1,k,temp_comp) + q_arr(i,j-2,k,temp_comp)) / (24.0_rt * dx[1]);
       }
 
     } else {
 
       if (is_avg) {
-        // we are working with the cell-center state
-        dTdx = (-q_arr(i,j,k+1,temp_comp) + 27*q_arr(i,j,k,temp_comp) -
-                27*q_arr(i,j,k-1,temp_comp) + q_arr(i,j,k-2,temp_comp)) / (24.0_rt * dx[2]);
-
-      } else {
         // we are working with the cell-average state
         dTdx = (-q_arr(i,j,k+1,temp_comp) + 15*q_arr(i,j,k,temp_comp) -
                 15*q_arr(i,j,k-1,temp_comp) + q_arr(i,j,k-2,temp_comp)) / (12.0_rt * dx[2]);
+
+      } else {
+        // we are working with the cell-center state
+        dTdx = (-q_arr(i,j,k+1,temp_comp) + 27*q_arr(i,j,k,temp_comp) -
+                27*q_arr(i,j,k-1,temp_comp) + q_arr(i,j,k-2,temp_comp)) / (24.0_rt * dx[2]);
       }
 
     }