fix for rank discontinuity in CG

src/EquationSystems/DoubleDiffusiveField.hpp

@@ -52,7 +52,7 @@ class DoubleDiffusiveField : public TokamakSystem
 
     void load_params() override;
 
-    bool v_PostIntegrate(int step);
+    bool v_PostIntegrate(int step) override;
 
     void v_ExtraFldOutput(std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
                           std::vector<std::string> &variables) override;

src/EquationSystems/SingleDiffusiveField.cpp

@@ -38,6 +38,17 @@ SingleDiffusiveField::SingleDiffusiveField(
 void SingleDiffusiveField::v_InitObject(bool DeclareFields)
 {
     TokamakSystem::v_InitObject(DeclareFields);
+
+    m_session->MatchSolverInfo("SpectralVanishingViscosity", "True",
+                               m_useSpecVanVisc, false);
+    m_session->LoadParameter("epsilon", m_epsilon, 1.0);
+    if (m_useSpecVanVisc)
+    {
+        m_session->LoadParameter("SVVCutoffRatio", m_sVVCutoffRatio, 0.75);
+        m_session->LoadParameter("SVVDiffCoeff", m_sVVDiffCoeff, 0.1);
+        m_factors[StdRegions::eFactorSVVCutoffRatio] = m_sVVCutoffRatio;
+        m_factors[StdRegions::eFactorSVVDiffCoeff] = m_sVVDiffCoeff / m_epsilon;
+    }
     CalcDiffTensor();
 
     // Setup diffusion object
@@ -94,23 +105,46 @@ void SingleDiffusiveField::ImplicitTimeIntCG(
     int npoints                          = m_fields[0]->GetNpoints();
     m_factors[StdRegions::eFactorLambda] = 1.0 / lambda / m_epsilon;
 
-    // if (m_useSpecVanVisc)
-    // {
-    //     m_factors[StdRegions::eFactorSVVCutoffRatio] = m_sVVCutoffRatio;
-    //     m_factors[StdRegions::eFactorSVVDiffCoeff] = m_sVVDiffCoeff /
-    //     m_epsilon;
-    // }
+    if (m_useSpecVanVisc)
+    {
+        m_factors[StdRegions::eFactorSVVCutoffRatio] = m_sVVCutoffRatio;
+        m_factors[StdRegions::eFactorSVVDiffCoeff] = m_sVVDiffCoeff / m_epsilon;
+    }
 
     // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
     // inarray = input: \hat{rhs} -> output: \hat{Y}
     // outarray = output: nabla^2 \hat{Y}
     // where \hat = modal coeffs
+    if (m_intScheme->GetIntegrationSchemeType() == LibUtilities::eImplicit)
+    {
+        for (int i = 0; i < nvariables; ++i)
+        {
+            Vmath::Zero(npoints, outarray[i], 1);
+        }
+
+        for (auto &x : m_forcing)
+        {
+            x->Apply(m_fields, inarray, outarray, time);
+        }
+    }
     CalcDiffTensor();
     for (int i = 0; i < nvariables; ++i)
     {
-        // Multiply 1.0/timestep/lambda
-        Vmath::Smul(npoints, -m_factors[StdRegions::eFactorLambda], inarray[i],
-                    1, outarray[i], 1);
+        if (m_intScheme->GetIntegrationSchemeType() == LibUtilities::eImplicit)
+        {
+            // Multiply forcing term by -1 for definition of HelmSolve function
+            Vmath::Smul(npoints, -1.0, outarray[i], 1, outarray[i], 1);
+
+            // Multiply 1.0/timestep/lambda
+            Vmath::Svtvp(npoints, -m_factors[StdRegions::eFactorLambda],
+                         inarray[i], 1, outarray[i], 1, outarray[i], 1);
+        }
+        else
+        {
+            // Multiply 1.0/timestep/lambda
+            Vmath::Smul(npoints, -m_factors[StdRegions::eFactorLambda],
+                        inarray[i], 1, outarray[i], 1);
+        }
 
         // Solve a system of equations with Helmholtz solver
         m_fields[i]->HelmSolve(outarray[i], m_fields[i]->UpdateCoeffs(),
@@ -171,18 +205,30 @@ void SingleDiffusiveField::DoOdeRhs(
     const Array<OneD, const Array<OneD, NekDouble>> &in_arr,
     Array<OneD, Array<OneD, NekDouble>> &out_arr, const NekDouble time)
 {
-    CalcDiffTensor();
-    size_t nvariables = in_arr.size();
-    m_diffusion->Diffuse(nvariables, m_fields, in_arr, out_arr);
+    if (m_explicitDiffusion)
+    {
+        CalcDiffTensor();
+        size_t nvariables = in_arr.size();
+        m_diffusion->Diffuse(nvariables, m_fields, in_arr, out_arr);
+    }
+    else
+    {
+        // RHS should be set to zero.
+        for (int i = 0; i < out_arr.size(); ++i)
+        {
+            Vmath::Zero(out_arr[i].size(), &out_arr[i][0], 1);
+        }
+    }
 
     if (this->particles_enabled)
     {
-        for (auto &s : particle_sys->get_species())
+        for (int i = 0; i < this->particle_sys->get_species().size(); ++i)
         {
             Vmath::Vadd(out_arr[0].size(), out_arr[0], 1,
-                        this->src_fields[s.first]->GetPhys(), 1, out_arr[0], 1);
+                        this->src_fields[i]->GetPhys(), 1, out_arr[0], 1);
         }
     }
+
     // Add forcing terms
     for (auto &x : m_forcing)
     {
@@ -230,11 +276,13 @@ void SingleDiffusiveField::v_ExtraFldOutput(
     TokamakSystem::v_ExtraFldOutput(fieldcoeffs, variables);
     const int nPhys   = m_fields[0]->GetNpoints();
     const int nCoeffs = m_fields[0]->GetNcoeffs();
-    for (auto s : this->particle_sys->get_species())
+    int i             = 0;
+    for (auto &[k, v] : this->particle_sys->get_species())
     {
-        variables.push_back(s.second.name + "_SOURCE_DENSITY");
+        variables.push_back(v.name + "_SOURCE_DENSITY");
         Array<OneD, NekDouble> SrcFwd(nCoeffs);
-        m_fields[0]->FwdTransLocalElmt(src_fields[s.first]->GetPhys(), SrcFwd);
+        m_fields[0]->FwdTransLocalElmt(this->src_fields[i++]->GetPhys(),
+                                       SrcFwd);
         fieldcoeffs.push_back(SrcFwd);
     }
 }

src/EquationSystems/SingleDiffusiveField.hpp

@@ -57,7 +57,11 @@ class SingleDiffusiveField : public TokamakSystem
 
     // For Diffusion
     StdRegions::ConstFactorMap m_factors;
-    /// Weight for average calculation of diffusion term
+    NekDouble m_epsilon;
+    bool m_useSpecVanVisc;
+    NekDouble
+        m_sVVCutoffRatio; // Cut-off ratio from which to start decaying modes
+    NekDouble m_sVVDiffCoeff; // Diffusion coefficient of SVV modes
 
     NekDouble m_k_B;
     NekDouble k_par;