3D HW

examples/H3LAPD/3D-hw/cuboid_periodic_5x5x10.geo

@@ -0,0 +1,37 @@
+//=============================== Parameters ==================================
+// Lengths and resolutions in each dimension
+xsize = 5;
+ysize = 5;
+zsize = 10;
+nx = 5;
+ny = 5;
+nz = 10;
+//=============================================================================
+
+// Create a line in the x-direction of length <xsize>, with <nx> divisions
+Point(1) = {-xsize/2, -ysize/2, 0, 0.01};
+Point(2) = {xsize/2, -ysize/2, 0, 0.01};
+Line(1) = {1, 2};
+Transfinite Line(1) = nx+1;
+
+// Extrude split line into meshed square/rectangle
+sq = Extrude {0,ysize,0} {Curve{1}; Layers{ny}; Recombine;};
+
+// Extrude square/rectangle into a cuboid
+cbd = Extrude {0,0,zsize} {Surface{sq[1]}; Layers{nz}; Recombine;};
+
+// Define physical volume, surfaces for BCs
+// Domain
+Physical Volume(0) = {cbd[1]};
+// Low-x side
+Physical Surface(1) = {cbd[5]};
+// High-x side
+Physical Surface(2) = {cbd[3]};
+// Low-y side
+Physical Surface(3) = {cbd[2]};
+// High-y side
+Physical Surface(4) = {cbd[4]};
+// Low-z side
+Physical Surface(5) = {sq[1]};
+// High-z side
+Physical Surface(6) = {cbd[0]};

examples/H3LAPD/3D-hw/hw.xml

@@ -0,0 +1,106 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+
+    <COLLECTIONS DEFAULT="MatrixFree" />
+
+    <!--
+        The composite index for the domain is expected to be 0 if the mesh was generated from the included .geo file
+     -->
+    <EXPANSIONS>
+       <E COMPOSITE="C[0]" NUMMODES="7" TYPE="MODIFIED" FIELDS="ne,w,phi" />
+    </EXPANSIONS>
+
+    <CONDITIONS>
+        <SOLVERINFO>
+            <I PROPERTY="EQTYPE"                   VALUE="3DHW"                 />
+            <I PROPERTY="AdvectionType"            VALUE="WeakDG"               />
+            <I PROPERTY="Projection"               VALUE="DisContinuous"        />
+            <I PROPERTY="TimeIntegrationMethod"    VALUE="ClassicalRungeKutta4" />
+            <I PROPERTY="UpwindType"               VALUE="Upwind"               />
+        </SOLVERINFO>
+
+        <GLOBALSYSSOLNINFO>
+            <V VAR="ne,w,phi">
+                <I PROPERTY="GlobalSysSoln" VALUE="IterativeStaticCond" />
+                <I PROPERTY="IterativeSolverTolerance" VALUE="1e-6"/>
+            </V>
+        </GLOBALSYSSOLNINFO>
+
+        <PARAMETERS>
+            <!-- Timestepping and output options -->
+            <P> TimeStep      = 0.00125           </P>
+            <P> NumSteps      = 32000             </P>
+            <P> TFinal        = NumSteps*TimeStep </P>
+            <P> IO_InfoSteps  = NumSteps/1600     </P>
+            <P> IO_CheckSteps = NumSteps/160      </P>
+            <!-- Magnetic field strength -->
+            <P> Bxy      = 1.0 </P>
+            <!-- d22 Coeff for Helmholtz solve -->
+            <P> d22      = 0.0 </P> 
+            <!-- HW params -->
+            <P> HW_omega_ce = 0.1 </P>
+            <P> HW_nu_ei    = 1.0 </P>
+            <P> HW_kappa    = 3.5 </P>
+            <!-- Scaling factor for ICs -->
+            <P> s        = 0.5 </P>
+            <!-- No particles -->
+            <P> num_particles_total = 0 </P>
+            <!-- Turn on energy, enstrophy output -->
+            <!-- <P> growth_rates_recording_step = 1 </P> -->
+        </PARAMETERS>
+
+        <VARIABLES>
+            <V ID="0"> ne  </V>
+            <V ID="1"> w   </V>
+            <V ID="2"> phi </V>
+        </VARIABLES>
+
+        <BOUNDARYREGIONS>
+            <B ID="0"> C[1] </B> <!-- Low x -->
+            <B ID="1"> C[2] </B> <!-- High x -->
+            <B ID="2"> C[3] </B> <!-- Low y -->
+            <B ID="3"> C[4] </B> <!-- High y -->
+            <B ID="4"> C[5] </B> <!-- Low-z end -->
+            <B ID="5"> C[6] </B> <!-- High-z end -->
+        </BOUNDARYREGIONS>
+
+        <!-- Periodic conditions for all fields on all boundaries -->
+        <BOUNDARYCONDITIONS>
+            <REGION REF="0">
+                <P VAR="ne"  VALUE="[1]" />
+                <P VAR="w"   VALUE="[1]" />
+                <P VAR="phi" VALUE="[1]" />
+            </REGION>
+            <REGION REF="1">
+                <P VAR="ne"  VALUE="[0]" />
+                <P VAR="w"   VALUE="[0]" />
+                <P VAR="phi" VALUE="[0]" />
+            </REGION>
+	        <REGION REF="2">
+                <P VAR="ne"  VALUE="[3]" />
+                <P VAR="w"   VALUE="[3]" />
+                <P VAR="phi" VALUE="[3]" />
+            </REGION>
+            <REGION REF="3">
+                <P VAR="ne"  VALUE="[2]" />
+                <P VAR="w"   VALUE="[2]" />
+                <P VAR="phi" VALUE="[2]" />
+            </REGION>
+            <REGION REF="4">
+                <P VAR="ne"  VALUE="[5]" />
+                <P VAR="w"   VALUE="[5]" />
+                <P VAR="phi" VALUE="[5]" />
+            </REGION>
+            <REGION REF="5">
+                <P VAR="ne"  VALUE="[4]" />
+                <P VAR="w"   VALUE="[4]" />
+                <P VAR="phi" VALUE="[4]" />
+            </REGION>
+        </BOUNDARYCONDITIONS>
+        <FUNCTION NAME="InitialConditions">
+            <E VAR="ne"  DOMAIN="0" VALUE="6*exp((-x*x-y*y)/(s*s))*sin(4*PI*z/10)" />
+            <E VAR="w"   DOMAIN="0" VALUE="(4*exp((-x*x-y*y)/(s*s))*(-s*s+x*x+y*y)/s^4)*sin(4*PI*z/10)" />
+            <E VAR="phi" DOMAIN="0" VALUE="exp(-(x*x+y*y)/(s*s))*sin(4*PI*z/10)" />
+        </FUNCTION>
+    </CONDITIONS>
+</NEKTAR>

