Fci: Add 3D MMS tests and fix neutral averaging

.github/workflows/tests.yml

@@ -56,7 +56,7 @@ jobs:
 
       - name: Install pip packages
         run: |
-          ./.pip_install_for_tests.sh 'netcdf4~=1.5' 'boutdata==0.3.0' boututils
+          ./.pip_install_for_tests.sh 'netcdf4~=1.5' 'boutdata==0.3.0' boututils 'git+https://github.com/boutproject/zoidberg'
           # Add the pip install location to the runner's PATH
           echo ~/.local/bin >> $GITHUB_PATH
       - name: Build and run tests

CMakeLists.txt

@@ -245,6 +245,14 @@ add_custom_command(TARGET setup-test PRE_BUILD
 if(HERMES_TESTS)
   enable_testing()
 
+  if (BOUT_ENABLE_METRIC_3D)
+    add_custom_command(TARGET setup-test PRE_BUILD
+      COMMAND python3 ./generate_axial_circular.py
+      WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/grids
+    )
+    file(CREATE_LINK ${CMAKE_CURRENT_SOURCE_DIR}/grids ${CMAKE_CURRENT_BINARY_DIR}/grids SYMBOLIC)
+  endif()
+
   # Integrated tests
   function(hermes_add_integrated_test TESTNAME)
     set(oneValueArgs DOWNLOAD DOWNLOAD_NAME)
@@ -298,6 +306,10 @@ if(HERMES_TESTS)
 	DEPENDS setup-test)
     endif()
   endfunction()
