Split uniform_cube and uniform_sphere into separate problems (#2528)

.github/workflows/uniform_cube.yml

@@ -0,0 +1,39 @@
+name: uniform_cube
+
+on: [pull_request]
+jobs:
+  uniform_cube:
+    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 g++>=9.3.0
+
+      - name: Compile uniform_cube
+        run: |
+          cd Exec/gravity_tests/uniform_cube
+          make USE_MPI=FALSE -j 2
+
+      - name: Run uniform_cube
+        run: |
+          cd Exec/gravity_tests/uniform_cube
+          ./convergence.sh &> cube_convergence.out
+
+      - name: Compare to the benchmark
+        run: |
+          cd Exec/gravity_tests/uniform_cube
+          diff cube_convergence.out ci-benchmarks/cube_convergence.out

.github/workflows/uniform_sphere.yml

@@ -0,0 +1,39 @@
+name: uniform_sphere
+
+on: [pull_request]
+jobs:
+  uniform_sphere:
+    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 g++>=9.3.0
+
+      - name: Compile uniform_sphere
+        run: |
+          cd Exec/gravity_tests/uniform_sphere
+          make USE_MPI=FALSE -j 2
+
+      - name: Run uniform_sphere
+        run: |
+          cd Exec/gravity_tests/uniform_sphere
+          ./convergence.sh &> sphere_convergence.out
+
+      - name: Compare to the benchmark
+        run: |
+          cd Exec/gravity_tests/uniform_sphere
+          diff sphere_convergence.out ci-benchmarks/sphere_convergence.out

Exec/gravity_tests/uniform_cube_sphere/GNUmakefile → Exec/gravity_tests/uniform_cube/GNUmakefile

@@ -6,7 +6,7 @@ CASTRO_HOME := ../../..
 # Location of this directory. Useful if
 # you're trying to compile this from another location.
 
-TEST_DIR = $(CASTRO_HOME)/Exec/gravity_tests/uniform_cube_sphere
+TEST_DIR = $(CASTRO_HOME)/Exec/gravity_tests/uniform_cube
 
 PRECISION   ?= DOUBLE
 PROFILE     ?= FALSE

Exec/gravity_tests/uniform_cube/Prob.cpp

@@ -0,0 +1,160 @@
+#include <Castro.H>
+#include <Castro_F.H>
+
+#include <Gravity.H>
+
+#include <AMReX_ParmParse.H>
+
+#include <AMReX_buildInfo.H>
+
+#include <prob_parameters.H>
+
+#include <fundamental_constants.H>
+
+using namespace amrex;
+
+void Castro::problem_post_init()
+{
+    BL_ASSERT(level == 0);
+
+    gravity->multilevel_solve_for_new_phi(0, parent->finestLevel());
+
+    const int norm_power = 2;
+
+    ReduceOps<ReduceOpSum, ReduceOpSum> reduce_op;
+    ReduceData<Real, Real> reduce_data(reduce_op);
+    using ReduceTuple = typename decltype(reduce_data)::Type;
+
+    for (int lev = 0; lev <= parent->finestLevel(); lev++)
+    {
+        const auto dx = getLevel(lev).geom.CellSizeArray();
+        const auto problo = getLevel(lev).geom.ProbLoArray();
+
+        const Real time = getLevel(lev).state[State_Type].curTime();
+
+        auto phiGrav = getLevel(lev).derive("phiGrav", time, 0);
+
+        if (lev < parent->finestLevel())
+        {
+            const MultiFab& mask = getLevel(lev+1).build_fine_mask();
+            MultiFab::Multiply(*phiGrav, mask, 0, 0, 1, 0);
+        }
+
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+        for (MFIter mfi(*phiGrav, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+        {
+            const Box& box = mfi.tilebox();
+
+            auto phi = (*phiGrav)[mfi].array();
+            auto vol = getLevel(lev).Volume()[mfi].array();
+
+            // Compute the norm of the difference between the calculated potential
+            // and the analytical solution.
+
+            reduce_op.eval(box, reduce_data,
+            [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) -> ReduceTuple
+            {
+                Real radius = 0.5_rt * problem::diameter;
+                Real mass = (4.0_rt / 3.0_rt) * M_PI * radius * radius * radius * problem::density;
+
+                Real xx = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt) - problem::center[0];
+
+#if AMREX_SPACEDIM >= 2
+                Real yy = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5_rt) - problem::center[1];
+#else
+                Real yy = 0.0_rt;
+#endif
+
+#if AMREX_SPACEDIM == 3
+                Real zz = problo[2] + dx[2] * (static_cast<Real>(k) + 0.5_rt) - problem::center[2];
+#else
+                Real zz = 0.0_rt;
+#endif
+
+                Real rr = std::sqrt(xx * xx + yy * yy + zz * zz);
+
+                Real phiExact = 0.0_rt;
+
+                Real x[2];
+                Real y[2];
+                Real z[2];
+
+                x[0] = radius + xx;
+                x[1] = radius - xx;
+                y[0] = radius + yy;
+                y[1] = radius - yy;
+                z[0] = radius + zz;
+                z[1] = radius - zz;
+
+                for (int ii = 0; ii <= 1; ++ii) {
+                    for (int jj = 0; jj <= 1; ++jj) {
+                        for (int kk = 0; kk <= 1; ++kk) {
+
+                            Real r = std::sqrt(x[ii] * x[ii] + y[jj] * y[jj] + z[kk] * z[kk]);
+
+                            // This is Equation 20 in Katz et al. (2016). There is a special case
+                            // where, e.g., x[ii] = y[jj] = 0, in which case z[kk] = r and the
+                            // atanh will evaluate to infinity, making the product ill-defined.
+                            // Handle this case by only doing the atanh evaluation away from those
+                            // points, since the product should be zero anyway. We also want to
+                            // avoid a divide by zero, so we'll skip all the cases where r = 0.
+
+                            if (r / radius > 1.e-6_rt) {
+
+                                if (std::abs(x[ii]) / radius > 1.e-6_rt && std::abs(y[jj]) / radius > 1.e-6_rt) {
+                                    phiExact -= C::Gconst * problem::density * (x[ii] * y[jj] * std::atanh(z[kk] / r));
+                                }
+
+                                if (std::abs(y[jj]) / radius > 1.e-6_rt && std::abs(z[kk]) / radius > 1.e-6_rt) {
+                                    phiExact -= C::Gconst * problem::density * (y[jj] * z[kk] * std::atanh(x[ii] / r));
+                                }
+
+                                if (std::abs(z[kk]) / radius > 1.e-6_rt && std::abs(x[ii]) / radius > 1.e-6_rt) {
+                                    phiExact -= C::Gconst * problem::density * (z[kk] * x[ii] * std::atanh(y[jj] / r));
+                                }
+
+                                // Also, avoid a divide-by-zero for the atan terms.
+
+                                if (std::abs(x[ii]) / radius > 1.e-6_rt) {
+                                    phiExact += C::Gconst * problem::density * (x[ii] * x[ii] / 2.0_rt * std::atan(y[jj] * z[kk] / (x[ii] * r)));
+                                }
+
+                                if (std::abs(y[jj]) / radius > 1.e-6_rt) {
+                                    phiExact += C::Gconst * problem::density * (y[jj] * y[jj] / 2.0_rt * std::atan(z[kk] * x[ii] / (y[jj] * r)));
+                                }
+
+                                if (std::abs(z[kk]) / radius > 1.e-6_rt) {
+                                    phiExact += C::Gconst * problem::density * (z[kk] * z[kk] / 2.0_rt * std::atan(x[ii] * y[jj] / (z[kk] * r)));
+                                }
+
+                            }
+                        }
+                    }
+                }
+
+                Real norm_diff = vol(i,j,k) * std::pow(phi(i,j,k) - phiExact, norm_power);
+                Real norm_exact = vol(i,j,k) * std::pow(phiExact, norm_power);
+
+                return {norm_diff, norm_exact};
+            });
+        }
+    }
+
+    ReduceTuple hv = reduce_data.value();
+    Real norm_diff = amrex::get<0>(hv);
+    Real norm_exact = amrex::get<1>(hv);
+
+    ParallelDescriptor::ReduceRealSum(norm_diff);
+    ParallelDescriptor::ReduceRealSum(norm_exact);
+
+    norm_diff = std::pow(norm_diff, 1.0 / norm_power);
+    norm_exact = std::pow(norm_exact, 1.0 / norm_power);
+
+    const Real error = norm_diff / norm_exact;
+
+    amrex::Print() << std::endl;
+    amrex::Print() << "Error = " << error << std::endl;
+    amrex::Print() << std::endl;
+}

Exec/gravity_tests/uniform_cube/README

@@ -0,0 +1,3 @@
+This is a simple test of Poisson gravity. It loads a cube of uniform density
+onto the grid. The goal is to determine whether the calculated potential
+converges to the analytical potential as resolution increases.

Exec/gravity_tests/uniform_cube_sphere/_prob_params → Exec/gravity_tests/uniform_cube/_prob_params

@@ -3,5 +3,3 @@ ambient_dens       real     1.e-8_rt     y
 density            real     1.0_rt       y
 
 diameter           real     1.0_rt       y
-
-problem            integer  1            y

Exec/gravity_tests/uniform_cube/ci-benchmarks/cube_convergence.out

@@ -0,0 +1,4 @@
+ncell = 16 error = 0.004276641412
+ncell = 32 error = 0.001071555867
+ncell = 64 error = 0.0002676508753
+Average convergence rate = 1.99932628187254717105

Exec/gravity_tests/uniform_cube/convergence.sh

@@ -0,0 +1,17 @@
+#!/bin/bash
+
+EXEC=./Castro3d.gnu.ex
+
+${EXEC} inputs amr.n_cell=16 16 16 amr.plot_file=cube_plt_16_ &> cube_16.out
+error_16=$(grep "Error" cube_16.out | awk '{print $3}')
+echo "ncell = 16 error =" $error_16
+
+${EXEC} inputs amr.n_cell=32 32 32 amr.plot_file=cube_plt_32_ &> cube_32.out
+error_32=$(grep "Error" cube_32.out | awk '{print $3}')
+echo "ncell = 32 error =" $error_32
+
+${EXEC} inputs amr.n_cell=64 64 64 amr.plot_file=cube_plt_64_ &> cube_64.out
+error_64=$(grep "Error" cube_64.out | awk '{print $3}')
+echo "ncell = 64 error =" $error_64
+
+echo "Average convergence rate =" $(echo "0.5 * (sqrt($error_16 / $error_32) + sqrt($error_32 / $error_64))" | bc -l)

Exec/gravity_tests/uniform_cube/inputs

@@ -0,0 +1,63 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 0
+
+# PROBLEM SIZE & GEOMETRY
+geometry.coord_sys   =  0
+geometry.is_periodic =  0    0    0
+geometry.prob_lo     = -1.6 -1.6 -1.6
+geometry.prob_hi     =  1.6  1.6  1.6
+amr.n_cell           =  16   16   16
+
+amr.max_level        = 0
+amr.ref_ratio        = 2 2 2 2 2 2 2 2 2 2 2
+# we are not doing hydro, so there is no reflux and we don't need an error buffer
+amr.n_error_buf      = 0 0 0 0 0 0 0 0 0 0 0
+amr.blocking_factor  = 8
+amr.max_grid_size    = 8
+
+amr.refinement_indicators = denerr
+
+amr.refine.denerr.value_greater = 1.0e0
+amr.refine.denerr.field_name = density
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+
+castro.lo_bc       =  2   2   2
+castro.hi_bc       =  2   2   2
+
+# WHICH PHYSICS
+castro.do_hydro = 0
+castro.do_grav  = 1
+
+# GRAVITY
+gravity.gravity_type = PoissonGrav # Full self-gravity with the Poisson equation
+gravity.max_multipole_order = 0    # Multipole expansion includes terms up to r**(-max_multipole_order)
+gravity.rel_tol = 1.e-12           # Relative tolerance for multigrid solver
+gravity.direct_sum_bcs = 1         # Calculate boundary conditions exactly
+
+# DIAGNOSTICS & VERBOSITY
+castro.sum_interval   = 1       # timesteps between computing integrals
+amr.data_log          = grid_diag.out
+
+# CHECKPOINT FILES
+amr.checkpoint_files_output = 1
+amr.check_file        = chk      # root name of checkpoint file
+amr.check_int         = 1        # timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_files_output = 1
+amr.plot_file         = plt      # root name of plotfile
+amr.plot_per          = 1        # timesteps between plotfiles
+amr.derive_plot_vars  = ALL
+
+# PROBLEM PARAMETERS
+problem.density      = 1.0e3
+problem.diameter     = 2.0e0
+problem.ambient_dens = 1.0e-8
+
+# EOS
+eos.eos_assume_neutral = 1

Exec/gravity_tests/uniform_cube_sphere/problem_initialize.H → Exec/gravity_tests/uniform_cube/problem_initialize.H

@@ -2,10 +2,8 @@
 #define problem_initialize_H
 
 #include <prob_parameters.H>
-#include <eos.H>
 
 AMREX_INLINE
-void problem_initialize ()
-{
-}
+void problem_initialize () {}
+
 #endif

Exec/gravity_tests/uniform_cube/problem_initialize_state_data.H

@@ -0,0 +1,39 @@
+#ifndef problem_initialize_state_data_H
+#define problem_initialize_state_data_H
+
+#include <prob_parameters.H>
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void problem_initialize_state_data (int i, int j, int k, Array4<Real> const& state, const GeometryData& geomdata)
+{
+    const Real* dx = geomdata.CellSize();
+    const Real* problo = geomdata.ProbLo();
+
+    Real xx = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt) - problem::center[0];
+    Real yy = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5_rt) - problem::center[1];
+    Real zz = problo[2] + dx[2] * (static_cast<Real>(k) + 0.5_rt) - problem::center[2];
+
+    // Establish the cube
+
+    if (std::abs(xx) < problem::diameter/2 &&
+        std::abs(yy) < problem::diameter/2 &&
+        std::abs(zz) < problem::diameter/2) {
+        state(i,j,k,URHO) = problem::density;
+    }
+    else {
+        state(i,j,k,URHO) = problem::ambient_dens;
+    }
+
+    // Establish the thermodynamic quantities. They don't have to be
+    // valid because this test will never do a hydro step.
+
+    state(i,j,k,UTEMP) = 1.0_rt;
+    state(i,j,k,UEINT) = 1.0_rt;
+    state(i,j,k,UEDEN) = 1.0_rt;
+
+    for (int n = 0; n < NumSpec; n++) {
+        state(i,j,k,UFS+n) = state(i,j,k,URHO) / static_cast<Real>(NumSpec);
+    }
+}
+
+#endif