Compute emf #642

AstroKriel · 2024-06-07T05:14:46Z

Description

Added infrastructure to compute EMF components on cell edges.

Still in progress:

Close the compute-loop and update magnetic field components on each face.
Post a flow chart of steps performed in ComputeEMF().
Implement an Alfven Wave test problem.

Related issues

#526

Checklist

I have added a description (see above).
I have added a link to any related issues see (see above).
I have read the Contributing Guide.
I have added tests for any new physics that this PR adds to the code.
I have tested this PR on my local computer and all tests pass.
I have manually triggered the GPU tests with the magic comment /azp run.
I have requested a reviewer for this PR.

… woohoo)

for more information, see https://pre-commit.ci

src/RadhydroSimulation.hpp

+		auto &stateNew_cc = state_inter_cc_;
+
+		auto const &stateOld_fc = state_old_fc_tmp;
+		auto &stateNew_fc = state_inter_fc_;


src/RadhydroSimulation.hpp

 		auto const &stateInter = state_inter_cc_;
 		auto &stateFinal = state_new_cc_[lev];
-		auto [fluxArrays, faceVel] = computeHydroFluxes(stateInter, ncompHydro_, lev);
+		auto [fluxArrays, faceVel, fspds] = computeHydroFluxes(stateInter, ncompHydro_, lev);


BenWibking · 2024-06-10T17:14:47Z

src/RadhydroSimulation.hpp

+		auto const &stateOld_cc = state_old_cc_tmp;
+		auto &stateNew_cc = state_inter_cc_;
+
+		auto const &stateOld_fc = state_old_fc_tmp;


When MHD is disabled, you should call amrex::ignore_unused(stateOld_fc) to suppress this warning.

src/RadhydroSimulation.hpp

+		auto const &stateOld_fc = state_old_fc_tmp;
+		auto &stateNew_fc = state_inter_fc_;
+
+		auto [fluxArrays, faceVel, fspds] = computeHydroFluxes(stateOld_cc, ncompHydro_, lev);


src/hydro_system.hpp

 			// bx = 0 for testing purposes
 			// TODO(Neco): pass correct bx value once magnetic fields are enabled
-			F_canonical = quokka::Riemann::HLLD<problem_t, nscalars_, nmscalars_, nvar_>(sL, sR, gamma_, 0.0);
+			auto [F_canonical, fspd_m, fspd_p] = quokka::Riemann::HLLD<problem_t, nscalars_, nmscalars_, nvar_>(sL, sR, gamma_, 0.0);


BenWibking · 2024-06-07T15:47:54Z

I've marked this as a draft PR since you've indicated in the description that there are a few items you're still working on. Ping me and @markkrumholz when it's ready.

BenWibking · 2024-06-07T16:24:11Z

Just FYI- I've added #526 to "Related Issues" in the PR description.

BenWibking

Overall it's looks like a good start.

My main comment is that the variable names are excessively short and unreadable. There are also several unnecessary variable renamings in existing code in RadhydroSimulation that should (IMO) be avoided.

Just FYI, it will be necessary to save the wavespeeds for LLF when MHD is used because LLF is always the fallback solver when first-order flux correction (FOFC) is activated. Some additional changes will be necessary to make the emf update consistent when FOFC is activated.

BenWibking · 2024-06-09T17:23:56Z

src/RadhydroSimulation.hpp

@@ -1550,13 +1602,13 @@ void RadhydroSimulation<problem_t>::hydroFOFluxFunction(amrex::MultiFab const &p
 	// donor-cell reconstruction
 	HydroSystem<problem_t>::template ReconstructStatesConstant<DIR>(primVar, leftState, rightState, ng_reconstruct, nvars);
 	// LLF solver
-	HydroSystem<problem_t>::template ComputeFluxes<RiemannSolver::LLF, DIR>(flux, faceVel, leftState, rightState, primVar, artificialViscosityK_);
+	HydroSystem<problem_t>::template ComputeFluxes<RiemannSolver::LLF, DIR>(flux, faceVel, leftState, rightState, primVar, artificialViscosityK_, nullptr);


It is actually necessary to save the wavespeeds from the LLF Riemann solver -- this is used whenever first-order flux correction is activated, regardless of whether MHD is used or not.

BenWibking · 2024-06-09T17:26:26Z

src/HLLD.hpp

@@ -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];


