Impurity test and electromagnetic component

CMakeLists.txt

@@ -209,6 +209,32 @@ configure_file(
   "${PROJECT_BINARY_DIR}/include/revision.hxx")
 
 
+
+
+# Once built, copy the data and test directories
+option(HERMES_COPY_EXAMPLES_TO_BUILD "Copy the examples to the build directory"
+       ON)
+
+option(HERMES_COPY_JSON_DATABASE_TO_BUILD
+       "Copy the json database to the build directory" ON)
+
+if(HERMES_COPY_EXAMPLES_TO_BUILD OR HERMES_COPY_JSON_DATABASE_TO_BUILD)
+  add_custom_target(setup-examples ALL)
+else()
+  add_custom_target(setup-examples)
+endif()
+
+
+add_custom_command(TARGET setup-examples PRE_BUILD
+                   COMMAND ${CMAKE_COMMAND} -E copy_directory
+                   ${CMAKE_SOURCE_DIR}/examples $<TARGET_FILE_DIR:${PROJECT_NAME}>/examples)
+
+add_custom_command(TARGET setup-examples PRE_BUILD
+                   COMMAND ${CMAKE_COMMAND} -E copy_directory
+                   ${CMAKE_SOURCE_DIR}/json_database $<TARGET_FILE_DIR:${PROJECT_NAME}>/json_database)
+
+
+
 # Tests
 option(HERMES_COPY_TESTS_TO_BUILD "Copy the tests to the build directory" ON)
 if(HERMES_COPY_TESTS_TO_BUILD)
@@ -225,23 +251,17 @@ else()
   add_custom_target(setup-test)
 endif()
 
-add_custom_target(setup-examples)
-
+# add_custom_target(setup-examples)
 
-# copy the data and test directories
-add_custom_command(TARGET setup-examples PRE_BUILD
-                   COMMAND ${CMAKE_COMMAND} -E copy_directory
-                   ${CMAKE_SOURCE_DIR}/examples $<TARGET_FILE_DIR:${PROJECT_NAME}>/examples)
-
-add_custom_command(TARGET setup-examples PRE_BUILD
-                   COMMAND ${CMAKE_COMMAND} -E copy_directory
-                   ${CMAKE_SOURCE_DIR}/json_database $<TARGET_FILE_DIR:${PROJECT_NAME}>/json_database)
 
 add_custom_command(TARGET setup-test PRE_BUILD
                    COMMAND ${CMAKE_COMMAND} -E copy_directory
                    ${CMAKE_SOURCE_DIR}/tests $<TARGET_FILE_DIR:${PROJECT_NAME}>/tests)
 
 
+
+
+
 if(HERMES_TESTS)
   enable_testing()
 
@@ -298,6 +318,8 @@ if(HERMES_TESTS)
 	DEPENDS setup-test)
     endif()
   endfunction()
+  hermes_add_integrated_test(1D-Cimpurities-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

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;
     }
 

include/electromagnetic.hxx

@@ -66,14 +66,14 @@ private:
   BoutReal beta_em; // Normalisation coefficient mu_0 e T n / B^2
 
   std::unique_ptr<Laplacian> aparSolver; // Laplacian solver in X-Z
-
+  bool use_normdensity;
   bool const_gradient; // Set neumann boundaries by extrapolation
   BoutReal apar_boundary_timescale; // Relaxation timescale
   BoutReal last_time;  // The last time the boundaries were updated
 
   bool magnetic_flutter; ///< Set the magnetic flutter term?
   Field3D Apar_flutter;
-
+  Field3D zeroes;
   bool diagnose; ///< Output additional diagnostics?
 };
 

include/evolve_density.hxx

@@ -64,7 +64,7 @@ private:
   BoutReal AA;          ///< Atomic mass e.g. proton = 1
 
   Field3D N;            ///< Species density (normalised, evolving)
-
+  bool magnetic_flutter;
   bool bndry_flux;      ///< Allow flows through boundaries?
   bool exb_advection;   ///< Include ExB advection?
   bool poloidal_flows;  ///< Include ExB flow in Y direction?

include/evolve_momentum.hxx

@@ -50,7 +50,7 @@ private:
   BoutReal pressure_floor;
   bool low_p_diffuse_perp; ///< Cross-field diffusion at low pressure?
   BoutReal scale_ExB;
-
+  bool magnetic_flutter;
   BoutReal hyper_z;  ///< Hyper-diffusion
   BoutReal hyper_nv;
   std::string Vname;

include/evolve_pressure.hxx

@@ -84,6 +84,8 @@ private:
   bool thermal_conduction;    ///< Include thermal conduction?
   BoutReal kappa_coefficient; ///< Leading numerical coefficient in parallel heat flux calculation
   BoutReal kappa_limit_alpha; ///< Flux limit if >0
+  bool kappa_limit_grillix;
+  BoutReal kappa_limit_Lpar;
   bool disable_ddt;
   bool p_div_v; ///< Use p*Div(v) form? False -> v * Grad(p)
   BoutReal T_lowsource;
