Improvements and bug fixes for LTE+Mixing-Length model

src/M2ulPhyS.cpp

@@ -1640,6 +1640,8 @@ void M2ulPhyS::initSolutionAndVisualizationVectors() {
   if (spaceVaryViscMult != NULL) paraviewColl->RegisterField("viscMult", spaceVaryViscMult);
 
   if (distance_ != NULL) paraviewColl->RegisterField("distance", distance_);
+  if (plasma_conductivity_ != NULL) paraviewColl->RegisterField("sigma", plasma_conductivity_);
+  if (joule_heating_ != NULL) paraviewColl->RegisterField("Sjoule", joule_heating_);
 
   paraviewColl->SetOwnData(true);
   // paraviewColl->Save();
@@ -1700,7 +1702,10 @@ void M2ulPhyS::projectInitialSolution() {
 
   // update plasma electrical conductivity
   if (tpsP->isFlowEMCoupled()) {
-    mixture->SetConstantPlasmaConductivity(plasma_conductivity_, Up);
+    ParGridFunction *coordsDof = new ParGridFunction(dfes);
+    mesh->GetNodes(*coordsDof);
+    mixture->SetConstantPlasmaConductivity(plasma_conductivity_, Up, coordsDof);
+    delete coordsDof;
   }
 
   if (config.GetRestartCycle() == 0 && !loadFromAuxSol) {
@@ -2394,6 +2399,7 @@ void M2ulPhyS::parseSolverOptions2() {
     parsePlasmaModels();
   } else {
     // parse options for other plasma presets.
+    config.const_plasma_conductivity_ = 50.0;
   }
 
   // species list.
@@ -2471,6 +2477,7 @@ void M2ulPhyS::parseFlowOptions() {
     tpsP->getInput("flow/mixing-length/max-mixing-length", config.mix_length_trans_input_.max_mixing_length_, 0.0);
     tpsP->getInput("flow/mixing-length/Pr_ratio", config.mix_length_trans_input_.Prt_, 1.0);
     tpsP->getInput("flow/mixing-length/Let", config.mix_length_trans_input_.Let_, 1.0);
+    tpsP->getInput("flow/mixing-length/bulk-multiplier", config.mix_length_trans_input_.bulk_multiplier_, 0.0);
   }
 
   assert(config.solOrder > 0);
@@ -3991,3 +3998,67 @@ void M2ulPhyS::updateVisualizationVariables() {
     }  // if (!isDryAir)
   }    // for (int n = 0; n < ndofs; n++)
 }
+
+void M2ulPhyS::evaluatePlasmaConductivityGF() {
+  assert(plasma_conductivity_ != NULL);
+  double *d_pc = plasma_conductivity_->Write();
+
+  const double *d_Up = Up->Read();
+  const double *d_U = U->Read();
+  const double *d_gradUp = gradUp->Read();
+
+  const double *d_distance = NULL;
+  if (distance_ != NULL) {
+    d_distance = distance_->Read();
+  }
+
+  TransportProperties *d_transport = transportPtr;
+
+  const int nnodes = vfes->GetNDofs();
+  const int _dim = dim;
+  const int _nvel = nvel;
+  const int _num_equation = num_equation;
+  const int _numSpecies = numSpecies;
+
+#if defined(_GPU_)
+  MFEM_FORALL(n, nnodes, {
+    double upn[gpudata::MAXEQUATIONS];
+    double Un[gpudata::MAXEQUATIONS];
+    double gradUpn[gpudata::MAXEQUATIONS * gpudata::MAXDIM];
+#else
+  double upn[gpudata::MAXEQUATIONS];
+  double Un[gpudata::MAXEQUATIONS];
+  double gradUpn[gpudata::MAXEQUATIONS * gpudata::MAXDIM];
+  for (int n = 0; n < nnodes; n++) {
+#endif
+    for (int eq = 0; eq < _num_equation; eq++) {
+      upn[eq] = d_Up[n + eq * nnodes];
+      Un[eq] = d_U[n + eq * nnodes];
+      for (int d = 0; d < _dim; d++)
+        gradUpn[eq + d * _num_equation] = d_gradUp[n + eq * nnodes + d * _num_equation * nnodes];
+    }
+    // TODO(kevin): update E-field with EM coupling.
+    // E-field can have azimuthal component.
+    double Efield[gpudata::MAXDIM];
+    for (int v = 0; v < _nvel; v++) Efield[v] = 0.0;
+
+    double globalTransport[gpudata::MAXSPECIES];
+    double speciesTransport[gpudata::MAXSPECIES * SpeciesTrns::NUM_SPECIES_COEFFS];
+    // NOTE: diffusion has nvel components, as E-field can have azimuthal component.
+    double diffusionVelocity[gpudata::MAXSPECIES * gpudata::MAXDIM];
+    for (int v = 0; v < _nvel; v++)
+      for (int sp = 0; sp < _numSpecies; sp++) diffusionVelocity[sp + v * _numSpecies] = 0.0;
+    double ns[gpudata::MAXSPECIES];
+
+    double dist = 0.0;
+    if (d_distance != NULL) dist = d_distance[n];
+    d_transport->ComputeSourceTransportProperties(Un, upn, gradUpn, Efield, dist, globalTransport, speciesTransport,
+                                                 diffusionVelocity, ns);
+
+    d_pc[n] = globalTransport[SrcTrns::ELECTRIC_CONDUCTIVITY];
+#if defined(_GPU_)
+  });  // MFEM_FORALL(n, nnodes, {
+#else
+  }  // for (int n = 0; n < nnodes; n++)
+#endif
+}

src/M2ulPhyS.hpp

@@ -430,6 +430,15 @@ class M2ulPhyS : public TPS::Solver {
 
   static int Check_NaN_GPU(ParGridFunction *U, int lengthU, Array<int> &loc_print);
 
+  void setConstantPlasmaConductivityGF() {
+    ParGridFunction *coordsDof = new ParGridFunction(dfes);
+    mesh->GetNodes(*coordsDof);
+    mixture->SetConstantPlasmaConductivity(plasma_conductivity_, Up, coordsDof);
+    delete coordsDof;
+  }
+
+  void evaluatePlasmaConductivityGF();
+
   // Exit code access
   void SetStatus(int code) {
     exit_status_ = code;

src/cycle_avg_joule_coupling.cpp

@@ -50,6 +50,7 @@ CycleAvgJouleCoupling::CycleAvgJouleCoupling(string &inputFileName, TPS::Tps *tp
   tps->getInput("cycle-avg-joule-coupled/axisymmetric", axisym, false);
   tps->getInput("cycle-avg-joule-coupled/input-power", input_power_, -1.);
   tps->getInput("cycle-avg-joule-coupled/initial-input-power", initial_input_power_, -1.);
+  tps->getInput("cycle-avg-joule-coupled/fixed-conductivity", fixed_conductivity_, false);
 
   if (axisym) {
     qmsa_solver_ = new QuasiMagnetostaticSolverAxiSym(em_opt_, tps);
@@ -307,17 +308,27 @@ void CycleAvgJouleCoupling::solveBegin() {
 void CycleAvgJouleCoupling::solveStep() {
   // Run the em solver when it is due
   if (current_iter_ % solve_em_every_n_ == 0) {
+    // update the power if necessary
     double delta_power = 0;
     if (input_power_ > 0) {
       delta_power = (input_power_ - initial_input_power_) * static_cast<double>(solve_em_every_n_) /
                     static_cast<double>(max_iters_);
     }
+
+    // evaluate electric conductivity and interpolate it to EM mesh
+    if (!fixed_conductivity_) flow_solver_->evaluatePlasmaConductivityGF();
     interpConductivityFromFlowToEM();
+
+    // solve the EM system
     qmsa_solver_->solve();
+
+    // report the "raw" Joule heating
     const double tot_jh = qmsa_solver_->totalJouleHeating();
     if (rank0_) {
       grvy_printf(GRVY_INFO, "The total input Joule heating = %.6e\n", tot_jh);
     }
+
+    // scale the Joule heating (if we are controlling the power input)
     if (input_power_ > 0) {
       const double target_power = initial_input_power_ + (current_iter_ / solve_em_every_n_ + 1) * delta_power;
       const double ratio = target_power / tot_jh;
@@ -327,6 +338,8 @@ void CycleAvgJouleCoupling::solveStep() {
         grvy_printf(GRVY_INFO, "The total input Joule heating after scaling = %.6e\n", upd_jh);
       }
     }
+
+    // interpolate the Joule heating to the flow mesh
     interpJouleHeatingFromEMToFlow();
   }
   // Run a step of the flow solver

src/cycle_avg_joule_coupling.hpp

@@ -61,6 +61,8 @@ class CycleAvgJouleCoupling : public TPS::Solver {
   int solve_em_every_n_;
   //! Report timings every timing_freq_ steps
   int timing_freq_;
+  //! Flag to specify the electrical conductivity should not be updated
+  bool fixed_conductivity_;
 
 #ifdef HAVE_GSLIB
   FindPointsGSLIB *interp_flow_to_em_;

src/dataStructures.hpp

@@ -473,6 +473,7 @@ struct mixingLengthTransportData {
   double max_mixing_length_;  // user-specifed mixing length
   double Prt_;                // eddy Prandtl number
   double Let_;                // eddy Lewis number
+  double bulk_multiplier_;    // bulk eddy visc = bulk_multiplier_ * mut
 };
 
 struct collisionInputs {

src/equation_of_state.cpp

@@ -44,19 +44,37 @@ MFEM_HOST_DEVICE GasMixture::GasMixture(WorkingFluid f, int _dim, int nvel, doub
   const_plasma_conductivity_ = pc;
 }
 
-void GasMixture::SetConstantPlasmaConductivity(ParGridFunction *pc, const ParGridFunction *Up) {
+void GasMixture::SetConstantPlasmaConductivity(ParGridFunction *pc, const ParGridFunction *Up,
+                                               const ParGridFunction *coords) {
   // quick return if pc is NULL (nothing to set)
   if (pc == NULL) return;
 
-  // otherwise, set plasma conductivity based on temperature
+  // otherwise, set plasma conductivity
   double *plasma_conductivity_gf = pc->HostWrite();
   const double *UpData = Up->HostRead();
 
+  // To a constant
   const int nnode = pc->FESpace()->GetNDofs();
-
   for (int n = 0; n < nnode; n++) {
     plasma_conductivity_gf[n] = const_plasma_conductivity_;
   }
+
+  // To use a spatially varying conductivity, uncomment the code
+  // below, and put in the function of space you wish to use.  Note
+  // that both the spatial coordinates and the primitive variables are
+  // available to construct this function.
+  //
+  // if (coords != NULL) {
+  //   for (int n = 0; n < nnode; n++) {
+  //     const double r0 = 0.005; // 5mm
+  //     const double x = (*coords)[n + 0 * nnode];
+  //     plasma_conductivity_gf[n] = 10.0 * const_plasma_conductivity_ * std::exp(-0.5 * (x / r0) * (x / r0));
+  //   }
+  // } else {
+  //   for (int n = 0; n < nnode; n++) {
+  //     plasma_conductivity_gf[n] = const_plasma_conductivity_;
+  //   }
+  // }
 }
 
 void GasMixture::UpdatePressureGridFunction(ParGridFunction *press, const ParGridFunction *Up) {

src/equation_of_state.hpp

@@ -267,7 +267,7 @@ class GasMixture {
     mfem_error("computeNumberDensities not implemented");
   }
 
-  void SetConstantPlasmaConductivity(ParGridFunction *pc, const ParGridFunction *Up);
+  void SetConstantPlasmaConductivity(ParGridFunction *pc, const ParGridFunction *Up, const ParGridFunction *coords);
 
   virtual void UpdatePressureGridFunction(ParGridFunction *press, const ParGridFunction *Up);
 

src/face_integrator.cpp

@@ -180,7 +180,7 @@ void FaceIntegrator::getDistanceDofs(FaceElementTransformations &Tr, const Finit
     dist2.SetSize(dofs2.Size());
     distance_->FaceNbrData().GetSubVector(dofs2, dist2);
   } else {
-    pfes->GetElementVDofs(no2, vdofs2);
+    pfes->GetElementVDofs(no2, dofs2);
     dist2.SetSize(dofs2.Size());
     distance_->GetSubVector(dofs2, dist2);
   }

src/forcing_terms.cpp

@@ -466,7 +466,9 @@ void JouleHeating::updateTerms(Vector &in) {
       const double heating = jh[h_index];
       // std::cout << "heating = " << heating << std::endl;
       const int e_index = h_index + (nvel + 1) * dof;
-      data[e_index] += heating;
+      if (heating > 0.) {
+        data[e_index] += heating;
+      }
     }
 
     // Add Joule heating to electron energy (assumes ion Joule heating is negligible)
@@ -475,7 +477,9 @@ void JouleHeating::updateTerms(Vector &in) {
         const int h_index = nodes[n];
         const double heating = jh[h_index];
         const int ee_index = h_index + (num_equation - 1) * dof;
-        data[ee_index] += heating;
+        if (heating > 0.) {
+          data[ee_index] += heating;
+        }
       }
     }
   }

src/io.cpp

@@ -359,13 +359,13 @@ void M2ulPhyS::restart_files_hdf5(string mode, string inputFileName) {
     // that code is duplicated here.
     // TODO(kevin): use mixture comptue primitive.
     double *dataUp = Up->HostReadWrite();
-    double *x = U->HostReadWrite();
+    const double *x = U->HostRead();
     for (int i = 0; i < vfes->GetNDofs(); i++) {
-      Vector iState(num_equation);
-      Vector conservedState(num_equation);
-      for (int eq = 0; eq < num_equation; eq++) iState[eq] = x[i + eq * vfes->GetNDofs()];
-      mixture->GetConservativesFromPrimitives(iState, conservedState);
-      for (int eq = 0; eq < num_equation; eq++) dataUp[i + eq * vfes->GetNDofs()] = conservedState[eq];
+      Vector conserved(num_equation);
+      Vector primitive(num_equation);
+      for (int eq = 0; eq < num_equation; eq++) conserved[eq] = x[i + eq * vfes->GetNDofs()];
+      mixture->GetPrimitivesFromConservatives(conserved, primitive);
+      for (int eq = 0; eq < num_equation; eq++) dataUp[i + eq * vfes->GetNDofs()] = primitive[eq];
     }
 
     // clean up aux data

src/lte_mixture.cpp

@@ -38,7 +38,7 @@
 #ifndef _GPU_  // this class only available for CPU currently
 
 LteMixture::LteMixture(RunConfiguration &_runfile, int _dim, int nvel)
-    : GasMixture(_runfile.lteMixtureInput.f, _dim, nvel) {
+    : GasMixture(_runfile.lteMixtureInput.f, _dim, nvel, _runfile.const_plasma_conductivity_) {
   numSpecies = 1;
   ambipolar = false;
   twoTemperature_ = false;

src/lte_transport_properties.cpp

@@ -91,6 +91,7 @@ void LteTransport::ComputeFluxTransportProperties(const double *state, const dou
   transportBuffer[FluxTrns::VISCOSITY] = mu_table_->eval(T, rho);
   transportBuffer[FluxTrns::BULK_VISCOSITY] = 0.0;  // bulk_visc_mult * viscosity;
   transportBuffer[FluxTrns::HEAVY_THERMAL_CONDUCTIVITY] = kappa_table_->eval(T, rho);
+  transportBuffer[FluxTrns::ELECTRON_THERMAL_CONDUCTIVITY] = 0.0;  // electron conductivity already accounted for
 }
 
 
@@ -111,14 +112,27 @@ void LteTransport::ComputeSourceTransportProperties(const Vector &state, const V
   for (int v = 0; v < nvel_; v++)
     for (int sp = 0; sp < numSpecies; sp++) diffusionVelocity(sp, v) = 0.0;
 
-  // TODO(trevilo): Compute conductivity!
-  globalTransport[SrcTrns::ELECTRIC_CONDUCTIVITY] = 0.;
+  const double rho = Up[0];
+  const double T = Up[1 + nvel_];
+
+  double sigma = sigma_table_->eval(T, rho);
+  if (sigma < 1.0) {
+    sigma = 1.0;
+  }
+  globalTransport[SrcTrns::ELECTRIC_CONDUCTIVITY] = sigma;
 }
 
 void LteTransport::ComputeSourceTransportProperties(const double *state, const double *Up, const double *gradUp,
                                                     const double *Efield, double distance, double *globalTransport,
                                                     double *speciesTransport, double *diffusionVelocity, double *n_sp) {
-  globalTransport[SrcTrns::ELECTRIC_CONDUCTIVITY] = 0.;
+  const double rho = Up[0];
+  const double T = Up[1 + nvel_];
+  double sigma = sigma_table_->eval(T, rho);
+
+  if (sigma < 1.0) {
+    sigma = 1.0;
+  }
+  globalTransport[SrcTrns::ELECTRIC_CONDUCTIVITY] = sigma;
 }
 
 void LteTransport::GetViscosities(const double *conserved, const double *primitive, double *visc) {