examples/H3LAPD/3D-hw/run_cmd_template.txt

@@ -0,0 +1 @@
+mpirun -np <NMPI> <SOLVER_EXEC> hw.xml cuboid.xml

solvers/H3LAPD/CMakeLists.txt

@@ -1,7 +1,8 @@
 # Identify source files
 set(H3LAPD_SRC_FILES
-    EquationSystems/DriftReducedSystem.cpp EquationSystems/HW2Din3DSystem.cpp
-    EquationSystems/LAPDSystem.cpp H3LAPD.cpp)
+    EquationSystems/DriftReducedSystem.cpp EquationSystems/HWSystem.cpp
+    EquationSystems/HW2Din3DSystem.cpp EquationSystems/HW3DSystem.cpp
+    EquationSystems/LAPDSystem.cpp)
 
 # ============================== Object library ===============================
 # Put solver specific source in an object library so that tests can use it

solvers/H3LAPD/Diagnostics/GrowthRatesRecorder.hpp

@@ -21,138 +21,159 @@ namespace NESO::Solvers::H3LAPD {
 template <typename T> class GrowthRatesRecorder {
 protected:
   /// Pointer to number density field
-  std::shared_ptr<T> m_n;
+  std::shared_ptr<T> n;
   /// Pointer to electric potential field
-  std::shared_ptr<T> m_phi;
+  std::shared_ptr<T> phi;
   /// Pointer to vorticity field
-  std::shared_ptr<T> m_w;
+  std::shared_ptr<T> w;
 
   /// HW α constant
-  double m_alpha;
+  double alpha;
   /// File handle for recording output
-  std::ofstream m_fh;
+  std::ofstream fh;
   /// HW κ constant
-  double m_kappa;
-  /// Number of quad points associated with fields m_n, m_phi and m_w
-  int m_npts;
+  double kappa;
+  /// Number of quad points associated with fields n, phi and w
+  int npts;
+  // Space dimension of the problem (not of the domain)
+  const int prob_ndims;
   /// MPI rank
-  int m_rank;
+  int rank;
   /// Sets recording frequency (value of 0 disables recording)
-  int m_recording_step;
+  int recording_step;
   /// Pointer to session object
-  const LU::SessionReaderSharedPtr m_session;
+  const LU::SessionReaderSharedPtr session;
 
 public:
-  GrowthRatesRecorder(const LU::SessionReaderSharedPtr session,
+  GrowthRatesRecorder(const LU::SessionReaderSharedPtr session, int prob_ndims,
                       std::shared_ptr<T> n, std::shared_ptr<T> w,
                       std::shared_ptr<T> phi, int npts, double alpha,
                       double kappa)
-      : m_session(session), m_n(n), m_w(w), m_phi(phi), m_alpha(alpha),
-        m_kappa(kappa), m_npts(npts) {
+      : session(session), prob_ndims(prob_ndims), n(n), w(w), phi(phi),
+        alpha(alpha), kappa(kappa), npts(npts) {
 
     // Store recording frequency for convenience
-    m_session->LoadParameter("growth_rates_recording_step", m_recording_step,
-                             0);
+    this->session->LoadParameter("growth_rates_recording_step",
+                                 this->recording_step, 0);
 
     // Store MPI rank for convenience
-    m_rank = m_session->GetComm()->GetRank();
+    this->rank = this->session->GetComm()->GetRank();
+
+    // Fail for anything but prob_ndims=2 or 3
+    NESOASSERT(this->prob_ndims == 2 || this->prob_ndims == 3,
+               "GrowthRatesRecorder: invalid problem dimensionality; expected "
+               "2 or 3.");
 
     // Write file header
-    if ((m_rank == 0) && (m_recording_step > 0)) {
-      m_fh.open("growth_rates.csv");
-      m_fh << "step,E,W,dEdt_exp,dWdt_exp\n";
+    if ((this->rank == 0) && (this->recording_step > 0)) {
+      this->fh.open("growth_rates.csv");
+      this->fh << "step,E,W,dEdt_exp,dWdt_exp\n";
     }
   };
 
   ~GrowthRatesRecorder() {
     // Close file on destruct
-    if ((m_rank == 0) && (m_recording_step > 0)) {
-      m_fh.close();
+    if ((this->rank == 0) && (this->recording_step > 0)) {
+      this->fh.close();
     }
   }
 
   /**
-   * Calculate Energy = 0.5 ∫ (n^2+|∇ϕ|^2) dV
+   * Calculate Energy = 0.5 ∫ (n^2+|∇⊥ϕ|^2) dV
    */
   inline double compute_energy() {
-    Array<OneD, NekDouble> integrand(m_npts);
+    Array<OneD, NekDouble> integrand(this->npts);
     // First, set integrand = n^2
-    Vmath::Vmul(m_npts, m_n->GetPhys(), 1, m_n->GetPhys(), 1, integrand, 1);
-
-    // Compute phi derivs, square them and add to integrand
-    Array<OneD, NekDouble> xderiv(m_npts), yderiv(m_npts), zderiv(m_npts);
-    m_phi->PhysDeriv(m_phi->GetPhys(), xderiv, yderiv, zderiv);
-    Vmath::Vvtvp(m_npts, xderiv, 1, xderiv, 1, integrand, 1, integrand, 1);
-    Vmath::Vvtvp(m_npts, yderiv, 1, yderiv, 1, integrand, 1, integrand, 1);
-    Vmath::Vvtvp(m_npts, zderiv, 1, zderiv, 1, integrand, 1, integrand, 1);
-
-    // integrand *= 0.5
-    Vmath::Smul(m_npts, 0.5, integrand, 1, integrand, 1);
-    return m_n->Integral(integrand);
+    Vmath::Vmul(this->npts, this->n->GetPhys(), 1, this->n->GetPhys(), 1,
+                integrand, 1);
+
+    // Compute ϕ derivs, square them and add to integrand
+    Array<OneD, NekDouble> xderiv(this->npts), yderiv(this->npts),
+        zderiv(this->npts);
+    this->phi->PhysDeriv(this->phi->GetPhys(), xderiv, yderiv, zderiv);
+    Vmath::Vvtvp(this->npts, xderiv, 1, xderiv, 1, integrand, 1, integrand, 1);
+    Vmath::Vvtvp(this->npts, yderiv, 1, yderiv, 1, integrand, 1, integrand, 1);
+    if (this->prob_ndims == 2) {
+      /* Should be ∇⊥ϕ^2, so ∂ϕ/∂z ought to be excluded in both 2D and
+       * 3D, but there's a small discrepancy in 2D without it. Energy 'leaking'
+       * into orthogonal dimension?!? */
+      Vmath::Vvtvp(this->npts, zderiv, 1, zderiv, 1, integrand, 1, integrand,
+                   1);
+    }
+
+    return 0.5 * this->n->Integral(integrand);
   }
 
   /**
    * Calculate Enstrophy = 0.5 ∫ (n-w)^2 dV
    */
   inline double compute_enstrophy() {
-    Array<OneD, NekDouble> integrand(m_npts);
+    Array<OneD, NekDouble> integrand(this->npts);
     // Set integrand = n-w
-    Vmath::Vsub(m_npts, m_n->GetPhys(), 1, m_w->GetPhys(), 1, integrand, 1);
+    Vmath::Vsub(this->npts, this->n->GetPhys(), 1, this->w->GetPhys(), 1,
+                integrand, 1);
     // Set integrand = (n-w)^2
-    Vmath::Vmul(m_npts, integrand, 1, integrand, 1, integrand, 1);
-    // Set integrand = 0.5*(n-w)^2
-    Vmath::Smul(m_npts, 0.5, integrand, 1, integrand, 1);
-    return m_n->Integral(integrand);
+    Vmath::Vmul(this->npts, integrand, 1, integrand, 1, integrand, 1);
+    return 0.5 * this->n->Integral(integrand);
   }
 
   /**
-   * Calculate Gamma_alpha = alpha ∫ (n-phi)^2 dV
+   * Calculate Γα
+   *  In 2D: Γα = α ∫ (n-ϕ)^2 dV​
+   *  In 3D: Γα = α ∫ [∂/∂z(n-ϕ)]^2 dV
    */
   inline double compute_Gamma_a() {
-    Array<OneD, NekDouble> integrand(m_npts);
-    // Set integrand = n - phi
-    Vmath::Vsub(m_npts, m_n->GetPhys(), 1, m_phi->GetPhys(), 1, integrand, 1);
-    // Set integrand = (n - phi)^2
-    Vmath::Vmul(m_npts, integrand, 1, integrand, 1, integrand, 1);
-    // Set integrand = alpha*(n - phi)^2
-    Vmath::Smul(m_npts, m_alpha, integrand, 1, integrand, 1);
-    return m_n->Integral(integrand);
+    Array<OneD, NekDouble> integrand(this->npts);
+    // Set integrand = n - ϕ
+    Vmath::Vsub(this->npts, this->n->GetPhys(), 1, this->phi->GetPhys(), 1,
+                integrand, 1);
+
+    switch (this->prob_ndims) {
+    case 2:
+      // Set integrand = (n - ϕ)^2
+      Vmath::Vmul(this->npts, integrand, 1, integrand, 1, integrand, 1);
+      break;
+    case 3:
+      // Set integrand = ∂/∂z(n-ϕ)
+      this->phi->PhysDeriv(2, integrand, integrand);
+      // Set integrand = [∂/∂z(n-ϕ)]^2
+      Vmath::Vmul(this->npts, integrand, 1, integrand, 1, integrand, 1);
+      break;
+    }
+    return this->alpha * this->n->Integral(integrand);
   }
 
   /**
-   * Calculate Gamma_n = -kappa ∫ n * dphi/dy dV
+   * Calculate Γn = -κ ∫ n * ∂ϕ/∂y dV
    */
   inline double compute_Gamma_n() {
-    Array<OneD, NekDouble> integrand(m_npts);
-
-    // Set integrand = n * dphi/dy
-    m_phi->PhysDeriv(1, m_phi->GetPhys(), integrand);
-    Vmath::Vmul(m_npts, integrand, 1, m_n->GetPhys(), 1, integrand, 1);
+    Array<OneD, NekDouble> integrand(this->npts);
 
-    // Set integrand = -kappa * n * dphi/dy
-    Vmath::Smul(m_npts, -1 * m_kappa, integrand, 1, integrand, 1);
-    return m_n->Integral(integrand);
+    // Set integrand = n * ∂ϕ/∂y
+    this->phi->PhysDeriv(1, this->phi->GetPhys(), integrand);
+    Vmath::Vmul(this->npts, integrand, 1, this->n->GetPhys(), 1, integrand, 1);
+    return -this->kappa * this->n->Integral(integrand);
   }
 
   /**
    * Compute energy, enstrophy and gamma values and output to file
    */
   inline void compute(int step) {
-    if (m_recording_step > 0) {
-      if (step % m_recording_step == 0) {
+    if (this->recording_step > 0) {
+      if (step % this->recording_step == 0) {
 
         const double energy = compute_energy();
         const double enstrophy = compute_enstrophy();
         const double Gamma_n = compute_Gamma_n();
         const double Gamma_a = compute_Gamma_a();
 
         // Write values to file. In Debug, print to stdout too.
-        if (m_rank == 0) {
+        if (this->rank == 0) {
           nprint(step, ",", energy, ",", enstrophy, ",", Gamma_n - Gamma_a, ",",
                  Gamma_n);
-          m_fh << step << "," << std::setprecision(9) << energy << ","
-               << enstrophy << "," << Gamma_n - Gamma_a << "," << Gamma_n
-               << "\n";
+          this->fh << step << "," << std::setprecision(9) << energy << ","
+                   << enstrophy << "," << Gamma_n - Gamma_a << "," << Gamma_n
+                   << "\n";
         }
       }
     }