@@ -103,7 +105,7 @@ private:
   Field3D source, final_source; ///< External pressure source
   Field3D Sp;     ///< Total pressure source
   FieldGeneratorPtr source_prefactor_function;
-
+  bool magnetic_flutter;
   BoutReal adapt_source;
 
   BoutReal hyper_z; ///< Hyper-diffusion

include/relax_potential.hxx

@@ -90,6 +90,9 @@ private:
   Field3D viscosity_core;
   Field3D viscosity_par;
   Field3D ones;
+  Field3D zeroes;
+  Field3D is_SOL;
+  bool core_dissipation;
   bool phi_dissipation; /// Parallel dissipation of potential
   BoutReal vort_timedissipation;
 

include/vorticity.hxx

@@ -120,20 +120,25 @@ private:
   BoutReal average_atomic_mass; // Weighted average atomic mass, for polarisaion current (Boussinesq approximation)
   bool poloidal_flows;   ///< Include poloidal ExB flow?
   bool bndry_flux;  ///< Allow flows through radial boundaries?
-
+  bool boussinesq;
   bool collisional_friction; ///< Damping of vorticity due to collisional friction
 
+  std::unique_ptr<Laplacian> phiSolver_zonalneumann;
+
+  bool zonal_neumann;
+
   bool sheath_boundary; ///< Set outer boundary to j=0?
   bool vort_dissipation; ///< Parallel dissipation of vorticity
   bool phi_dissipation;  ///< Parallel dissipation of potential
+  BoutReal phi_diss_factor;
   bool phi_sheath_dissipation; ///< Dissipation at the sheath if phi < 0
   bool damp_core_vorticity; ///< Damp axisymmetric component of vorticity
 
   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?
-
+  bool magnetic_flutter;
   bool split_n0; // Split phi into n=0 and n!=0 components
   LaplaceXY* laplacexy; // Laplacian solver in X-Y (n=0)
   Field3D logB;

src/anomalous_diffusion_3d.cxx

