From 953ea500302737e11b5a6d29bf822da90a93ff47 Mon Sep 17 00:00:00 2001 From: Max Katz Date: Sun, 9 Jul 2023 21:51:35 -0400 Subject: [PATCH] Split uniform_cube and uniform_sphere into separate problems --- .github/workflows/uniform_cube.yml | 39 +++ .github/workflows/uniform_sphere.yml | 39 +++ .../GNUmakefile | 2 +- .../Make.package | 0 Exec/gravity_tests/uniform_cube/Prob.cpp | 160 ++++++++++ .../Problem.H | 0 Exec/gravity_tests/uniform_cube/README | 3 + .../_prob_params | 2 - .../ci-benchmarks/cube_convergence.out | 4 + .../gravity_tests/uniform_cube/convergence.sh | 17 ++ Exec/gravity_tests/uniform_cube/inputs | 63 ++++ .../problem_initialize.H | 6 +- .../problem_initialize_state_data.H | 39 +++ .../uniform_cube_sphere/Prob.cpp | 275 ------------------ Exec/gravity_tests/uniform_cube_sphere/README | 9 - Exec/gravity_tests/uniform_sphere/GNUmakefile | 35 +++ .../gravity_tests/uniform_sphere/Make.package | 2 + Exec/gravity_tests/uniform_sphere/Prob.cpp | 204 +++++++++++++ Exec/gravity_tests/uniform_sphere/Problem.H | 7 + Exec/gravity_tests/uniform_sphere/README | 3 + .../gravity_tests/uniform_sphere/_prob_params | 5 + .../ci-benchmarks/sphere_convergence.out | 4 + .../uniform_sphere/convergence.sh | 17 ++ .../inputs | 0 .../uniform_sphere/problem_initialize.H | 9 + .../problem_initialize_state_data.H | 30 +- 26 files changed, 660 insertions(+), 314 deletions(-) create mode 100644 .github/workflows/uniform_cube.yml create mode 100644 .github/workflows/uniform_sphere.yml rename Exec/gravity_tests/{uniform_cube_sphere => uniform_cube}/GNUmakefile (91%) rename Exec/gravity_tests/{uniform_cube_sphere => uniform_cube}/Make.package (100%) create mode 100644 Exec/gravity_tests/uniform_cube/Prob.cpp rename Exec/gravity_tests/{uniform_cube_sphere => uniform_cube}/Problem.H (100%) create mode 100644 Exec/gravity_tests/uniform_cube/README rename Exec/gravity_tests/{uniform_cube_sphere => uniform_cube}/_prob_params (74%) create mode 100644 Exec/gravity_tests/uniform_cube/ci-benchmarks/cube_convergence.out create mode 100755 Exec/gravity_tests/uniform_cube/convergence.sh create mode 100644 Exec/gravity_tests/uniform_cube/inputs rename Exec/gravity_tests/{uniform_cube_sphere => uniform_cube}/problem_initialize.H (69%) create mode 100644 Exec/gravity_tests/uniform_cube/problem_initialize_state_data.H delete mode 100644 Exec/gravity_tests/uniform_cube_sphere/Prob.cpp delete mode 100644 Exec/gravity_tests/uniform_cube_sphere/README create mode 100644 Exec/gravity_tests/uniform_sphere/GNUmakefile create mode 100755 Exec/gravity_tests/uniform_sphere/Make.package create mode 100644 Exec/gravity_tests/uniform_sphere/Prob.cpp create mode 100644 Exec/gravity_tests/uniform_sphere/Problem.H create mode 100644 Exec/gravity_tests/uniform_sphere/README create mode 100644 Exec/gravity_tests/uniform_sphere/_prob_params create mode 100644 Exec/gravity_tests/uniform_sphere/ci-benchmarks/sphere_convergence.out create mode 100755 Exec/gravity_tests/uniform_sphere/convergence.sh rename Exec/gravity_tests/{uniform_cube_sphere => uniform_sphere}/inputs (100%) create mode 100644 Exec/gravity_tests/uniform_sphere/problem_initialize.H rename Exec/gravity_tests/{uniform_cube_sphere => uniform_sphere}/problem_initialize_state_data.H (57%) diff --git a/.github/workflows/uniform_cube.yml b/.github/workflows/uniform_cube.yml new file mode 100644 index 0000000000..33d7e8b5df --- /dev/null +++ b/.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 diff --git a/.github/workflows/uniform_sphere.yml b/.github/workflows/uniform_sphere.yml new file mode 100644 index 0000000000..1dbd01f90b --- /dev/null +++ b/.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 diff --git a/Exec/gravity_tests/uniform_cube_sphere/GNUmakefile b/Exec/gravity_tests/uniform_cube/GNUmakefile similarity index 91% rename from Exec/gravity_tests/uniform_cube_sphere/GNUmakefile rename to Exec/gravity_tests/uniform_cube/GNUmakefile index 48b2dcb35f..e5c37c4da5 100644 --- a/Exec/gravity_tests/uniform_cube_sphere/GNUmakefile +++ b/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 diff --git a/Exec/gravity_tests/uniform_cube_sphere/Make.package b/Exec/gravity_tests/uniform_cube/Make.package similarity index 100% rename from Exec/gravity_tests/uniform_cube_sphere/Make.package rename to Exec/gravity_tests/uniform_cube/Make.package diff --git a/Exec/gravity_tests/uniform_cube/Prob.cpp b/Exec/gravity_tests/uniform_cube/Prob.cpp new file mode 100644 index 0000000000..cb7bce3401 --- /dev/null +++ b/Exec/gravity_tests/uniform_cube/Prob.cpp @@ -0,0 +1,160 @@ +#include +#include + +#include + +#include + +#include + +#include + +#include + +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 reduce_op; + ReduceData 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(i) + 0.5_rt) - problem::center[0]; + +#if AMREX_SPACEDIM >= 2 + Real yy = problo[1] + dx[1] * (static_cast(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(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; +} diff --git a/Exec/gravity_tests/uniform_cube_sphere/Problem.H b/Exec/gravity_tests/uniform_cube/Problem.H similarity index 100% rename from Exec/gravity_tests/uniform_cube_sphere/Problem.H rename to Exec/gravity_tests/uniform_cube/Problem.H diff --git a/Exec/gravity_tests/uniform_cube/README b/Exec/gravity_tests/uniform_cube/README new file mode 100644 index 0000000000..4537997ffa --- /dev/null +++ b/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. diff --git a/Exec/gravity_tests/uniform_cube_sphere/_prob_params b/Exec/gravity_tests/uniform_cube/_prob_params similarity index 74% rename from Exec/gravity_tests/uniform_cube_sphere/_prob_params rename to Exec/gravity_tests/uniform_cube/_prob_params index 95f7f71ede..6e0c70d602 100644 --- a/Exec/gravity_tests/uniform_cube_sphere/_prob_params +++ b/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 diff --git a/Exec/gravity_tests/uniform_cube/ci-benchmarks/cube_convergence.out b/Exec/gravity_tests/uniform_cube/ci-benchmarks/cube_convergence.out new file mode 100644 index 0000000000..7c839879c7 --- /dev/null +++ b/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 diff --git a/Exec/gravity_tests/uniform_cube/convergence.sh b/Exec/gravity_tests/uniform_cube/convergence.sh new file mode 100755 index 0000000000..c3ef40cb14 --- /dev/null +++ b/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) diff --git a/Exec/gravity_tests/uniform_cube/inputs b/Exec/gravity_tests/uniform_cube/inputs new file mode 100644 index 0000000000..4717fb4cc1 --- /dev/null +++ b/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 diff --git a/Exec/gravity_tests/uniform_cube_sphere/problem_initialize.H b/Exec/gravity_tests/uniform_cube/problem_initialize.H similarity index 69% rename from Exec/gravity_tests/uniform_cube_sphere/problem_initialize.H rename to Exec/gravity_tests/uniform_cube/problem_initialize.H index f6e0e77004..cabea0283d 100644 --- a/Exec/gravity_tests/uniform_cube_sphere/problem_initialize.H +++ b/Exec/gravity_tests/uniform_cube/problem_initialize.H @@ -2,10 +2,8 @@ #define problem_initialize_H #include -#include AMREX_INLINE -void problem_initialize () -{ -} +void problem_initialize () {} + #endif diff --git a/Exec/gravity_tests/uniform_cube/problem_initialize_state_data.H b/Exec/gravity_tests/uniform_cube/problem_initialize_state_data.H new file mode 100644 index 0000000000..c98d3e6456 --- /dev/null +++ b/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 + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void problem_initialize_state_data (int i, int j, int k, Array4 const& state, const GeometryData& geomdata) +{ + const Real* dx = geomdata.CellSize(); + const Real* problo = geomdata.ProbLo(); + + Real xx = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; + Real yy = problo[1] + dx[1] * (static_cast(j) + 0.5_rt) - problem::center[1]; + Real zz = problo[2] + dx[2] * (static_cast(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(NumSpec); + } +} + +#endif diff --git a/Exec/gravity_tests/uniform_cube_sphere/Prob.cpp b/Exec/gravity_tests/uniform_cube_sphere/Prob.cpp deleted file mode 100644 index b76bfb2492..0000000000 --- a/Exec/gravity_tests/uniform_cube_sphere/Prob.cpp +++ /dev/null @@ -1,275 +0,0 @@ -#include -#include - -#include - -#include - -#include - -#include - -#include - -using namespace amrex; - -void Castro::problem_post_init() -{ - BL_ASSERT(level == 0); - - // If we're doing problem 2, the normalized sphere, - // add up the mass on the domain and then update - // the density in the sphere so that it has the 'correct' - // amount of total mass. - - if (problem::problem == 2) { - - Real actual_mass = 0.0; - - bool local_flag = true; - Real time = state[State_Type].curTime(); - - for (int lev = 0; lev <= parent->finestLevel(); ++lev) - actual_mass += getLevel(lev).volWgtSum("density", time, local_flag); - - ParallelDescriptor::ReduceRealSum(actual_mass); - - // The correct amount of mass is the mass of a sphere - // with the given diameter and density. - - Real target_mass = problem::density * (1.0e0 / 6.0e0) * M_PI * std::pow(problem::diameter, 3); - - Real update_factor = target_mass / actual_mass; - - // Now update the density given this factor. - - amrex::Print() << "\n"; - amrex::Print() << " Updating density by the factor " << update_factor << " to ensure total mass matches target mass.\n"; - amrex::Print() << "\n"; - - problem::density = problem::density * update_factor; - - for (int lev = 0; lev <= parent->finestLevel(); lev++) - { - MultiFab& state = getLevel(lev).get_new_data(State_Type); - - const auto dx = getLevel(lev).geom.CellSizeArray(); - const auto problo = getLevel(lev).geom.ProbLoArray(); - -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(state, TilingIfNotGPU()); mfi.isValid(); ++mfi) { - - const Box& box = mfi.tilebox(); - - auto u = state[mfi].array(); - - // Update the density field. This ensures that the sum of the - // mass on the domain is what we intend it to be. - - amrex::ParallelFor(box, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) - { - if (problem::problem != 2) return; - - Real xx = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; - -#if AMREX_SPACEDIM >= 2 - Real yy = problo[1] + dx[1] * (static_cast(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(k) + 0.5_rt) - problem::center[2]; -#else - Real zz = 0.0_rt; -#endif - - if (std::sqrt(xx * xx + yy * yy + zz * zz) < problem::diameter / 2) { - u(i,j,k,URHO) *= update_factor; - for (int n = 0; n < NumSpec; ++n) { - u(i,j,k,UFS+n) *= update_factor; - } - } - }); - } - - } - - // Do a final check, for sanity purposes. - - actual_mass = 0.0; - - for (int lev = 0; lev <= parent->finestLevel(); lev++) - actual_mass += getLevel(lev).volWgtSum("density", time, local_flag); - - ParallelDescriptor::ReduceRealSum(actual_mass); - - if (std::abs( (actual_mass - target_mass) / target_mass ) > 1.0e-6) { - amrex::Print() << "\n"; - amrex::Print() << "Actual mass: " << actual_mass << "\n"; - amrex::Print() << "Target mass: " << target_mass << "\n"; - amrex::Print() << "\n"; - amrex::Abort("Sphere does not have the right amount of mass."); - } - - } - - gravity->multilevel_solve_for_new_phi(0, parent->finestLevel()); - - const int norm_power = 2; - - ReduceOps reduce_op; - ReduceData 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(i) + 0.5_rt) - problem::center[0]; - -#if AMREX_SPACEDIM >= 2 - Real yy = problo[1] + dx[1] * (static_cast(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(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; - - if (problem::problem == 1 || problem::problem == 2) { - - if (rr <= radius) { - phiExact = -C::Gconst * mass * (3 * radius * radius - rr * rr) / (2 * radius * radius * radius); - } else { - phiExact = -C::Gconst * mass / rr; - } - - } - else if (problem::problem == 3) { - - 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; - -} diff --git a/Exec/gravity_tests/uniform_cube_sphere/README b/Exec/gravity_tests/uniform_cube_sphere/README deleted file mode 100644 index 6a99f7e5ae..0000000000 --- a/Exec/gravity_tests/uniform_cube_sphere/README +++ /dev/null @@ -1,9 +0,0 @@ -This is a simple test of Poisson gravity. It loads a sphere (problem = 1, 2) or -cube (problem = 3) of uniform density onto the grid. The goal is to determine -whether the calculated potential converges to the analytical potential -as resolution increases. Problems 1 and 2 are the same setup except for one -detail: in Problem 2, instead of using the density requested, we modify -the density on the grid so that the total amount of mass in the domain -is equal to the mass of a sphere of the requested diameter. Problem 1 -uses the density requested by the user, and so it will not get the right -mass: the object will not be exactly spherical due to Cartesian grid effects. diff --git a/Exec/gravity_tests/uniform_sphere/GNUmakefile b/Exec/gravity_tests/uniform_sphere/GNUmakefile new file mode 100644 index 0000000000..60fb5ac5bf --- /dev/null +++ b/Exec/gravity_tests/uniform_sphere/GNUmakefile @@ -0,0 +1,35 @@ +# Define the location of the CASTRO top directory, +# if not already defined by an environment variable. + +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_sphere + +PRECISION ?= DOUBLE +PROFILE ?= FALSE + +DEBUG ?= FALSE + +DIM ?= 3 + +COMP ?= gcc + +USE_MPI ?= FALSE +USE_OMP ?= FALSE + +USE_GRAV ?= TRUE + +# This sets the EOS directory in $(MICROPHYSICS_HOME)/EOS +EOS_DIR := gamma_law + +# This sets the network directory in $(MICROPHSICS_HOME)/Networks +NETWORK_DIR ?= general_null +NETWORK_INPUTS = gammalaw.net + +Bpack += $(TEST_DIR)/Make.package +Blocs += $(TEST_DIR) + +include $(CASTRO_HOME)/Exec/Make.Castro diff --git a/Exec/gravity_tests/uniform_sphere/Make.package b/Exec/gravity_tests/uniform_sphere/Make.package new file mode 100755 index 0000000000..139597f9cb --- /dev/null +++ b/Exec/gravity_tests/uniform_sphere/Make.package @@ -0,0 +1,2 @@ + + diff --git a/Exec/gravity_tests/uniform_sphere/Prob.cpp b/Exec/gravity_tests/uniform_sphere/Prob.cpp new file mode 100644 index 0000000000..fd6cd8af4e --- /dev/null +++ b/Exec/gravity_tests/uniform_sphere/Prob.cpp @@ -0,0 +1,204 @@ +#include +#include + +#include + +#include + +#include + +#include + +#include + +using namespace amrex; + +void Castro::problem_post_init() +{ + BL_ASSERT(level == 0); + + // Add up the mass on the domain and then update the density + // in the sphere so that it has the 'correct' amount of total mass. + + Real actual_mass = 0.0; + + bool local_flag = true; + Real time = state[State_Type].curTime(); + + for (int lev = 0; lev <= parent->finestLevel(); ++lev) { + actual_mass += getLevel(lev).volWgtSum("density", time, local_flag); + } + + ParallelDescriptor::ReduceRealSum(actual_mass); + + // The correct amount of mass is the mass of a sphere + // with the given diameter and density. + + Real target_mass = problem::density * (1.0e0 / 6.0e0) * M_PI * std::pow(problem::diameter, 3); + + Real update_factor = target_mass / actual_mass; + + // Now update the density given this factor. + + amrex::Print() << "\n"; + amrex::Print() << " Updating density by the factor " << update_factor << " to ensure total mass matches target mass.\n"; + amrex::Print() << "\n"; + + problem::density = problem::density * update_factor; + + for (int lev = 0; lev <= parent->finestLevel(); lev++) + { + MultiFab& state = getLevel(lev).get_new_data(State_Type); + + const auto dx = getLevel(lev).geom.CellSizeArray(); + const auto problo = getLevel(lev).geom.ProbLoArray(); + +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi(state, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const Box& box = mfi.tilebox(); + + auto u = state[mfi].array(); + + // Update the density field. This ensures that the sum of the + // mass on the domain is what we intend it to be. + + amrex::ParallelFor(box, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + Real xx = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; + +#if AMREX_SPACEDIM >= 2 + Real yy = problo[1] + dx[1] * (static_cast(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(k) + 0.5_rt) - problem::center[2]; +#else + Real zz = 0.0_rt; +#endif + + if (std::sqrt(xx * xx + yy * yy + zz * zz) < problem::diameter / 2) { + u(i,j,k,URHO) *= update_factor; + for (int n = 0; n < NumSpec; ++n) { + u(i,j,k,UFS+n) *= update_factor; + } + } + }); + } + + } + + // Do a final check to ensure we got what we intended. + + actual_mass = 0.0; + + for (int lev = 0; lev <= parent->finestLevel(); lev++) { + actual_mass += getLevel(lev).volWgtSum("density", time, local_flag); + } + + ParallelDescriptor::ReduceRealSum(actual_mass); + + if (std::abs( (actual_mass - target_mass) / target_mass ) > 1.0e-6) { + amrex::Print() << "\n"; + amrex::Print() << "Actual mass: " << actual_mass << "\n"; + amrex::Print() << "Target mass: " << target_mass << "\n"; + amrex::Print() << "\n"; + amrex::Abort("Sphere does not have the right amount of mass."); + } + + gravity->multilevel_solve_for_new_phi(0, parent->finestLevel()); + + const int norm_power = 2; + + ReduceOps reduce_op; + ReduceData 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(i) + 0.5_rt) - problem::center[0]; + +#if AMREX_SPACEDIM >= 2 + Real yy = problo[1] + dx[1] * (static_cast(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(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; + + if (rr <= radius) { + phiExact = -C::Gconst * mass * (3 * radius * radius - rr * rr) / (2 * radius * radius * radius); + } + else { + phiExact = -C::Gconst * mass / rr; + } + + 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; +} diff --git a/Exec/gravity_tests/uniform_sphere/Problem.H b/Exec/gravity_tests/uniform_sphere/Problem.H new file mode 100644 index 0000000000..f4bc212035 --- /dev/null +++ b/Exec/gravity_tests/uniform_sphere/Problem.H @@ -0,0 +1,7 @@ +// Preprocessor directive for allowing us to do a post-initialization update. + +#ifndef DO_PROBLEM_POST_INIT +#define DO_PROBLEM_POST_INIT +#endif + +void problem_post_init(); diff --git a/Exec/gravity_tests/uniform_sphere/README b/Exec/gravity_tests/uniform_sphere/README new file mode 100644 index 0000000000..566885de8e --- /dev/null +++ b/Exec/gravity_tests/uniform_sphere/README @@ -0,0 +1,3 @@ +This is a simple test of Poisson gravity. It loads a sphere of uniform density +onto the grid. The goal is to determine whether the calculated potential +converges to the analytical potential as resolution increases. diff --git a/Exec/gravity_tests/uniform_sphere/_prob_params b/Exec/gravity_tests/uniform_sphere/_prob_params new file mode 100644 index 0000000000..6e0c70d602 --- /dev/null +++ b/Exec/gravity_tests/uniform_sphere/_prob_params @@ -0,0 +1,5 @@ +ambient_dens real 1.e-8_rt y + +density real 1.0_rt y + +diameter real 1.0_rt y diff --git a/Exec/gravity_tests/uniform_sphere/ci-benchmarks/sphere_convergence.out b/Exec/gravity_tests/uniform_sphere/ci-benchmarks/sphere_convergence.out new file mode 100644 index 0000000000..585bf5294c --- /dev/null +++ b/Exec/gravity_tests/uniform_sphere/ci-benchmarks/sphere_convergence.out @@ -0,0 +1,4 @@ +ncell = 16 error = 0.05274082586 +ncell = 32 error = 0.008316536152 +ncell = 64 error = 0.001270740323 +Average convergence rate = 2.53825936367487233285 diff --git a/Exec/gravity_tests/uniform_sphere/convergence.sh b/Exec/gravity_tests/uniform_sphere/convergence.sh new file mode 100755 index 0000000000..c3ef40cb14 --- /dev/null +++ b/Exec/gravity_tests/uniform_sphere/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) diff --git a/Exec/gravity_tests/uniform_cube_sphere/inputs b/Exec/gravity_tests/uniform_sphere/inputs similarity index 100% rename from Exec/gravity_tests/uniform_cube_sphere/inputs rename to Exec/gravity_tests/uniform_sphere/inputs diff --git a/Exec/gravity_tests/uniform_sphere/problem_initialize.H b/Exec/gravity_tests/uniform_sphere/problem_initialize.H new file mode 100644 index 0000000000..cabea0283d --- /dev/null +++ b/Exec/gravity_tests/uniform_sphere/problem_initialize.H @@ -0,0 +1,9 @@ +#ifndef problem_initialize_H +#define problem_initialize_H + +#include + +AMREX_INLINE +void problem_initialize () {} + +#endif diff --git a/Exec/gravity_tests/uniform_cube_sphere/problem_initialize_state_data.H b/Exec/gravity_tests/uniform_sphere/problem_initialize_state_data.H similarity index 57% rename from Exec/gravity_tests/uniform_cube_sphere/problem_initialize_state_data.H rename to Exec/gravity_tests/uniform_sphere/problem_initialize_state_data.H index 81e475f838..9ae21249f6 100644 --- a/Exec/gravity_tests/uniform_cube_sphere/problem_initialize_state_data.H +++ b/Exec/gravity_tests/uniform_sphere/problem_initialize_state_data.H @@ -14,30 +14,13 @@ void problem_initialize_state_data (int i, int j, int k, Array4 const& sta Real yy = problo[1] + dx[1] * (static_cast(j) + 0.5_rt) - problem::center[1]; Real zz = problo[2] + dx[2] * (static_cast(k) + 0.5_rt) - problem::center[2]; - // Establish the cube or sphere + // Establish the sphere - if (problem::problem == 1 || problem::problem == 2) { - - if (std::pow(xx * xx + yy * yy + zz * zz, 0.5_rt) < problem::diameter / 2) { - state(i,j,k,URHO) = problem::density; - } else { - state(i,j,k,URHO) = problem::ambient_dens; - } - - } else if (problem::problem == 3) { - - 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; - } - -#ifndef AMREX_USE_GPU - } else { - amrex::Error("Problem not defined."); -#endif + if (std::pow(xx * xx + yy * yy + zz * zz, 0.5_rt) < 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 @@ -51,4 +34,5 @@ void problem_initialize_state_data (int i, int j, int k, Array4 const& sta state(i,j,k,UFS+n) = state(i,j,k,URHO) / static_cast(NumSpec); } } + #endif