Fci: Neutral model for curvilinear coorinates

CMakeLists.txt

@@ -58,6 +58,7 @@ set(HERMES_SOURCES
     src/div_ops.cxx
     src/fci_gridcheck.cxx
     src/neutral_mixed.cxx
+    src/neutral_full_velocity_curv.cxx
     src/full_velocity.cxx
     src/electromagnetic.cxx
     src/electron_force_balance.cxx
@@ -149,6 +150,7 @@ set(HERMES_SOURCES
     include/ionisation.hxx
     include/neutral_boundary.hxx
     include/neutral_mixed.hxx
+    include/neutral_full_velocity_curv.hxx
     include/neutral_parallel_diffusion.hxx
     include/solkit_neutral_parallel_diffusion.hxx
     include/noflow_boundary.hxx
@@ -298,6 +300,10 @@ if(HERMES_TESTS)
 	DEPENDS setup-test)
     endif()
   endfunction()
+  hermes_add_integrated_test(2D-neutral-full-density-fci
+    REQUIRES BOUT_ENABLE_METRIC_3D)
+  hermes_add_integrated_test(2D-neutral-full-momentum-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

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