diff --git a/src/RadForce/test_radiation_force.cpp b/src/RadForce/test_radiation_force.cpp index f7fc0ba0a..02560b182 100644 --- a/src/RadForce/test_radiation_force.cpp +++ b/src/RadForce/test_radiation_force.cpp @@ -7,6 +7,7 @@ /// \brief Defines a test problem for radiation force terms. /// +#include #include #include "AMReX.H" @@ -38,7 +39,6 @@ constexpr double tau = 1.0e-6; // optical depth (dimensionless) constexpr double rho0 = 1.0e5 * mu; // g cm^-3 constexpr double Mach0 = 1.1; // Mach number at wind base constexpr double Mach1 = 2.128410288469465339; -constexpr double rho1 = (Mach0 / Mach1) * rho0; constexpr double Frad0 = rho0 * a0 * c_light_cgs_ / tau; // erg cm^-2 s^-1 constexpr double g0 = kappa0 * Frad0 / c_light_cgs_; // cm s^{-2} @@ -60,7 +60,7 @@ template <> struct Physics_Traits { // face-centred static constexpr bool is_mhd_enabled = false; // number of radiation groups - static constexpr int nGroups = 2; + static constexpr int nGroups = 1; }; template <> struct RadSystem_Traits { @@ -69,7 +69,6 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr double energy_unit = C::ev2erg; - static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{0., 13.6, inf}; // eV static constexpr int beta_order = 1; }; @@ -89,68 +88,18 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const -> quokka::valarray::nGroups> { quokka::valarray::nGroups> kappaFVec{}; - kappaFVec[0] = kappa0 * 1.5; - kappaFVec[1] = kappa0 * 0.5; - return kappaFVec; -} - -// declare global variables -// initial conditions read from file -amrex::Gpu::HostVector x_arr; -amrex::Gpu::HostVector rho_arr; -amrex::Gpu::HostVector Mach_arr; - -amrex::Gpu::DeviceVector x_arr_g; -amrex::Gpu::DeviceVector rho_arr_g; -amrex::Gpu::DeviceVector Mach_arr_g; - -template <> void RadhydroSimulation::preCalculateInitialConditions() -{ - std::string filename = "../extern/pressure_tube/optically_thin_wind.txt"; - std::ifstream fstream(filename, std::ios::in); - AMREX_ALWAYS_ASSERT(fstream.is_open()); - std::string header; - std::getline(fstream, header); - - for (std::string line; std::getline(fstream, line);) { - std::istringstream iss(line); - std::vector values; - for (double value = NAN; iss >> value;) { - values.push_back(value); - } - auto x = values.at(0); // position - auto rho = values.at(1); // density - auto Mach = values.at(2); // Mach number - - x_arr.push_back(x); - rho_arr.push_back(rho); - Mach_arr.push_back(Mach); + for (int g = 0; g < nGroups_; ++g) { + kappaFVec[g] = kappa0; } - - // copy to device - x_arr_g.resize(x_arr.size()); - rho_arr_g.resize(rho_arr.size()); - Mach_arr_g.resize(Mach_arr.size()); - - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, x_arr.begin(), x_arr.end(), x_arr_g.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, rho_arr.begin(), rho_arr.end(), rho_arr_g.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, Mach_arr.begin(), Mach_arr.end(), Mach_arr_g.begin()); - amrex::Gpu::streamSynchronizeAll(); + return kappaFVec; } template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) { // extract variables required from the geom object - amrex::GpuArray dx = grid_elem.dx_; - amrex::GpuArray prob_lo = grid_elem.prob_lo_; const amrex::Box &indexRange = grid_elem.indexRange_; const amrex::Array4 &state_cc = grid_elem.array_; - auto const &x_ptr = x_arr_g.dataPtr(); - auto const &rho_ptr = rho_arr_g.dataPtr(); - auto const &Mach_ptr = Mach_arr_g.dataPtr(); - int x_size = static_cast(x_arr_g.size()); - // calculate radEnergyFractions quokka::valarray::nGroups> radEnergyFractions{}; for (int g = 0; g < Physics_Traits::nGroups; ++g) { @@ -159,12 +108,7 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid(quo // loop over the grid and set the initial condition amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - amrex::Real const x = (prob_lo[0] + (i + amrex::Real(0.5)) * dx[0]) / Lx; - amrex::Real const D = interpolate_value(x, x_ptr, rho_ptr, x_size); - amrex::Real const Mach = interpolate_value(x, x_ptr, Mach_ptr, x_size); - - amrex::Real const rho = D * rho0; - amrex::Real const vel = Mach * a0; + amrex::Real const rho = rho0; for (int g = 0; g < Physics_Traits::nGroups; ++g) { state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = @@ -175,7 +119,7 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid(quo } state_cc(i, j, k, RadSystem::gasDensity_index) = rho; - state_cc(i, j, k, RadSystem::x1GasMomentum_index) = rho * vel; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = 0; state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0; state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0; state_cc(i, j, k, RadSystem::gasEnergy_index) = 0; @@ -204,7 +148,6 @@ AMRSimulation::setCustomBoundaryConditions(const amrex::IntVect &iv amrex::Box const &box = geom.Domain(); amrex::GpuArray lo = box.loVect3d(); - amrex::GpuArray hi = box.hiVect3d(); amrex::Real const Erad = Frad0 / c_light_cgs_; amrex::Real const Frad = Frad0; @@ -220,13 +163,6 @@ AMRSimulation::setCustomBoundaryConditions(const amrex::IntVect &iv // left side rho = rho0; vel = Mach0 * a0; - } else if (i > hi[0]) { - // right side - rho = rho1; - vel = Mach1 * a0; - } - - if ((i < lo[0]) || (i > hi[0])) { // Dirichlet for (int g = 0; g < Physics_Traits::nGroups; ++g) { consVar(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad * radEnergyFractions[g]; @@ -249,6 +185,7 @@ auto problem_main() -> int // Problem parameters // const int nx = 128; constexpr double CFL_number = 0.4; + double max_dt = 1.0e10; constexpr double tmax = 10.0 * (Lx / a0); constexpr int max_timesteps = 1e6; @@ -258,7 +195,7 @@ auto problem_main() -> int for (int n = 0; n < nvars; ++n) { // for x-axis: BCs_cc[n].setLo(0, amrex::BCType::ext_dir); - BCs_cc[n].setHi(0, amrex::BCType::ext_dir); + BCs_cc[n].setHi(0, amrex::BCType::foextrap); // for y-, z- axes: for (int i = 1; i < AMREX_SPACEDIM; ++i) { // periodic @@ -267,6 +204,10 @@ auto problem_main() -> int } } + // Read max_dt from parameter file + amrex::ParmParse const pp; + pp.query("max_dt", max_dt); + // Problem initialization RadhydroSimulation sim(BCs_cc); @@ -277,50 +218,74 @@ auto problem_main() -> int sim.radiationCflNumber_ = CFL_number; sim.maxTimesteps_ = max_timesteps; sim.plotfileInterval_ = -1; + sim.maxDt_ = max_dt; // initialize sim.setInitialConditions(); - auto [position0, values0] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0); // evolve sim.evolve(); // read output variables auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0); - const int nx = static_cast(position0.size()); + const int nx = static_cast(position.size()); // compute error norm std::vector xs(nx); + std::vector xs_norm(nx); std::vector rho_arr(nx); - std::vector rho_exact_arr(nx); - std::vector rho_err(nx); - std::vector vx_arr(nx); - std::vector vx_exact_arr(nx); - std::vector Frad_err(nx); + std::vector Mach_arr(nx); for (int i = 0; i < nx; ++i) { xs.at(i) = position[i]; - double rho_exact = values0.at(RadSystem::gasDensity_index)[i]; - double x1GasMom_exact = values0.at(RadSystem::x1GasMomentum_index)[i]; - double rho = values.at(RadSystem::gasDensity_index)[i]; - double Frad = values.at(RadSystem::x1RadFlux_index)[i]; - double x1GasMom = values.at(RadSystem::x1GasMomentum_index)[i]; - double vx = x1GasMom / rho; - double vx_exact = x1GasMom_exact / rho_exact; - - vx_arr.at(i) = vx / a0; - vx_exact_arr.at(i) = vx_exact / a0; - Frad_err.at(i) = (Frad - Frad0) / Frad0; - rho_err.at(i) = (rho - rho_exact) / rho_exact; - rho_exact_arr.at(i) = rho_exact; + xs_norm.at(i) = position[i] / Lx; + + double const rho = values.at(RadSystem::gasDensity_index)[i]; + double const x1GasMom = values.at(RadSystem::x1GasMomentum_index)[i]; + double const vx = x1GasMom / rho; + rho_arr.at(i) = rho; + Mach_arr.at(i) = vx / a0; } + // read in exact solution + std::vector x_exact; + std::vector x_exact_scaled; + std::vector rho_exact; + std::vector Mach_exact; + + std::string const filename = "../extern/pressure_tube/optically_thin_wind.txt"; + std::ifstream fstream(filename, std::ios::in); + AMREX_ALWAYS_ASSERT(fstream.is_open()); + std::string header; + std::getline(fstream, header); + + for (std::string line; std::getline(fstream, line);) { + std::istringstream iss(line); + std::vector values; + for (double value = NAN; iss >> value;) { + values.push_back(value); + } + auto x = values.at(0); // position + auto rho = values.at(1); // density + auto Mach = values.at(2); // Mach number + + x_exact.push_back(x); + x_exact_scaled.push_back(x * Lx); + rho_exact.push_back(rho * rho0); + Mach_exact.push_back(Mach); + } + + // interpolate exact solution to simulation grid + std::vector Mach_interp(nx); + + interpolate_arrays(xs_norm.data(), Mach_interp.data(), nx, x_exact.data(), Mach_exact.data(), static_cast(x_exact.size())); + double err_norm = 0.; double sol_norm = 0.; for (int i = 0; i < nx; ++i) { - err_norm += std::abs(rho_arr[i] - rho_exact_arr[i]); - sol_norm += std::abs(rho_exact_arr[i]); + err_norm += std::abs(Mach_arr[i] - Mach_interp[i]); + sol_norm += std::abs(Mach_interp[i]); } const double rel_err_norm = err_norm / sol_norm; @@ -333,12 +298,14 @@ auto problem_main() -> int #ifdef HAVE_PYTHON // Plot density - std::map rho_args; - std::unordered_map rhoexact_args; - rho_args["label"] = "simulation"; + std::map rhoexact_args; + std::unordered_map rho_args; rhoexact_args["label"] = "exact solution"; - matplotlibcpp::plot(xs, rho_arr, rho_args); - matplotlibcpp::scatter(xs, rho_exact_arr, 1.0, rhoexact_args); + rhoexact_args["color"] = "C0"; + rho_args["label"] = "simulation"; + rho_args["color"] = "C1"; + matplotlibcpp::plot(x_exact_scaled, rho_exact, rhoexact_args); + matplotlibcpp::scatter(xs, rho_arr, 1.0, rho_args); matplotlibcpp::legend(); matplotlibcpp::title(fmt::format("t = {:.4g} s", sim.tNew_[0])); matplotlibcpp::xlabel("x (cm)"); @@ -347,18 +314,17 @@ auto problem_main() -> int matplotlibcpp::save("./radiation_force_tube.pdf"); // plot velocity - int s = 4; // stride - std::map vx_args; - std::unordered_map vxexact_args; - vxexact_args["label"] = "exact solution"; + int const s = nx / 64; // stride + std::map vx_exact_args; + std::unordered_map vx_args; + vx_exact_args["label"] = "exact solution"; + vx_exact_args["color"] = "C0"; vx_args["label"] = "simulation"; - vx_args["color"] = "C3"; - vxexact_args["color"] = "C3"; - vxexact_args["marker"] = "o"; - // vxexact_args["edgecolors"] = "k"; + vx_args["marker"] = "o"; + vx_args["color"] = "C1"; matplotlibcpp::clf(); - matplotlibcpp::plot(xs, vx_arr, vx_args); - matplotlibcpp::scatter(strided_vector_from(xs, s), strided_vector_from(vx_exact_arr, s), 10.0, vxexact_args); + matplotlibcpp::plot(x_exact_scaled, Mach_exact, vx_exact_args); + matplotlibcpp::scatter(strided_vector_from(xs, s), strided_vector_from(Mach_arr, s), 10.0, vx_args); matplotlibcpp::legend(); // matplotlibcpp::title(fmt::format("t = {:.4g} s", sim.tNew_[0])); matplotlibcpp::xlabel("length x (cm)"); diff --git a/tests/RadForce.in b/tests/RadForce.in index 1a09d1fb3..c5cb6d09d 100644 --- a/tests/RadForce.in +++ b/tests/RadForce.in @@ -5,6 +5,9 @@ geometry.prob_lo = 0.0 0.0 0.0 geometry.prob_hi = 1.0263747986171498e16 1.0263747986171498e16 1.0263747986171498e16 geometry.is_periodic = 0 1 1 +# for convergence test, start with max_dt = 1.0e8 +max_dt = 1.0e10 + # ***************************************************************** # VERBOSITY # ***************************************************************** @@ -20,4 +23,7 @@ amr.max_grid_size_x = 128 do_reflux = 0 do_subcycle = 0 -suppress_output = 0 \ No newline at end of file +suppress_output = 1 + +cfl = 0.4 +radiation.cfl = 0.4