I would strongly prefer to avoid Fortran-style four-character variable names. In my opinion, there is no reason in the year 2024 not to use descriptive variables names that use full words.

BenWibking · 2024-06-09T17:27:33Z

src/RadhydroSimulation.hpp

 	for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-		auto ba_face = amrex::convert(ba, amrex::IntVect::TheDimensionVector(idim));
+		auto ba_fc = amrex::convert(ba_cc, amrex::IntVect::TheDimensionVector(idim));


Is it necessary to rename ba_face -> ba_fc? ba_face is readable and ba_fc is shortens "face" to "fc" for no good reason.

BenWibking · 2024-06-09T17:28:34Z

src/RadhydroSimulation.hpp

-		rightState[idim] = amrex::MultiFab(ba_face, dm, nvars, reconstructGhost);
-		flux[idim] = amrex::MultiFab(ba_face, dm, nvars, 0);
-		facevel[idim] = amrex::MultiFab(ba_face, dm, 1, 0);
+		auto ba_fc = amrex::convert(ba_cc, amrex::IntVect::TheDimensionVector(idim));


Changing ba_face to ba_fc makes the variable name less understandable. There is no good reason for this change.

I had only changed it so that the variable naming was consistent across the codebase, but I am happy to change it back if you like.

Ah, okay. Making it consistent is a good reason.

BenWibking · 2024-06-09T17:29:03Z

src/hydro_system.hpp

@@ -116,7 +118,8 @@ template <typename problem_t> class HydroSystem : public HyperbolicSystem<proble

 	template <RiemannSolver RIEMANN, FluxDir DIR>
 	static void ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::MultiFab &x1FaceVel_mf, amrex::MultiFab const &x1LeftState_mf,
-				  amrex::MultiFab const &x1RightState_mf, amrex::MultiFab const &primVar_mf, amrex::Real K_visc);
+				  amrex::MultiFab const &x1RightState_mf, amrex::MultiFab const &primVar_mf, amrex::Real K_visc,
+				  amrex::MultiFab *x1FSpds_mf = nullptr);


It is necessary to save the wavespeeds for LLF when MHD is used.

BenWibking · 2024-06-09T17:30:28Z

src/mhd_system.hpp

