77// / \brief Defines a test problem for pressureless spherical collapse.
88// /
99#include " hydro/hydro_system.hpp"
10- #include " math/interpolate.hpp"
11- #include < limits>
1210
1311#include " AMReX_BLassert.H"
1412#include " AMReX_MultiFab.H"
1513#include " AMReX_ParmParse.H"
1614#include " QuokkaSimulation.hpp"
17- #include " util/fextract.hpp"
18- #ifdef HAVE_PYTHON
19- #include " util/matplotlibcpp.h"
20- #endif
15+ #include " util/BC.hpp"
2116
2217struct GlobalConfig {
2318 static int num_particles;
@@ -143,13 +138,6 @@ template <> void QuokkaSimulation<CollapseProblem>::ComputeDerivedVar(int lev, s
143138 }
144139}
145140
146- // A GPU helper function to set up the density floor
147- auto localDensityFloor (amrex::Real x, amrex::Real y, amrex::Real z) -> amrex::Real
148- {
149- // density_floor_expr = "0.1 * (1.0 - sqrt(x*x + y*y + z*z))"
150- return 0.1 * (1.0 - std::sqrt (x * x + y * y + z * z));
151- }
152-
153141auto problem_main () -> int
154142{
155143 amrex::ParmParse const pp (" problem" );
@@ -158,64 +146,13 @@ auto problem_main() -> int
158146
159147 // Problem initialization
160148 QuokkaSimulation<CollapseProblem> sim;
161- sim.densityFloor_ = 1.0e-5 ;
162149
163150 // initialize
164151 sim.setInitialConditions ();
165152
166- // read output variables
167- auto [position, values] = fextract (sim.state_new_cc_ [0 ], sim.Geom (0 ), 2 , 0.0 , true ); // z direction
168- const int nz = static_cast <int >(position.size ());
169-
170- // extract density and check floor
171- std::vector<double > zs (nz);
172- std::vector<double > rho_z (nz);
173- std::vector<double > custom_floor_z (nz);
174- amrex::Real min_density = std::numeric_limits<amrex::Real>::max ();
175- amrex::Real min_density_ratio = std::numeric_limits<amrex::Real>::max ();
176- for (int i = 0 ; i < nz; ++i) {
177- amrex::Real const z = position[i];
178- custom_floor_z.at (i) = localDensityFloor (0.0 , 0.0 , z); // note that the real x and y are 0.5 * delta_x
179- amrex::Real const rho = values.at (HydroSystem<CollapseProblem>::density_index)[i];
180- zs.at (i) = z;
181- rho_z.at (i) = rho;
182- min_density = std::min (min_density, rho);
183- min_density_ratio = std::min (min_density_ratio, rho / custom_floor_z.at (i));
184- }
185-
186- // Check that custom floor is working: min_density_ratio should not be smaller than 0.98
187- int status = 0 ;
188- if (min_density_ratio > 0.99 ) {
189- amrex::Print () << " Custom density floor test PASSED: min density ratio = " << min_density_ratio << " > 0.99\n " ;
190- } else {
191- amrex::Print () << " Custom density floor test FAILED: min density ratio = " << min_density_ratio << " < 0.99\n " ;
192- status = 1 ;
193- }
194-
195- #ifdef HAVE_PYTHON
196- // Plot results
197- matplotlibcpp::clf ();
198-
199- std::map<std::string, std::string> rho_args;
200- std::map<std::string, std::string> custom_floor_args;
201- rho_args[" label" ] = " rho" ;
202- rho_args[" linestyle" ] = " -" ;
203- rho_args[" color" ] = " C0" ;
204- custom_floor_args[" label" ] = " custom floor" ;
205- custom_floor_args[" linestyle" ] = " --" ;
206- custom_floor_args[" color" ] = " C1" ;
207- matplotlibcpp::plot (zs, rho_z, rho_args);
208- matplotlibcpp::plot (zs, custom_floor_z, custom_floor_args);
209- matplotlibcpp::ylim (1.0e-3 , 11.0 );
210- matplotlibcpp::yscale (" log" );
211-
212- matplotlibcpp::legend ();
213- matplotlibcpp::title (" Spherical Collapse with custom floor" );
214- matplotlibcpp::save (" ./spherical_collapse_density_floor_z.pdf" );
215- #endif // HAVE_PYTHON
216-
217153 // evolve
218154 sim.evolve ();
219155
156+ int const status = 0 ;
220157 return status;
221158}
0 commit comments