Split uniform_cube and uniform_sphere into separate problems

.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