Add 3D example, read 3D HW params.

ExCALIBUR-NEPTUNE · Dec 4, 2023 · ea72b3a · ea72b3a
1 parent 1ce10cd
commit ea72b3a
Show file tree

Hide file tree

Showing 6 changed files with 162 additions and 18 deletions.
diff --git a/examples/H3LAPD/3D-hw/cuboid_periodic_5x5x10.geo b/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]};
diff --git a/examples/H3LAPD/3D-hw/hw.xml b/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="2Din3DHW"                 />
+            <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,ne_src">
+                <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> omega_ce = 1.0 </P>
+            <P> 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>
diff --git a/examples/H3LAPD/3D-hw/run_cmd_template.txt b/examples/H3LAPD/3D-hw/run_cmd_template.txt
@@ -0,0 +1 @@
+mpirun -np <NMPI> <SOLVER_EXEC> hw.xml cuboid.xml
diff --git a/solvers/H3LAPD/CMakeLists.txt b/solvers/H3LAPD/CMakeLists.txt
@@ -1,7 +1,7 @@
 # Identify source files
 set(H3LAPD_SRC_FILES
     EquationSystems/DriftReducedSystem.cpp EquationSystems/HW2Din3DSystem.cpp
-    EquationSystems/LAPDSystem.cpp H3LAPD.cpp)
+    EquationSystems/HW3DSystem.cpp EquationSystems/LAPDSystem.cpp H3LAPD.cpp)
 
 # ============================== Object library ===============================
 # Put solver specific source in an object library so that tests can use it

diff --git a/solvers/H3LAPD/EquationSystems/HW3DSystem.cpp b/solvers/H3LAPD/EquationSystems/HW3DSystem.cpp
@@ -85,14 +85,6 @@ void HW3DSystem::explicit_time_int(
   add_adv_terms({"ne"}, m_adv_elec, m_adv_vel_elec, in_arr, out_arr, time);
   add_adv_terms({"w"}, m_adv_vort, m_ExB_vel, in_arr, out_arr, time);
 
-  // Add \alpha*(\phi-n_e) to RHS
-  Array<OneD, NekDouble> HWterm_2D_alpha(npts);
-  Vmath::Vsub(npts, m_fields[phi_idx]->GetPhys(), 1,
-              m_fields[ne_idx]->GetPhys(), 1, HWterm_2D_alpha, 1);
-  Vmath::Smul(npts, m_omega_ce, HWterm_2D_alpha, 1, HWterm_2D_alpha, 1);
-  Vmath::Vadd(npts, out_arr[w_idx], 1, HWterm_2D_alpha, 1, out_arr[w_idx], 1);
-  Vmath::Vadd(npts, out_arr[ne_idx], 1, HWterm_2D_alpha, 1, out_arr[ne_idx], 1);
-
   // Add \kappa*\dpartial\phi/\dpartial y to RHS
   Array<OneD, NekDouble> HWterm_2D_kappa(npts);
   m_fields[phi_idx]->PhysDeriv(1, m_fields[phi_idx]->GetPhys(),
@@ -124,10 +116,14 @@ void HW3DSystem::get_phi_solve_rhs(
 void HW3DSystem::load_params() {
   DriftReducedSystem::load_params();
 
-  // alpha
-  m_session->LoadParameter("HW_alpha", m_alpha, 2);
-  // kappa
-  m_session->LoadParameter("HW_kappa", m_kappa, 1);
+  // kappa (required)
+  m_session->LoadParameter("HW_kappa", m_kappa);
+
+  // ω_ce (required)
+  m_session->LoadParameter("HW_omega_ce", m_omega_ce);
+
+  // ν_ei (required)
+  m_session->LoadParameter("HW_nu_ei", m_nu_ei);
 }
 
 /**

diff --git a/solvers/H3LAPD/EquationSystems/HW3DSystem.hpp b/solvers/H3LAPD/EquationSystems/HW3DSystem.hpp
@@ -76,16 +76,20 @@ class HW3DSystem : virtual public DriftReducedSystem {
   virtual bool v_PreIntegrate(int step) override;
 
 private:
-  /// Hasegawa-Wakatani α
-  NekDouble m_alpha;
   /// Bool to enable/disable growth rate recordings
   bool m_diag_growth_rates_recording_enabled;
   /// Bool to enable/disable mass recordings
   bool m_diag_mass_recording_enabled;
-  /// Hasegawa-Wakatani κ - not clear whether this should be included in 3D or
-  /// not
+  /**
+   * Hasegawa-Wakatani κ - not clear whether this should be included in 3D or
+   * not
+   */
   NekDouble m_kappa;
+  /// Electron-ion collision frequency
+  NekDouble m_nu_ei;
+  /// Cyclotron frequency for electrons
+  NekDouble m_omega_ce;
 };
 
 } // namespace NESO::Solvers::H3LAPD
-#endif // H3LAPD_HW2DIN3D_SYSTEM_H
+#endif // H3LAPD_HW3D_SYSTEM_H
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1 @@
		mpirun -np <NMPI> <SOLVER_EXEC> hw.xml cuboid.xml