@@ -97,7 +97,7 @@ void AnomalousDiffusion3D::transform(Options& state) {
                         ? GET_NOBOUNDARY(Field3D, species["temperature"])
                         : 0.0;
   Field3D V =
-      species.isSet("velocity") ? GET_NOBOUNDARY(Field3D, species["velocity"]) : 0.0;
+    (species.isSet("velocity") && include_nu)? GET_NOBOUNDARY(Field3D, species["velocity"]) : 0.0;
 
   Field3D flow_xlow, flow_zlow; // Flows through cell faces
 

src/electromagnetic.cxx

@@ -38,6 +38,10 @@ Electromagnetic::Electromagnetic(std::string name, Options &alloptions, Solver*
     .doc("Extrapolate gradient of Apar into all radial boundaries?")
     .withDefault<bool>(false);
 
+  use_normdensity = options["use_normdensity"]
+    .doc("Use the normalized density instead of the full density in the Apar equation?")
+    .withDefault<bool>(false);
+
   // Give Apar an initial value because we solve Apar by iteration
   // starting from the previous solution
   // Note: On restart the value is restored (if available) in restartVars
@@ -73,6 +77,11 @@ Electromagnetic::Electromagnetic(std::string name, Options &alloptions, Solver*
     .doc("Output additional diagnostics?")
     .withDefault<bool>(false);
 
+  zeroes = 0.0;
+  zeroes.applyBoundary("neumann");
+  bout::globals::mesh->communicate(zeroes);
+  zeroes.applyParallelBoundary("parallel_neumann_o1");
+
   magnetic_flutter = options["magnetic_flutter"]
     .doc("Set magnetic flutter terms (Apar_flutter)?")
     .withDefault<bool>(false);
@@ -127,15 +136,21 @@ void Electromagnetic::transform(Options &state) {
     const BoutReal A = get<BoutReal>(species["AA"]);
 
     // Coefficient in front of A_||
-    alpha_em += floor(N, 1e-5) * (SQ(Z) / A);
 
+    if (use_normdensity) {
+      alpha_em += (SQ(Z) / A); 
+    } else {
+      alpha_em += floor(N, 1e-5) * (SQ(Z) / A);
+    }
+
     // Right hand side
     Ajpar += mom * (Z / A);
   }
 
   // Invert Helmholtz equation for Apar
   aparSolver->setCoefA((-beta_em) * alpha_em);
-
+  aparSolver->setCoefC(1.0);
+  //aparSolver->setCoefA(0.0);
   if (const_gradient) {
     // Set gradient boundary condition from gradient inside boundary
     Field3D rhs = (-beta_em) * Ajpar;
@@ -173,7 +188,7 @@ void Electromagnetic::transform(Options &state) {
     // Use previous value of Apar as initial guess
     Apar = aparSolver->solve(rhs, Apar);
   } else {
-    Apar = aparSolver->solve((-beta_em) * Ajpar, Apar);
+    Apar = aparSolver->solve((-beta_em) * Ajpar, zeroes);
   }
 
   // Save in the state
@@ -194,50 +209,60 @@ void Electromagnetic::transform(Options &state) {
     const Field3D N = GET_NOBOUNDARY(Field3D, species["density"]);
 
     Field3D nv = getNonFinal<Field3D>(species["momentum"]);
-    nv -= Z * N * Apar;
+    if (use_normdensity) {
+      nv -= Z * Apar;
+    } else {
+      nv -= Z * N * Apar;
+    }
     // Note: velocity is momentum / (A * N)
     Field3D v = getNonFinal<Field3D>(species["velocity"]);
-    v -= (Z / A) * N * Apar / floor(N, 1e-5);
-    // Need to update the guard cells
+    if (use_normdensity) {
+      v -= (Z / A) * Apar / floor(N, 1e-5);
+    } else {
+      v -= (Z / A) * N * Apar / floor(N, 1e-5);
+    }  
+
     nv.applyBoundary("neumann");
     v.applyBoundary("neumann");
     bout::globals::mesh->communicate(nv, v);
     v.applyParallelBoundary("parallel_neumann_o1");
     nv.applyParallelBoundary("parallel_neumann_o1");
-
 
     set(species["momentum"], nv);
     set(species["velocity"], v);
   }
 
   if (magnetic_flutter) {
     // Magnetic flutter terms
-    Apar_flutter = Apar - DC(Apar);
-
+    if (bout::globals::mesh->isFci()){
+      Apar_flutter = Apar;
+      set(state["fields"]["Apar_flutter"], Apar);
+    } else {
+      Apar_flutter = Apar - DC(Apar);
     // Ensure that guard cells are communicated
-    Apar.getMesh()->communicate(Apar_flutter);
+      Apar.getMesh()->communicate(Apar_flutter);
+      set(state["fields"]["Apar_flutter"], Apar_flutter);
 
-    set(state["fields"]["Apar_flutter"], Apar_flutter);
-
-#if 0
-    // Create a vector A from covariant components
-    // (A^x, A^y, A^z)
-    // Note: b = e_y / (JB)
+      #if 0
+    // Create a vector A from covariant components                                                                                                                                                         
+    // (A^x, A^y, A^z)                                                                                                                                                                                     
+    // Note: b = e_y / (JB)                                                                                                                                                                                
     const auto* coords = Apar.getCoordinates();
     Vector3D A;
     A.covariant = true;
     A.x = A.z = 0.0;
     A.y = Apar_flutter * (coords->J * coords->Bxy);
 
-    // Perturbed magnetic field vector
-    // Note: Contravariant components (dB_x, dB_y, dB_z)
+    // Perturbed magnetic field vector                                                                                                                                                                     
+    // Note: Contravariant components (dB_x, dB_y, dB_z)                                                                                                                                                   
     Vector3D delta_B = Curl(A);
 
-    // Set components of the perturbed unit vector
-    // Note: Options can't (yet) contain vectors
+    // Set components of the perturbed unit vector                                                                                                                                                         
+    // Note: Options can't (yet) contain vectors                                                                                                                                                           
     set(state["fields"]["deltab_flutter_x"], delta_B.x / coords->Bxy);
     set(state["fields"]["deltab_flutter_z"], delta_B.z / coords->Bxy);
 #endif
+    }
   }
 }
 

src/evolve_density.cxx

@@ -51,6 +51,7 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv
                            .doc("Perpendicular diffusion at low density")
                            .withDefault<bool>(false);
 
+
   pressure_floor = density_floor * (1./get<BoutReal>(alloptions["units"]["eV"]));
 
 
@@ -131,6 +132,11 @@ EvolveDensity::EvolveDensity(std::string name, Options& alloptions, Solver* solv
     .withDefault(source)
     / source_normalisation;
 
+  magnetic_flutter=
+      n_options["magnetic_flutter"]
+          .doc("Use flutter terms??")
+          .withDefault<bool>(true);
+
   disable_ddt = n_options["disable_ddt"]
     .withDefault<bool>(false);
 
@@ -299,12 +305,13 @@ void EvolveDensity::finally(const Options& state) {
     flow_ylow = 0.0;
     ddt(N) -= FV::Div_par_mod<hermes::Limiter>(N, V, fastest_wave, flow_ylow, false, dissipative);
 
-    if (state.isSection("fields") and state["fields"].isSet("Apar_flutter")) {
+    if (state.isSection("fields") and state["fields"].isSet("Apar_flutter") and magnetic_flutter) {
       // Magnetic flutter term
       const Field3D Apar_flutter = get<Field3D>(state["fields"]["Apar_flutter"]);
       // Note: Using -Apar_flutter rather than reversing sign in front,
       //       so that upwinding is handled correctly
-      ddt(N) -= Div_n_g_bxGrad_f_B_XZ(N, V, -Apar_flutter);
+      auto* coord = N.getCoordinates();
+      ddt(N) -= coord->Bxy * bracket(N*V/coord->Bxy, Apar_flutter, BRACKET_ARAKAWA) * bracket_factor;
     }
   }