+					 const int reconstructionOrder)
+{
+	if (reconstructionOrder == 3) {
+		switch (dir) {


This switch statement is not necessary if you make dir a template parameter.

BenWibking · 2024-06-09T17:30:41Z

src/mhd_system.hpp

+				break;
+		}
+	} else if (reconstructionOrder == 2) {
+		switch (dir) {


This switch statement is not necessary. Make dir a template parameter instead.

BenWibking · 2024-06-09T17:30:54Z

src/mhd_system.hpp

+				break;
+		}
+	} else if (reconstructionOrder == 1) {
+		switch (dir) {


This switch statement is not necessary. Make dir a template parameter instead.

BenWibking · 2024-06-09T17:35:43Z

src/mhd_system.hpp

+								iquad1 = (iface == 0) ? 3 : 2;
+							}
+							ec_fabs_Ui_q[icomp][iquad0].atomicAdd(ec_fabs_U_ieside[0], 0, 0, 1);
+							ec_fabs_Ui_q[icomp][iquad1].atomicAdd(ec_fabs_U_ieside[1], 0, 0, 1);


Why is atomicAdd necessary here?

BenWibking · 2024-06-09T17:38:26Z

src/mhd_system.hpp

+void MHDSystem<problem_t>::ComputeEMF(std::array<std::array<amrex::MultiFab, 2>, AMREX_SPACEDIM> &ec_mf_emf_comps, amrex::MultiFab const &cc_mf_cVars,
+				      std::array<amrex::MultiFab, AMREX_SPACEDIM> const &fcx_mf_cVars,
+				      std::array<amrex::MultiFab, AMREX_SPACEDIM> const &fcx_mf_fspds, const int nghost_fc)
+{


There are a lot of very small (one-line or similar) GPU kernels in this function. This will have a significant negative performance impact. Once test problems are working, I would strongly recommend seeing if it's possible to combine this kernels in this function into a much smaller number of ParallelFors.

BenWibking · 2024-06-09T22:00:50Z

src/RadhydroSimulation.hpp

+				emf_comps_ec[idim][1].define(
+				    amrex::convert(ba_cc, amrex::IntVect::TheDimensionVector(idim) + amrex::IntVect::TheDimensionVector((idim + 2) % 3)), dm, 1,
+				    nghost_fc_);
+				fspds[idim].FillBoundary(geom[lev].periodicity());


Just a note for the future: this will have to be modified for AMR to use FillPatchTwoLevels instead (since the wavespeeds may be at a coarse-fine boundary).

BenWibking · 2024-06-10T17:13:33Z

src/RadhydroSimulation.hpp

+				    nghost_fc_);
+				fspds[idim].FillBoundary(geom[lev].periodicity());
+			}
+			MHDSystem<problem_t>::ComputeEMF(emf_comps_ec, stateOld_cc, stateOld_fc, fspds, nghost_fc_);


WeightedSync should be called here (#526).

BenWibking · 2024-06-10T17:16:43Z

src/RadhydroSimulation.hpp

@@ -1212,8 +1258,8 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
 		redoFlag.setVal(quokka::redoFlag::none);



It looks like the RK stage 2 update is missing. MHDSystem<problem_t>::ComputeEMF should be called here, followed by WeightedSync.

BenWibking · 2024-06-10T17:25:31Z

src/RadhydroSimulation.hpp

+			std::array<std::array<amrex::MultiFab, 2>, AMREX_SPACEDIM> emf_comps_ec;
+			for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+				emf_comps_ec[idim][0].define(
+				    amrex::convert(ba_cc, amrex::IntVect::TheDimensionVector(idim) + amrex::IntVect::TheDimensionVector((idim + 1) % 3)), dm, 1,


I don't think this is correct. There should only be 1 emf MultiFab per spatial dimension. The emf multifabs should be identical to the E multifabs here: https://github.com/WeiqunZhang/amrex-devtests/blob/8ceb0e5e9ba61f553499bf19b21fcac21e25e7c7/edgefluxregister/main.cpp#L48

… with different permutations can live in the same multifab

…tiple boxes

… and computing fc-reference solution

…itialise additive quantities to 0). still broken: MHD fast wave values are NaN in ghost_zones.

sonarcloud · 2024-11-18T03:23:04Z

Quality Gate failed

Failed conditions
B Maintainability Rating on New Code (required ≥ A)

See analysis details on SonarQube Cloud

Catch issues before they fail your Quality Gate with our IDE extension SonarQube for IDE

AstroKriel and others added 24 commits May 8, 2024 13:44

add call to solve induction eqn

703537c

draft func to solve induction eqn

d79c692

indicate box-array centering

ae7a364

initialise fc state variables

a8fe9fd

focus on setup to reconstruct

44e1ded

debug: add evolve step to check induction eqn

db4f3d1

fc-rhs is only used for mhd

b3fdd3a

don't index comp (=1)

4ba7bc7

fix index-range of reconstruction (can compile and run w.o. errors --…

33a7b0e

… woohoo)

fix arg parsing to induction soln

4309caa

recuce domain size for terminal debugging

d2c3bd4

specify state centering

65f747a

reinstate all reconstruct options

03e6001

add alfven wave test problem (not implemented yet)

ec40a30

compute emf computes things correctly, but output still needs work

df25a22

propogate fast mhd wave speeds to compute emf func

83af98a

remove debug evolve step

6afe3fa

add alfven wave test problem (not implemented yet)

7552e65

output fast mhd wave speed

90b5417

code to save mhd wavespeeds

1cba5a5

return emf components

40231d1

make a view of the bfield data / don't make an unnecessary copy

9593208

remove unused code

17fbf88

[pre-commit.ci] auto fixes from pre-commit.com hooks

c2dc7c4

for more information, see https://pre-commit.ci

github-advanced-security bot found potential problems Jun 7, 2024

View reviewed changes

BenWibking marked this pull request as draft June 7, 2024 15:47

BenWibking reviewed Jun 9, 2024

View reviewed changes

BenWibking reviewed Jun 10, 2024

View reviewed changes

BenWibking added the MHD label Jun 20, 2024

AstroKriel added 15 commits July 9, 2024 09:00

store edge-centered emf components more compactly: data reconstructed…

bffd5b3

… with different permutations can live in the same multifab

more meaningful variable names

b81e07c

store edge-centered emf components more compactly: data reconstructed…

6bfa86e

… with different permutations can live in the same multifab

synchronise emf values that are on identical edges shared between mul…

b49a1f1

…tiple boxes

add mask for the synching emf components

23aaa22

rename IC funcs for consistency and to indicate centering

5b05c4d

compute momentum for ICs

4dedeb1

rename vars for consistency

bd79eb1

(in-progress) setup of b-fields for Riemann solver

41622e7

we have an extra cc-ghost for MHD problems

0227b20

being pedantic about syntax

9b2509b

update backend for mhd

dd70e56

check alfven wave test passes

9e7a72c

implemented induction eqn solved by line integral of emf around face,…

cd9c937

… and computing fc-reference solution

debugging: NaNs produced by MHD solver. fixed: EMF producing NaNs (in…

8fa82fa

…itialise additive quantities to 0). still broken: MHD fast wave values are NaN in ghost_zones.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Compute emf #642

Compute emf #642

AstroKriel commented Jun 7, 2024 •

edited by BenWibking

Loading

BenWibking Jun 10, 2024

BenWibking commented Jun 7, 2024

BenWibking commented Jun 7, 2024

BenWibking left a comment

BenWibking Jun 9, 2024

BenWibking Jun 9, 2024

BenWibking Jun 9, 2024

BenWibking Jun 9, 2024

AstroKriel Jun 28, 2024

BenWibking Jun 28, 2024 •

edited

Loading

BenWibking Jun 9, 2024

BenWibking Jun 9, 2024

BenWibking Jun 9, 2024

BenWibking Jun 9, 2024

BenWibking Jun 9, 2024

BenWibking Jun 9, 2024

BenWibking Jun 9, 2024

BenWibking Jun 10, 2024

BenWibking Jun 10, 2024 •

edited

Loading

BenWibking Jun 10, 2024

sonarcloud bot commented Nov 18, 2024

		@@ -1212,8 +1258,8 @@ auto RadhydroSimulation<problem_t>::advanceHydroAtLevel(amrex::MultiFab &state_o
		redoFlag.setVal(quokka::redoFlag::none);

Compute emf #642

Are you sure you want to change the base?

Compute emf #642

Conversation

AstroKriel commented Jun 7, 2024 • edited by BenWibking Loading

Description

Related issues

Checklist

Choose a reason for hiding this comment

BenWibking commented Jun 7, 2024

BenWibking commented Jun 7, 2024

BenWibking left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

BenWibking Jun 28, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

BenWibking Jun 10, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sonarcloud bot commented Nov 18, 2024

Quality Gate failed

AstroKriel commented Jun 7, 2024 •

edited by BenWibking

Loading

BenWibking Jun 28, 2024 •

edited

Loading

BenWibking Jun 10, 2024 •

edited

Loading