Fci: Add zonal neumann BC for potential

.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")

include/vorticity.hxx

@@ -129,11 +129,14 @@ private:
   bool phi_sheath_dissipation; ///< Dissipation at the sheath if phi < 0
   bool damp_core_vorticity; ///< Damp axisymmetric component of vorticity
 
+  std::unique_ptr<Laplacian> phiSolver_zonalneumann;
+  bool zonal_neumann;
+
   bool phi_boundary_relax; ///< Relax boundary to zero-gradient
   BoutReal phi_boundary_timescale; ///< Relaxation timescale [normalised]
   BoutReal phi_boundary_last_update; ///< Time when last updated
   bool phi_core_averagey; ///< Average phi core boundary in Y?
-
+  Field3D zeroes;
   bool split_n0; // Split phi into n=0 and n!=0 components
   LaplaceXY* laplacexy; // Laplacian solver in X-Y (n=0)
   Field3D logB;

src/vorticity.cxx

@@ -1,6 +1,9 @@
 
 #include "../include/vorticity.hxx"
 #include "../include/div_ops.hxx"
+#include "../include/hermes_build_config.hxx"
+#include "../include/hermes_utils.hxx"
+
 
 #include <bout/constants.hxx>
 #include <bout/derivs.hxx>
@@ -14,11 +17,6 @@
 using bout::globals::mesh;
 
 namespace {
-BoutReal floor(BoutReal value, BoutReal min) {
-  if (value < min)
-    return min;
-  return value;
-}
 
 Ind3D indexAt(const Field3D& f, int x, int y, int z) {
   int ny = f.getNy();
@@ -216,6 +214,12 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) {
         Curlb_B = 0.0;
       }
     }
