Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implicit implementation of 2Din3D Hasegawa-Wakatani #232

Closed
wants to merge 18 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
c218cfe
Add skeleton for 3D equation system and example; just copied from 2D…
oparry-ukaea Dec 4, 2023
5eccd35
Working 3D implementation and example.
oparry-ukaea Dec 6, 2023
061d6d3
Remove spurious field in 3D e.g
oparry-ukaea Dec 7, 2023
924b5c1
Adapt GrowthRateRecorder for true 3D HW.
oparry-ukaea Dec 7, 2023
4ff1528
Refactor HW equation classes.
oparry-ukaea Dec 7, 2023
1cd6318
Add test for energy, enstrophy growth rates in 3D. Energy slightly of…
oparry-ukaea Dec 7, 2023
f2a00d7
Simplify GrowthRatesRecorder by applying scalar factors after integrals.
oparry-ukaea Dec 8, 2023
ee31225
Clearer comments and docstrings in GrowthRatesRecorder header.
oparry-ukaea Dec 8, 2023
d7c46dc
Allow HW alpha and kappa to be specified via physical params.
oparry-ukaea Feb 21, 2024
acba2d9
Report parameter values in DriftReducedSystem and HWSystem.
oparry-ukaea Jan 30, 2024
c753113
Remove member var m_alpha from 2Din3D system (was masking m_alpha in …
oparry-ukaea Feb 27, 2024
4ea6bc2
Tweak existing (explicit) HW configs to match nektar-driftwave.
oparry-ukaea Apr 4, 2024
05766b8
Add implicit 2Din3DHW eg.
oparry-ukaea Apr 4, 2024
d411449
Use recent nektar version - required for some functionality used in i…
oparry-ukaea Apr 5, 2024
7ef3381
Implicit version of 2Din3D, copied more-or-less wholesale from nektar…
oparry-ukaea Apr 5, 2024
8e9434a
Minor bugfix.
oparry-ukaea Apr 5, 2024
109f920
Make implicit 2Din3D config closer to the working nektar-driftwave ex…
oparry-ukaea Apr 5, 2024
63e9523
Correct access for do_ode_projection
oparry-ukaea Apr 8, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 7 additions & 6 deletions examples/H3LAPD/2Din3D-hw/hw.xml
Original file line number Diff line number Diff line change
Expand Up @@ -12,17 +12,18 @@

<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" />
<I PROPERTY="EQTYPE" VALUE="2Din3DHW" />
<I PROPERTY="AdvectionType" VALUE="WeakDG" />
<I PROPERTY="Projection" VALUE="DisContinuous" />
<I PROPERTY="TimeIntegrationMethod" VALUE="RungeKutta4" />
<I PROPERTY="AdvectionAdvancement" VALUE="Explicit" />
<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"/>
<I PROPERTY="IterativeSolverTolerance" VALUE="1e-6" />
</V>
</GLOBALSYSSOLNINFO>

Expand Down
Original file line number Diff line number Diff line change
@@ -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]};
111 changes: 111 additions & 0 deletions examples/H3LAPD/2Din3D-hw_fluid-only-implicit/hw.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
<?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="4" 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="DIRKOrder2" />
<I PROPERTY="AdvectionAdvancement" VALUE="Implicit" />
<I PROPERTY="LinSysIterSolverTypeInNonlin" VALUE="GMRES" />
</SOLVERINFO>

<GLOBALSYSSOLNINFO>
<V VAR="ne,w,phi,ne_src">
<I PROPERTY="GlobalSysSoln" VALUE="IterativeStaticCond" />
</V>
</GLOBALSYSSOLNINFO>

<PARAMETERS>
<P> NumSteps = 30000 </P>
<P> TimeStep = 0.005 </P>
<P> IO_InfoSteps = NumSteps/1000 </P>
<P> IO_CheckSteps = NumSteps/100 </P>
<P> s = 2.0 </P>
<P> kappa = 1.0 </P>
<P> alpha = 2.0 </P>
<P> IterativeSolverTolerance = 1.0E-6 </P>
<P> NekNonlinSysMaxIterations = 50 </P>
<P> NewtonRelativeIteTol = 1.0E-6 </P>
<P> NonlinIterTolRelativeL2 = 1.0E-6 </P>
<P> LinSysRelativeTolInNonlin = 1.0E-02 </P>
<!-- Magnetic field strength -->
<P> Bxy = 1.0 </P>
<!-- d22 Coeff for Helmholtz solve -->
<P> d22 = 0.0 </P>
<!-- HW params -->
<P> HW_alpha = 0.1 </P>
<P> HW_kappa = 0.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>
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
mpirun -np <NMPI> <SOLVER_EXEC> hw.xml cuboid.xml
13 changes: 7 additions & 6 deletions examples/H3LAPD/2Din3D-hw_fluid-only/hw.xml
Original file line number Diff line number Diff line change
Expand Up @@ -12,17 +12,18 @@

<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" />
<I PROPERTY="EQTYPE" VALUE="2Din3DHW" />
<I PROPERTY="AdvectionType" VALUE="WeakDG" />
<I PROPERTY="Projection" VALUE="DisContinuous" />
<I PROPERTY="TimeIntegrationMethod" VALUE="RungeKutta4" />
<I PROPERTY="AdvectionAdvancement" VALUE="Explicit" />
<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"/>
<I PROPERTY="IterativeSolverTolerance" VALUE="1e-6" />
</V>
</GLOBALSYSSOLNINFO>

Expand Down
37 changes: 37 additions & 0 deletions examples/H3LAPD/3D-hw/cuboid_periodic_5x5x10.geo
Original file line number Diff line number Diff line change
@@ -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]};
107 changes: 107 additions & 0 deletions examples/H3LAPD/3D-hw/hw.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
<?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="RungeKutta4" />
<I PROPERTY="AdvectionAdvancement" VALUE="Explicit" />
<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>
1 change: 1 addition & 0 deletions examples/H3LAPD/3D-hw/run_cmd_template.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
mpirun -np <NMPI> <SOLVER_EXEC> hw.xml cuboid.xml
2 changes: 1 addition & 1 deletion nektar
Submodule nektar updated from 2e0fb8 to 8bc4b3
3 changes: 2 additions & 1 deletion solvers/H3LAPD/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
# Identify source files
set(H3LAPD_SRC_FILES
EquationSystems/DriftReducedSystem.cpp EquationSystems/HW2Din3DSystem.cpp
EquationSystems/DriftReducedSystem.cpp EquationSystems/HWSystem.cpp
EquationSystems/HW2Din3DSystem.cpp EquationSystems/HW3DSystem.cpp
EquationSystems/LAPDSystem.cpp H3LAPD.cpp)

# ============================== Object library ===============================
Expand Down
Loading