+  hermes_add_integrated_test(3D-Axial-circular-conduction-fci
+    REQUIRES BOUT_ENABLE_METRIC_3D)
+  hermes_add_integrated_test(2D-axial-circular-fci
+    REQUIRES BOUT_ENABLE_METRIC_3D)
   hermes_add_integrated_test(2D-vorticity-inversion-fci
     REQUIRES BOUT_ENABLE_METRIC_3D)
   hermes_add_integrated_test(1D-sheathtest-recycling-fci

grids/generate_axial_circular.py

@@ -0,0 +1,92 @@
+import zoidberg 
+from zoidberg.field import Slab
+import numpy as np
+
+import os
+import sys
+
+class Axial_circular(Slab):
+    def __init__(self,rho_1,rho_2,Lz,By,R0, q0 = 3.0, q1 = 0.2):
+        self.rho_1 = float(rho_1)
+        self.rho_2 = float(rho_2)
+        self.Lz = float(Lz)
+        self.By = float(By)
+        self.R0 = R0
+
+    def qprofile(self,x,z,phi):
+        theta=np.arctan2(z,x-self.R0)
+        r = np.sqrt(np.power(x-self.R0,2)+np.power(z,2))
+        return q0 + q1 * r
+
+    def Bxfunc(self,x,z,phi):
+        theta=np.arctan2(z,x-self.R0)
+        r = np.sqrt(np.power(x-self.R0,2)+np.power(z,2))
+        q = self.qprofile(x,z,phi)
+        return -r * np.sin(theta) / q
+
+    def Bzfunc(self,x,z,phi):
+        theta=np.arctan2(z,x-self.R0)
+        r = np.sqrt(np.power(x-self.R0,2)+np.power(z,2))
+        q = self.qprofile(x,z,phi)
+        return r * np.cos(theta) / q   
+
+    def Byfunc(self,x,z,phi):
+        return np.full(x.shape,self.By)
+
+
+
+
+R0 = 1.5
+rho_1 = 0.4
+rho_2 = 0.6
+B0 = 1.0
+Ly = 2.0 * np.pi
+q0 = 3.0
+q1 = 0.0
+field = Axial_circular(rho_1, rho_2, Ly, B0,R0, q0=q0, q1=q1)
+
+folder = "mms_axial_circular_fci/"
+
+if not os.path.exists(folder):
+    os.mkdir(folder)
+
+force = "-f" in sys.argv or "--force" in sys.argv
+
+
+
+for scale in [1,2]:
+    nx = scale * 16
+    nz = scale * 64
+    ny = scale * 8
+
+    filename = f"Axial_circular_q_{nx}_{nz}_{rho_1}_{rho_2}_{q0}_{q1}.fci.grid.nc"
+    fn = folder + filename
+
+    if os.path.exists(fn) and not force:
+        print(fn, " exists")
+        continue
+
+
+    inner = zoidberg.rzline.shaped_line(R0=R0,a=rho_1,n=nz)
+    outer = zoidberg.rzline.shaped_line(R0=R0,a=rho_2,n=nz)
+
+    pol_grid = zoidberg.poloidal_grid.grid_annulus(inner,outer,nx+4,nz)
+    pol_grids = []
+    for i in range(ny):
+        pol_grids.append(pol_grid)
+
+    grid = zoidberg.grid.Grid(pol_grids, np.linspace(0.0, 2.0 * np.pi, ny, endpoint=False),Ly,yperiodic=True)
+    maps = zoidberg.zoidberg.make_maps(grid, field)
+
+    maps["forward_xt_prime"][2,:,:] = 2.0
+    maps["backward_xt_prime"][2,:,:] = 2.0
+
+    maps["forward_xt_prime"][-3,:,:] = float(nx+4-3)
+    maps["backward_xt_prime"][-3,:,:] = float(nx+4-3)
+
+    with zoidberg.zoidberg.MapWriter(fn) as mw:
+        mw.add_grid_field(grid, field)
+        mw.add_maps(maps)
+        mw.add_dagp()
+
+    print(fn, " done")

hermes-3.cxx

@@ -273,7 +273,7 @@ int Hermes::init(bool restarting) {
 
 
 	// try loading J from the grid, otherwise use the one calculated from the metric coefficients
-	Field3D Jtmp = 0.0;
+	Coordinates::FieldMetric Jtmp = 0.0;
 	if (mesh->get(Jtmp, "J_new")==0){
 	  mesh->communicate(Jtmp);
 	  coord->J = Jtmp;
@@ -291,7 +291,7 @@ int Hermes::init(bool restarting) {
 	coord->J_perp = sqrt(coord->g_11 * coord->g_33 - coord->g_13 * coord->g_13);
 
 	// try loading g_22 at the lower and upper cell interface from the grid, otherwise caluculate from the mean of the two cellcentered ones                                                           
-        Field3D loadtmp = 0.0;
+        Coordinates::FieldMetric loadtmp = 0.0;
         if (mesh->get(loadtmp, "g_22_cell_ylow")==0) {
           coord->g_22_cell_ylow = loadtmp / SQ(rho_s0);
         } else {

include/amjuel_reaction.hxx

@@ -102,16 +102,9 @@ protected:
       BoutReal avgNe = 0.0;
       BoutReal avgN1 = 0.0;
       BoutReal avgTe = 0.0;
-      if (Ne.isFci()) {
-	avgNe = (4.0 * Ne[i] + Ne.ydown()[iym] + Ne.yup()[iyp]) / 6.0;
-	avgN1 = (4.0 * N1[i] + N1.ydown()[iym] + N1.yup()[iyp]) / 6.0;
-	avgTe = (4.0 * Te[i] + Te.ydown()[iym] + Te.yup()[iyp]) / 6.0;
-      } else {
-	avgNe = Ne[i];
-	avgN1 = N1[i];
-	avgTe = Te[i];
-      }
-
+      avgNe = Ne[i];
+      avgN1 = N1[i];
+      avgTe = Te[i];      
       reaction_rate[i] = avgNe * avgN1 * evaluate(rate_coefs, avgTe * Tnorm, avgNe * Nnorm) * Nnorm / FreqNorm * rate_multiplier;
     }
 
@@ -179,15 +172,11 @@ protected:
       BoutReal avgNe = 0.0;
       BoutReal avgN1 = 0.0;
       BoutReal avgTe = 0.0;
-      if (Ne.isFci()) {
-        avgNe = (4.0 * Ne[i] + Ne.ydown()[iym] + Ne.yup()[iyp]) / 6.0;
-        avgN1 = (4.0 * N1[i] + N1.ydown()[iym] + N1.yup()[iyp]) / 6.0;
-        avgTe = (4.0 * Te[i] + Te.ydown()[iym] + Te.yup()[iyp]) / 6.0;
-      }	else {
-	avgNe =	Ne[i];
-        avgN1 =	N1[i];
-        avgTe =	Te[i];
-      }
+
+      avgNe =	Ne[i];
+      avgN1 =	N1[i];
+      avgTe =	Te[i];
+
       energy_loss[i] = avgNe * avgN1 * evaluate(radiation_coefs, avgTe * Tnorm, avgNe * Nnorm) * Nnorm / (Tnorm * FreqNorm) * radiation_multiplier;
     }
 

tests/integrated/2D-axial-circular-fci/axial.m

@@ -0,0 +1,72 @@
+(* ::Package:: *)
+
+BeginPackage["axial`"]
+
+(* To run the package, execute the following line in a notebook
+<<axial.m
+ *)
+
+(* 
+"Provides MMS solution and source for the anisotropic diffusion model in circular geometry
+- Only Dirichlet boundary conditions in the perpendicular direction is considered
+- Parallel boundaries (Penalisation) are not taken into account 
+- MMS solution is firstly is prescribed in cylindrical coordinates 
+(r,p,z) and afterwards a transformed to Cartesian coordinates (x,y,phi) used in GRILLIX 
+r = sqrt(x^2+y^2);  p = arctan(y/x);  z = phi, t = time
+- Paramters of model are: chi_par, chi_perp, safety factor q, and limiting flux surfaces rmin and rmax"
+*)
+
+
+Print["Computing MMS Terms"];
+
+(*
+define x,y and parallel derivatives in terms of r,p,z derivative \
+and normalised radial coordinate rn
+*)
+absb[x_] = Sqrt[1 + x^2/q[x]^2];
+pgrad[f_, x_, z_, y_, t_] = (D[f[x,z,y,t],y] + 1/q[x]*D[f[x,z,y,t],z])/absb[x];
+ddx[f_, r_, p_, z_, t_] = D[f[r, p, z, t], r];
+ddy[f_, r_, p_, z_, t_] = D[f[r, p, z, t], r]*Sin[p] + D[f[r, p, z, t], p]*Cos[p]/r;
+d2dx2[f_, r_, p_, z_, t_] = D[ddx[f, r, p, z, t], r];
+d2dy2[f_, r_, p_, z_, t_] = D[ddy[f, r, p, z, t], r]*Sin[p] +  D[ddy[f, r, p, z, t], p]*Cos[p]/r;
+
+
+LaplacePerpMmsSol[f_,r_, p_, z_, t_] =  d2dx2[f, r, p, z, t] + d2dy2[f, r, p, z, t];
+
+LaplacePerp[f_,r_,p_,z_,t_] = D[f[r,p,z,t],{r,2}] + D[f[r,p,z,t],r]/r + 1.0/(r*r)*D[f[r,p,z,t],{p,2}];
+
+
+
+xn[x_] = (x-xmin)/(xmax-xmin);  
+
+d2dpar2[f_, x_, z_, y_, t_] = (D[D[f[x, z, y, t], y], y] + 2/q[x]*D[D[f[x, z, y, t], y], z] +  1/q[x]^2*D[D[f[x, z, y, t], z], z])/absb[x]^2;
+
+
+
+(*
+Define normalised rho and MMS solution in terms of mode numbers \
+given above
+*)
+MmsDens[x_, z_, y_, t_] = amp*Sin[2.0*Pi*kx*xn[x]]*Sin[kz*z - phz]*Cos[ky*y- phy]+offset;
+MmsUpar[x_, z_, y_, t_]=1;
+rhos = 0.0002284697436697996
+Bnorm = 1.0
+qe = 1.60217663*^-19
+Me = 9.1093837*^-31
+e0 = 8.85418781*^-12
+Omegaci = qe * Bnorm / (1836.0*Me)
+
+
+pflux[x_, z_, y_, t_]=MmsDens[x, z, y, t];
+(*Smms[x_, z_, y_, t_]=D[MmsDens[x,z,y,t],t]-d2dpar2[MmsDens,x,z,y,t]-Laplaceperpe[MmsDens,x,z,y,t];*)
+Smms[x_, z_, y_, t_]=D[MmsDens[x,z,y,t],t]-Dperp/(rhos*rhos*Omegaci) * LaplacePerp[MmsDens,x,z,y,t] * (rhos*rhos)-
+					Dpar/(rhos*rhos*Omegaci)*d2dpar2[MmsDens,x,z,y,t]* (rhos*rhos);
+
+
+Print["Finished MMS Terms"];
+Print[Omegaci]
+(* Set a dummy return value *)
+1
+
+
+EndPackage[ ]