+
+    zeroes = 0.0;
+    zeroes.applyBoundary("neumann");
+    mesh->communicate(zeroes);
+    zeroes.applyParallelBoundary("parallel_neumann_o1");
+
   }
 
   if (Options::root()["mesh"]["paralleltransform"]["type"].as<std::string>()
@@ -260,7 +264,14 @@ Vorticity::Vorticity(std::string name, Options& alloptions, Solver* solver) {
   mesh->communicate(logB);
   logB.applyParallelBoundary("parallel_neumann_o2");
 
-
+  zonal_neumann = options["zonal_neumann"]
+      .doc("Do a second laplace solve for the zonal neumann?")
+      .withDefault<bool>(false);
+
+  if (zonal_neumann) {
+    phiSolver_zonalneumann = Laplacian::create(&options["laplacian_zonalneumann"]);
+  }
+
 
 }
 
@@ -420,40 +431,13 @@ void Vorticity::transform(Options& state) {
       Te = 0.0;
     }
 
-    // Sheath multiplier Te -> phi (2.84522 for Deuterium if Ti = 0)
-    if ( mesh->firstX()) {
-      for (int j = mesh->ystart; j <= mesh->yend; j++) {
-        BoutReal teavg = 0.0; // Average Te in Z
-
-	for (int k = 0; k < mesh->LocalNz; k++) {
-	  teavg += Te(mesh->xstart, j, k);
-	}
-	teavg /= mesh->LocalNz;
-	BoutReal phivalue = sheathmult * teavg;
-
-        // Set midpoint (boundary) value
-        for (int k = 0; k < mesh->LocalNz; k++) {
-          phi(mesh->xstart - 1, j, k) = 2. * phivalue - phi(mesh->xstart, j, k);
-
-          // Note: This seems to make a difference, but don't know why.
-          // Without this, get convergence failures with no apparent instability
-          // (all fields apparently smooth, well behaved)
-          phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k);
-        }
-      }
-    }
 
     if ( mesh->lastX()) {
       for (int j = mesh->ystart; j <= mesh->yend; j++) {
-        BoutReal teavg = 0.0; // Average Te in Z
-
-	for (int k = 0; k < mesh->LocalNz; k++) {
-	  teavg += Te(mesh->xend, j, k);
-	}
-	teavg /= mesh->LocalNz;
-	BoutReal phivalue = sheathmult * teavg;
         // Set midpoint (boundary) value
         for (int k = 0; k < mesh->LocalNz; k++) {
+	  BoutReal phivalue = sheathmult * 0.5 * (Te(mesh->xend, j, k) + Te(mesh->xend + 1, j, k));
+
           phi(mesh->xend + 1, j, k) = 2. * phivalue - phi(mesh->xend, j, k);
 
           // Note: This seems to make a difference, but don't know why.
@@ -538,6 +522,23 @@ void Vorticity::transform(Options& state) {
     }
   }
 
+  if (zonal_neumann) {
+    auto coord = mesh->getCoordinates();
+    Field2D avg_phi = DC(phi);
+    if ( mesh->firstX() ) {
+      for (int j = mesh->ystart; j <= mesh->yend; j++) {
+	for (int k = 0; k < mesh->LocalNz; k++) {
+	  phi_plus_pi(mesh->xstart - 1, j, k) = 0.5 * ( avg_phi(mesh->xstart - 1 ,j) + avg_phi(mesh->xstart ,j) ) +
+	    0.5 * (Pi_hat(mesh->xstart - 1, j, k) + Pi_hat(mesh->xstart, j, k));
+	  phi_plus_pi(mesh->xstart - 2, j, k) = phi_plus_pi(mesh->xstart - 1, j, k);
+	}
+      }
+    }
+    phiSolver_zonalneumann->setCoefC(average_atomic_mass / SQ(coord->Bxy));
+    const auto tosolve = Vort * (Bsq / average_atomic_mass);
+    phi = phiSolver_zonalneumann->solve(tosolve, phi_plus_pi) - Pi_hat;
+  }
+
   // Ensure that potential is set in the communication guard cells
   mesh->communicate(phi);
 
@@ -774,11 +775,13 @@ void Vorticity::finally(const Options& state) {
 
     const Field3D N = get<Field3D>(species["density"]);
     const Field3D NV = get<Field3D>(species["momentum"]);
+    const Field3D V = get<Field3D>(species["velocity"]);
     const BoutReal A = get<BoutReal>(species["AA"]);
 
-    // Note: Using NV rather than N*V so that the cell boundary flux is correct
-    const Field3D jpar = (Z / A) * NV;
-    ddt(Vort) += Div_par(jpar);
+
+    Field3D flow_ylow = 0.0;
+    ddt(Vort) += Z * FV::Div_par_mod<hermes::Limiter>(N, V, zeroes, flow_ylow,  false,
+						   false, true);
 
     if (state["fields"].isSet("Apar_flutter")) {
       // Magnetic flutter term
@@ -787,7 +790,7 @@ void Vorticity::finally(const Options& state) {
       // Div_par(jpar) = B * Grad_par(jpar / B)
       // Using the approximation for small delta-B/B
       // b dot Grad(jpar) = Grad_par(jpar) + [jpar, Apar]
-      ddt(Vort) += coord->Bxy * bracket(jpar / coord->Bxy, Apar_flutter, BRACKET_ARAKAWA);
+      ddt(Vort) += coord->Bxy * bracket(Z*NV / coord->Bxy, Apar_flutter, BRACKET_ARAKAWA);
     }
   }
 
@@ -812,7 +815,13 @@ void Vorticity::finally(const Options& state) {
     // Adds dissipation term like in other equations, but depending on gradient of
     // potential
     Field3D sound_speed = get<Field3D>(state["sound_speed"]);
-    ddt(Vort) -= FV::Div_par(-phi, 0.0, sound_speed);
+    Field3D dummy1;
+    zeroes = 0.0;
+    zeroes.applyBoundary("neumann");
+    mesh->communicate(zeroes);
+    zeroes.applyParallelBoundary("parallel_neumann_o1");
+    ddt(Vort) -= FV::Div_par_mod<hermes::Limiter>(-phi, zeroes, sound_speed, dummy1,  false,
+						  false, true);
   }
 
   if (hyper > 0) {

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[ ]