Add initial pertubation for cluster pgen #67

pgrete · 2023-06-02T20:41:41Z

From the docs:

Initial perturbations for both the velocity field and the magnetic field are
supported.

Note that these intial perturbation are currently incompatible with
other initial conditions that modify the velocity or magnetic field, e.g.,
an initial magnetic dipole or a uniform field.
This restriction could be lifted if required/desired but the normalization
of the fields would need to be adjusted.

In general, the perturbations will be seeded by an inverse parabolic
spectral profile centered on a peak wavenumber (in normalized units, i.e.,
k_peak = 2 mean half the box size) with a finite number of modes (default 40)
randomly chosen between k_peak/2 and 2*k_peak.

Pertubations are controlled by the following parameters:
<problem/cluster/init_perturb>
# for the velocity field
sigma_v = 0.0      # in code velocity; default: 0.0 (disabled)
k_peak_v = ???     # wavenumber in normalized units where the velocity spectrum peaks. No default value.
num_modes_v = 40   # (optional) number of wavemodes in spectral space; default: 40
sol_weight_v = 1.0 # (optional) power in solenoidal (rotational) modes of the perturbation. Range between 0 (fully compressive) and 1.0 (default, fully solenoidal).
rseed_v = 1        # (optional) integer seed for RNG for wavenumbers and amplitudes
# for the magnetic field
sigma_b = 0.0      # in code magnetic; default: 0.0 (disabled)
k_peak_b = ???     # wavenumber in normalized units where the magnetic field spectrum peaks. No default value.
num_modes_b = 40   # (optional) number of wavemodes in spectral space; default: 40
rseed_b = 2        # (optional) integer seed for RNG for wavenumbers and amplitudes

In addition, this PR also

separates the few modes Fourier transform so that it can be used for multiple pgens
allows turbulence sims with mesh refinement (as the driving is now calculated consistently across the mesh)

pgrete · 2023-06-16T16:44:32Z

src/pgen/cluster.cpp

+      a_jb.e += 1;
+      IndexRange a_kb = kb;
+      a_kb.s -= 1;
+      a_kb.e += 1;
+
+      /************************************************************
+       * Initialize an initial magnetic tower
+       ************************************************************/
+      const auto &magnetic_tower = hydro_pkg->Param<MagneticTower>("magnetic_tower");
+
+      magnetic_tower.AddInitialFieldToPotential(pmb.get(), a_kb, a_jb, a_ib, A);
+
+      /************************************************************
+       * Add dipole magnetic field to the magnetic potential
+       ************************************************************/
+      const auto &init_dipole_b_field = hydro_pkg->Param<bool>("init_dipole_b_field");
+      if (init_dipole_b_field) {
+        const Real mx = hydro_pkg->Param<Real>("dipole_b_field_mx");
+        const Real my = hydro_pkg->Param<Real>("dipole_b_field_my");
+        const Real mz = hydro_pkg->Param<Real>("dipole_b_field_mz");
+        parthenon::par_for(
+            DEFAULT_LOOP_PATTERN, "MagneticTower::AddInitialFieldToPotential",
+            parthenon::DevExecSpace(), a_kb.s, a_kb.e, a_jb.s, a_jb.e, a_ib.s, a_ib.e,
+            KOKKOS_LAMBDA(const int &k, const int &j, const int &i) {
+              // Compute and apply potential
+              const Real x = coords.Xc<1>(i);
+              const Real y = coords.Xc<2>(j);
+              const Real z = coords.Xc<3>(k);
+
+              const Real r3 = pow(SQR(x) + SQR(y) + SQR(z), 3. / 2);
+
+              const Real m_cross_r_x = my * z - mz * y;
+              const Real m_cross_r_y = mz * x - mx * z;
+              const Real m_cross_r_z = mx * y - mx * y;
+
+              A(0, k, j, i) += m_cross_r_x / (4 * M_PI * r3);
+              A(1, k, j, i) += m_cross_r_y / (4 * M_PI * r3);
+              A(2, k, j, i) += m_cross_r_z / (4 * M_PI * r3);
+            });
+      }
+
+      /************************************************************
+       * Apply the potential to the conserved variables
+       ************************************************************/
      parthenon::par_for(
-          DEFAULT_LOOP_PATTERN, "MagneticTower::AddInitialFieldToPotential",
-          parthenon::DevExecSpace(), a_kb.s, a_kb.e, a_jb.s, a_jb.e, a_ib.s, a_ib.e,
+          DEFAULT_LOOP_PATTERN, "cluster::ProblemGenerator::ApplyMagneticPotential",
+          parthenon::DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
          KOKKOS_LAMBDA(const int &k, const int &j, const int &i) {
-            // Compute and apply potential
-            const Real x = coords.Xc<1>(i);
-            const Real y = coords.Xc<2>(j);
-            const Real z = coords.Xc<3>(k);
+            u(IB1, k, j, i) =
+                (A(2, k, j + 1, i) - A(2, k, j - 1, i)) / coords.Dxc<2>(j) / 2.0 -
+                (A(1, k + 1, j, i) - A(1, k - 1, j, i)) / coords.Dxc<3>(k) / 2.0;
+            u(IB2, k, j, i) =
+                (A(0, k + 1, j, i) - A(0, k - 1, j, i)) / coords.Dxc<3>(k) / 2.0 -
+                (A(2, k, j, i + 1) - A(2, k, j, i - 1)) / coords.Dxc<1>(i) / 2.0;
+            u(IB3, k, j, i) =
+                (A(1, k, j, i + 1) - A(1, k, j, i - 1)) / coords.Dxc<1>(i) / 2.0 -
+                (A(0, k, j + 1, i) - A(0, k, j - 1, i)) / coords.Dxc<2>(j) / 2.0;
+
+            u(IEN, k, j, i) += 0.5 * (SQR(u(IB1, k, j, i)) + SQR(u(IB2, k, j, i)) +
+                                      SQR(u(IB3, k, j, i)));
+          });

