Compute emf

src/AlfvenWave/CMakeLists.txt

@@ -0,0 +1,7 @@
+add_executable(test_alfven_wave test_alfven_wave.cpp ../interpolate.cpp ${QuokkaObjSources})
+
+if(AMReX_GPU_BACKEND MATCHES "CUDA")
+    setup_target_for_cuda_compilation(test_alfven_wave)
+endif()
+
+add_test(NAME AlfvenWave COMMAND test_alfven_wave alfven_wave.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)

src/AlfvenWave/test_alfven_wave.cpp

@@ -0,0 +1,174 @@
+//==============================================================================
+// Copyright 2022 Neco Kriel.
+// Released under the MIT license. See LICENSE file included in the GitHub repo.
+//==============================================================================
+/// \file test_fc_quantities.cpp
+/// \brief Defines a test problem to make sure face-centred quantities are created correctly.
+///
+
+#include <cassert>
+#include <ostream>
+#include <stdexcept>
+#include <valarray>
+
+#include "AMReX_Array.H"
+#include "AMReX_Array4.H"
+#include "AMReX_Print.H"
+#include "AMReX_REAL.H"
+
+#include "RadhydroSimulation.hpp"
+#include "fextract.hpp"
+#include "grid.hpp"
+#include "physics_info.hpp"
+#include "test_alfven_wave.hpp"
+
+struct AlfvenWave {
+};
+
+template <> struct quokka::EOS_Traits<AlfvenWave> {
+	static constexpr double gamma = 5. / 3.;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
+};
+
+template <> struct Physics_Traits<AlfvenWave> {
+	// cell-centred
+	static constexpr bool is_hydro_enabled = true;
+	static constexpr int numMassScalars = 0;		     // number of mass scalars
+	static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars
+	static constexpr bool is_radiation_enabled = false;
+	// face-centred
+	static constexpr bool is_mhd_enabled = true;
+	static constexpr int nGroups = 1; // number of radiation groups
+};
+
+constexpr double rho0 = 1.0;					   // background density
+constexpr double P0 = 1.0 / quokka::EOS_Traits<AlfvenWave>::gamma; // background pressure
+constexpr double v0 = 0.;					   // background velocity
+constexpr double amp = 1.0e-6;					   // perturbation amplitude
+
+AMREX_GPU_DEVICE void computeWaveSolution(int i, int j, int k, amrex::Array4<amrex::Real> const &state, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx,
+					  amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_lo)
+{
+	const amrex::Real x_L = prob_lo[0] + (i + amrex::Real(0.0)) * dx[0];
+	const amrex::Real x_R = prob_lo[0] + (i + amrex::Real(1.0)) * dx[0];
+	const amrex::Real A = amp;
+
+	const quokka::valarray<double, 3> R = {1.0, -1.0, 1.5}; // right eigenvector of sound wave
+	const quokka::valarray<double, 3> U_0 = {rho0, rho0 * v0, P0 / (quokka::EOS_Traits<AlfvenWave>::gamma - 1.0) + 0.5 * rho0 * std::pow(v0, 2)};
+	const quokka::valarray<double, 3> dU = (A * R / (2.0 * M_PI * dx[0])) * (std::cos(2.0 * M_PI * x_L) - std::cos(2.0 * M_PI * x_R));
+
+	double rho = U_0[0] + dU[0];
+	double xmom = U_0[1] + dU[1];
+	double Etot = U_0[2] + dU[2];
+	double Eint = Etot - 0.5 * (xmom * xmom) / rho;
+
+	state(i, j, k, HydroSystem<AlfvenWave>::density_index) = rho;
+	state(i, j, k, HydroSystem<AlfvenWave>::x1Momentum_index) = xmom;
+	state(i, j, k, HydroSystem<AlfvenWave>::x2Momentum_index) = 0;
+	state(i, j, k, HydroSystem<AlfvenWave>::x3Momentum_index) = 0;
+	state(i, j, k, HydroSystem<AlfvenWave>::energy_index) = Etot;
+	state(i, j, k, HydroSystem<AlfvenWave>::internalEnergy_index) = Eint;
+}
+
+template <> void RadhydroSimulation<AlfvenWave>::setInitialConditionsOnGrid(quokka::grid grid_elem)
+{
+	// extract grid information
+	amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx = grid_elem.dx_;
+	amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> prob_lo = grid_elem.prob_lo_;
+	const amrex::Array4<double> &state = grid_elem.array_;
+	const amrex::Box &indexRange = grid_elem.indexRange_;
+
+	const int ncomp_cc = Physics_Indices<AlfvenWave>::nvarTotal_cc;
+	// loop over the grid and set the initial condition
+	amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
+		for (int n = 0; n < ncomp_cc; ++n) {
+			state(i, j, k, n) = 0; // fill unused quantities with zeros
+		}
+		computeWaveSolution(i, j, k, state, dx, prob_lo);
+	});
+}
+
+template <> void RadhydroSimulation<AlfvenWave>::setInitialConditionsOnGridFaceVars(quokka::grid grid_elem)
+{
+	// extract grid information
+	const amrex::Array4<double> &state = grid_elem.array_;
+	const amrex::Box &indexRange = grid_elem.indexRange_;
+	const quokka::direction dir = grid_elem.dir_;
+
+	if (dir == quokka::direction::x) {
+		amrex::ParallelFor(indexRange,
+				   [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { state(i, j, k, MHDSystem<AlfvenWave>::bfield_index) = 1.0 + (i % 2); });
+	} else if (dir == quokka::direction::y) {
+		amrex::ParallelFor(indexRange,
+				   [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { state(i, j, k, MHDSystem<AlfvenWave>::bfield_index) = 2.0 + (j % 2); });
+	} else if (dir == quokka::direction::z) {
+		amrex::ParallelFor(indexRange,
+				   [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { state(i, j, k, MHDSystem<AlfvenWave>::bfield_index) = 3.0 + (k % 2); });
+	}
+}
+
+void checkMFs(amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> const &state1,
+	      amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> const &state2)
+{
+	double err = 0.0;
+	for (int level = 0; level < state1.size(); ++level) {
+		for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+			// initialise MF
+			const amrex::BoxArray &ba = state1[level][idim].boxArray();
+			const amrex::DistributionMapping &dm = state1[level][idim].DistributionMap();
+			int ncomp = state1[level][idim].nComp();
+			int ngrow = state1[level][idim].nGrow();
+			amrex::MultiFab mf_diff(ba, dm, ncomp, ngrow);
+			// compute difference between two MFs (at level)
+			amrex::MultiFab::Copy(mf_diff, state1[level][idim], 0, 0, ncomp, ngrow);
+			amrex::MultiFab::Subtract(mf_diff, state2[level][idim], 0, 0, ncomp, ngrow);
+			// compute error (summed over each component)
+			for (int icomp = 0; icomp < Physics_Indices<AlfvenWave>::nvarPerDim_fc; ++icomp) {
+				err += mf_diff.norm1(icomp);
+			}
+		}
+	}
+	amrex::Print() << "Accumulated error in MFs read from chk-file: " << err << "\n";
+	amrex::Print() << "\n";
+	AMREX_ALWAYS_ASSERT(std::abs(err) == 0.0);
+}
+
+auto problem_main() -> int
+{
+	// Problem initialization
+	const int ncomp_cc = Physics_Indices<AlfvenWave>::nvarTotal_cc;
+	amrex::Vector<amrex::BCRec> BCs_cc(ncomp_cc);
+	for (int n = 0; n < ncomp_cc; ++n) {
+		for (int i = 0; i < AMREX_SPACEDIM; ++i) {
+			BCs_cc[n].setLo(i, amrex::BCType::int_dir); // periodic
+			BCs_cc[n].setHi(i, amrex::BCType::int_dir);
+		}
+	}
+
+	const int nvars_fc = Physics_Indices<AlfvenWave>::nvarTotal_fc;
+	amrex::Vector<amrex::BCRec> BCs_fc(nvars_fc);
+	for (int n = 0; n < nvars_fc; ++n) {
+		for (int i = 0; i < AMREX_SPACEDIM; ++i) {
+			BCs_fc[n].setLo(i, amrex::BCType::int_dir); // periodic
+			BCs_fc[n].setHi(i, amrex::BCType::int_dir);
+		}
+	}
+
+	RadhydroSimulation<AlfvenWave> sim_write(BCs_cc, BCs_fc);
+	sim_write.setInitialConditions();
+	sim_write.evolve();
+	amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> const &state_new_fc_write = sim_write.getNewMF_fc();
+	amrex::Print() << "\n";
+
+	RadhydroSimulation<AlfvenWave> sim_restart(BCs_cc, BCs_fc);
+	sim_restart.setChkFile("chk00000");
+	sim_restart.setInitialConditions();
+	amrex::Vector<amrex::Array<amrex::MultiFab, AMREX_SPACEDIM>> const &state_new_fc_restart = sim_restart.getNewMF_fc();
+	amrex::Print() << "\n";
+
+	amrex::Print() << "Checking new FC MFs...\n";
+	checkMFs(state_new_fc_write, state_new_fc_restart);
+
+	return 0;
+}

src/AlfvenWave/test_alfven_wave.hpp

@@ -0,0 +1,20 @@
+#ifndef TEST_ALFVEN_WAVE_HPP_ // NOLINT
+#define TEST_ALFVEN_WAVE_HPP_
+//==============================================================================
+// Copyright 2024 Neco Kriel.
+// Released under the MIT license. See LICENSE file included in the GitHub repo.
+//==============================================================================
+/// \file test_alfven_wave.hpp
+/// \brief Defines a test problem to check alfven waves propogate correctly.
+///
+
+// external headers
+#include <fmt/format.h>
+
+// internal headers
+#include "hydro_system.hpp"
+#include "mhd_system.hpp"
+
+// function definitions
+
+#endif // TEST_ALFVEN_WAVE_HPP_

src/CMakeLists.txt

@@ -138,6 +138,9 @@ add_subdirectory(Advection)
 add_subdirectory(Advection2D)
 add_subdirectory(AdvectionSemiellipse)
 
+add_subdirectory(FCQuantities)
+add_subdirectory(AlfvenWave)
+
 add_subdirectory(HydroBlast2D)
 add_subdirectory(HydroBlast3D)
 add_subdirectory(HydroContact)
@@ -180,7 +183,6 @@ add_subdirectory(RadhydroPulseMG)
 
 add_subdirectory(BinaryOrbitCIC)
 add_subdirectory(Cooling)
-add_subdirectory(FCQuantities)
 add_subdirectory(NSCBC)
 add_subdirectory(ODEIntegration)
 add_subdirectory(PassiveScalar)

src/HLLD.hpp

@@ -42,7 +42,7 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto FastMagnetoSonicSpeed(double gamma, quo
 // HLLD solver following Miyoshi and Kusano (2005), hereafter MK5.
 template <typename problem_t, int N_scalars, int N_mscalars, int fluxdim>
 AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_mscalars> const &sL, quokka::HydroState<N_scalars, N_mscalars> const &sR,
-					      const double gamma, const double bx) -> quokka::valarray<double, fluxdim>
+					      const double gamma, const double bx) -> std::tuple<quokka::valarray<double, fluxdim>, double, double>
 {
 	//--- Step 1. Compute L/R states
 
@@ -96,6 +96,8 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
 	const double cfs_R = FastMagnetoSonicSpeed(gamma, sR, bx);
 	spds[0] = std::min(sL.u - cfs_L, sR.u - cfs_R);
 	spds[4] = std::max(sL.u + cfs_L, sR.u + cfs_R);
+	const double fspd_m = spds[0];
+	const double fspd_p = spds[0];
 
 	//--- Step 3. Compute L/R fluxes
 
@@ -328,7 +330,7 @@ AMREX_FORCE_INLINE AMREX_GPU_DEVICE auto HLLD(quokka::HydroState<N_scalars, N_ms
 
 	// TODO(neco): Eint=0 for now; pscalars will also be needed in the future.
 	quokka::valarray<double, fluxdim> F_hydro = {f_x.rho, f_x.mx, f_x.my, f_x.mz, f_x.E, 0.0};
-	return F_hydro;
+	return std::make_tuple(std::move(F_hydro), fspd_m, fspd_p);
 }
 } // namespace quokka::Riemann