-            const Real r3 = pow(SQR(x) + SQR(y) + SQR(z), 3. / 2);
+      /************************************************************
+       * Add uniform magnetic field to the conserved variables
+       ************************************************************/
+      const auto &init_uniform_b_field = hydro_pkg->Param<bool>("init_uniform_b_field");
+      if (init_uniform_b_field) {
+        const Real bx = hydro_pkg->Param<Real>("uniform_b_field_bx");
+        const Real by = hydro_pkg->Param<Real>("uniform_b_field_by");
+        const Real bz = hydro_pkg->Param<Real>("uniform_b_field_bz");
+        parthenon::par_for(
+            DEFAULT_LOOP_PATTERN, "cluster::ProblemGenerator::ApplyUniformBField",
+            parthenon::DevExecSpace(), kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+            KOKKOS_LAMBDA(const int &k, const int &j, const int &i) {
+              const Real bx_i = u(IB1, k, j, i);
+              const Real by_i = u(IB2, k, j, i);
+              const Real bz_i = u(IB3, k, j, i);
+
+              u(IB1, k, j, i) += bx;
+              u(IB2, k, j, i) += by;
+              u(IB3, k, j, i) += bz;
+
+              // Old magnetic energy is b_i^2, new Magnetic energy should be 0.5*(b_i +
+              // b)^2, add b_i*b + 0.5b^2  to old energy to accomplish that
+              u(IEN, k, j, i) +=
+                  bx_i * bx + by_i * by + bz_i * bz + 0.5 * (SQR(bx) + SQR(by) + SQR(bz));
+            });
+        // end if(init_uniform_b_field)
+      }
+
+    } // END if(hydro_pkg->Param<Fluid>("fluid") == Fluid::glmmhd)
+  }



Note that I didn't touch any of this. It's just marked because it was intended due to the explicit loop over blocks (at the ProblemGenerator is not the Mesh one and not the MeshBlock one)

pgrete · 2023-06-16T16:46:35Z

tst/regression/test_suites/cluster_tabular_cooling/cluster_tabular_cooling.py

-                prim[:, prim_col_dict["density"]], "code_mass*code_length**-3"
+                prim[prim_col_dict["density"], :], "code_mass*code_length**-3"
            )
            pres = unyt.unyt_array(
-                prim[:, prim_col_dict["pressure"]],
+                prim[prim_col_dict["pressure"], :],


bugfix from when we updated the output format in Parthenon

BenWibking · 2023-06-28T20:10:38Z

@pgrete Just wanted to check on this. I want to run some driving runs with the precipitator next week if possible.

forrestglines

Been a while since I reviewed most of this PR, I think it looks good. Good idea to test the perturbations on top of the HSE since the two don't interact at initialization.

forrestglines · 2023-06-22T07:42:40Z

docs/cluster.md

+<problem/cluster/init_perturb>
+# for the velocity field
+sigma_v = 0.0      # in code velocity; default: 0.0 (disabled)
+k_peak_v = ???     # wavenumber in normalized units where the velocity spectrum peaks. No default value.


Can we change to peak wavelength for these parameters? Or provide that as an alternative option?

It's awkward for a problem parameter to have a different action depending on the box size, especially when we'll be testing a number of box sizes.

Could also normalize k by code length instead of box length but l seems more intuitive.

I don't have a strong preference here. I guess I'm just used to normalized units.
I agree that l (in code_length) would be the most intuitive parameter given that all other params of this pgen are also in code units.
I'll add it.

forrestglines · 2023-06-22T07:44:25Z

src/pgen/cluster.cpp

+  // This could be more optimized, but require a refactor of init routines being called.
+  // However, given that it's just called during initial setup, this should not be a
+  // performance concern.
+  for (int b = 0; b < md->NumBlocks(); b++) {


I think this actually might be a performance concern. Start up is painfully slow but that might also be writing ASCENT slices and particularly profiles at the same time.

I'd be surprised if this is a big performance concern, but I'll double check.

I suspect it's Ascent or I/O. Did we ever get Darshan to work with cce/15.0.1?

I did not experience any slow startups over the past weeks when running without Ascent. So I don't think this is a concern at the moment.

src/pgen/cluster.cpp

docs/cluster.md

src/pgen/cluster.cpp

BenWibking · 2023-07-19T17:37:46Z

Just want to check- what is remaining to do on this PR?

) * Velocity and temperature ceilings * Added ceilings to velocity and temperature within AGN region * Fixed to apply only to kinetic jet injection zone * Fixed missing magnetic energy from temperature ceil * Fixed typos in ceilings * Fixed arguments for EOS constructor * cpp-py-formatter * Added V, V_a, T, D clipping within defined radius * Check kinetic feedback in release builds * Fixed Va ceiling typo * Fixed cluster/clips -> problem/cluster/clips * Update src/pgen/cluster.cpp Co-authored-by: Philipp Grete <[email protected]> * Fix comments * Fix formattin --------- Co-authored-by: Forrest Glines <[email protected]> Co-authored-by: Forrest Glines <[email protected]> Co-authored-by: par-hermes <[email protected]> Co-authored-by: Philipp Grete <[email protected]>

pgrete · 2023-08-07T15:17:51Z

The contents of this PR have been tested over the past weeks for cluster and turbulence runs and nothing surprising came up.
So in order to keep reviewing easier (I just merged the ceilings/clipping PR in here), I'm going to merge this now into main.
The remaining comments I kept track of in #46 and will address shortly (in a separate PR).

pgrete added 2 commits June 2, 2023 19:05

Move cluster pgen to MeshPgen to allow for reductions

68de864

Prepare hse test case for init pert

abc0822

pgrete changed the title ~~Add initial pertubation for cluster pgen~~ WIP Add initial pertubation for cluster pgen Jun 2, 2023

pgrete added 13 commits June 3, 2023 12:27

Git add sigma_v plus test for cluster pgen

3905cc0

Fix index error in cooling test

196afe9

Move iFT to utils from turb driver

8d05808

Isolate fmft construtor

e046a66

Separate OU and FT parts from turbulence pgen

fbbe2dc

Make FMFT smooth and consistent for mesh refinement runs

cd2420c

Fix update of internal RNG state

83de454

Remove restriction on unit domain dims

8c74b1d

Fix k_vec check on device

6ed8170

Add init perturb for velocity field

2571e1b

Move construction of modes to shared util

429d1b1

Allow ghost zones filling in FMFT

a616d93

Add initial magnetic field perturbations

639fcc1

pgrete changed the title ~~WIP Add initial pertubation for cluster pgen~~ Add initial pertubation for cluster pgen Jun 16, 2023

Add doc for cluster init perturb

b31f88e

pgrete force-pushed the pgrete/init-pert branch from cce1711 to b31f88e Compare June 16, 2023 15:27

Fix typo for using host view

3ebd164

pgrete commented Jun 16, 2023

View reviewed changes

Forrest Glines added 2 commits June 27, 2023 15:55

Fixed entropy for ASCENT output

9b114e6

Clarified feedback efficiency error

4a5010f

forrestglines approved these changes Jul 11, 2023

View reviewed changes

Merge branch 'main' into pgrete/init-pert

08ab35c

forrestglines and others added 2 commits August 7, 2023 16:08

Merge branch 'main' into pgrete/init-pert

4c51506

pgrete added 3 commits August 7, 2023 16:37

Enable first order flux correct for non-vl2 integrators

bc98afa

Fix comments

0c6596a

Fix weird missing header for non-mpi non-hdf5 builds

8057a63

pgrete enabled auto-merge (squash) August 7, 2023 15:23

pgrete merged commit ea838bf into main Aug 7, 2023

pgrete deleted the pgrete/init-pert branch March 18, 2024 10:44

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add initial pertubation for cluster pgen #67

Add initial pertubation for cluster pgen #67

pgrete commented Jun 2, 2023 •

edited

Loading

pgrete Jun 16, 2023

pgrete Jun 16, 2023

BenWibking commented Jun 28, 2023

forrestglines left a comment

forrestglines Jun 22, 2023

pgrete Jul 12, 2023

forrestglines Jun 22, 2023

pgrete Jul 12, 2023

BenWibking Jul 12, 2023

pgrete Aug 7, 2023

BenWibking commented Jul 19, 2023

pgrete commented Aug 7, 2023

Add initial pertubation for cluster pgen #67

Add initial pertubation for cluster pgen #67

Conversation

pgrete commented Jun 2, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

BenWibking commented Jun 28, 2023

forrestglines 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

Choose a reason for hiding this comment

BenWibking commented Jul 19, 2023

pgrete commented Aug 7, 2023

pgrete commented Jun 2, 2023 •

edited

Loading