Implicit drag for dust species

doc/source/modules/dust.rst

@@ -39,15 +39,45 @@ which guarantees total momentum conservation.
 
 
 
-Drag CFL condition
--------------------
-*Idefix* computes the drag terms with a time-explicit scheme. Hence, an addition CFL constraint arrises because of the drag:
+Time Integration
+----------------
+
+*Idefix* can compute the drag with a 2nd order time-explicit (default) or 1st order time-implicit scheme. Each scheme has its pros and cons. The choice
+is made by the user in the input file.
+
+Explicit scheme: drag CFL condition
++++++++++++++++++++++++++++++++++++
+
+When using the time explicit scheme, an addition CFL constraint arrises because of the drag:
 
 .. math::
 
     dt < \min(\frac{1}{\sum_i\gamma_i(\rho_i+\rho)})
 
-*Idefix* automatically adjusts the CFL to satisfy this inequality, in addition to the usual CFL condition.
+*Idefix* automatically adjusts the CFL to satisfy this inequality, in addition to the usual CFL condition. This can severly limit the time step,
+especially when the gas density is high.
+
+Implicit scheme
++++++++++++++++
+
+When the drag is applied with a 1st order implicit scheme, the drag force is applied at the end of each step. In order to avoid
+the complete inversion of the system, we follow a simplified inversion procedure, where we first update the gas momentum as:
+
+.. math::
+
+    v_g^{(n+1)}=\left(v_g^{(n)}+\sum_i\frac{\rho_i\gamma_i dt}{1+\rho \gamma_i dt}v_i^{(n)}\right)/\left(1+\sum_i\frac{\rho_i\gamma_i dt}{1+\rho\gamma_idt}\right)
+
+And then update the dust momentum as:
+
+.. math::
+
+    v_i^{(n+1)}=\left(v_i^{(n)}+\rho\gamma_i v_g^{(n+1)}dt\right)/(1+\rho\gamma_i dt)
+
+Note that the latter equation relies on the *updated* gas velocity.
+
+.. warning::
+  While the implicit scheme is more stable than the explicit one, and it does not require any additional CFL condition, it is less accurate and
+  possibly lead to inacurrate dust velocities when :math:`dt\gg (\gamma_i\rho)^{-1}`. Use it at your own risk.
 
 Dust parameters
 ---------------
@@ -66,6 +96,9 @@ The dust module can be enabled adding a block `[Dust]` in your input .ini file.
 +----------------+-------------------------+---------------------------------------------------------------------------------------------+
 | drag_feedback  | bool                    | | (optionnal) whether the gas feedback is enabled (default true).                           |
 +----------------+-------------------------+---------------------------------------------------------------------------------------------+
+| drag_implicit  | bool                    | | (optionnal) whether the drag uses a 1st order implicit method. Otherwise use the          |
+|                |                         | | 2nd order time-explicit scheme (default is false=time explicit)                           |
++----------------+-------------------------+---------------------------------------------------------------------------------------------+
 
 The drag parameter :math:`\beta_i` above sets the functional form of :math:`\gamma_i(\rho, \rho_i, c_s)` depending on the drag type:
 

src/dataBlock/evolveStage.cpp

@@ -19,6 +19,16 @@ void DataBlock::EvolveStage() {
     for(int i = 0 ; i < dust.size() ; i++) {
       dust[i]->EvolveStage(this->t,this->dt);
     }
+    // Add implicit term for dust drag
+    if(dust[0]->drag->IsImplicit()) {
+      for(int i = 0 ; i < dust.size() ; i++) {
+        dust[i]->drag->AddImplicitBackReaction(this->dt,dust[0]->drag->implicitFactor);
+      }
+      dust[0]->drag->NormalizeImplicitBackReaction(this->dt);
+      for(int i = 0 ; i < dust.size() ; i++) {
+        dust[i]->drag->AddImplicitFluidMomentum(this->dt);
+      }
+    }
   }
 
   idfx::popRegion();

src/fluid/drag.cpp

@@ -5,7 +5,9 @@
 // Licensed under CeCILL 2.1 License, see COPYING for more information
 // ***********************************************************************************
 #include "drag.hpp"
+#include <string>
 #include "physics.hpp"
+
 void Drag::AddDragForce(const real dt) {
   idfx::pushRegion("Drag::AddDragForce");
 
@@ -15,83 +17,182 @@ void Drag::AddDragForce(const real dt) {
   auto VcDust = this->VcDust;
   auto InvDt = this->InvDt;
 
-  const Type type = this->type;
-  real dragCoeff = this->dragCoeff;
   bool feedback = this->feedback;
 
-  EquationOfState eos = *(this->eos);
-
-  auto userGammai = this->gammai;
-  if(type == Type::Userdef) {
-    if(userDrag != NULL) {
-      idfx::pushRegion("Drag::UserDrag");
-      userDrag(data, dragCoeff, userGammai);
-      idfx::popRegion();
-    } else {
-      IDEFIX_ERROR("No User-defined drag function has been enrolled");
-    }
+  if(implicit) {
+    IDEFIX_ERROR("Add DragForce should not be called when drag is implicit");
   }
+
+  auto gammaDrag = this->gammaDrag;
+  gammaDrag.RefreshUserDrag(data);
+
   // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg)
   // Where gamma is computed according to the choice of drag type
   idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR],
     KOKKOS_LAMBDA (int k, int j, int i) {
-      real gamma;  // The drag coefficient
-      if(type ==  Type::Gamma) {
-        gamma = dragCoeff;
-
-      } else if(type == Type::Tau) {
-        // In this case, the coefficient is the stopping time (assumed constant)
-        gamma = 1/(dragCoeff*VcGas(RHO,k,j,i));
-      } else if(type == Type::Size) {
-        real cs;
-        // Assume a fixed size, hence for both Epstein or Stokes, gamma~1/rho_g/cs
-        // Get the sound speed
-        #if HAVE_ENERGY == 1
-          cs = std::sqrt(eos.GetGamma(VcGas(PRS,k,j,i),VcGas(RHO,k,j,i)
-                         *VcGas(PRS,k,j,i)/VcGas(RHO,k,j,i)));
-        #else
-          cs = eos.GetWaveSpeed(k,j,i);
-        #endif
-        gamma = cs/dragCoeff;
-      } else if(type == Type::Userdef) {
-        gamma = userGammai(k,j,i);
-      }
+      real gamma = gammaDrag.GetGamma(k,j,i);  // The drag coefficient
 
       real dp = dt * gamma * VcDust(RHO,k,j,i) * VcGas(RHO,k,j,i);
       for(int n = MX1 ; n < MX1+COMPONENTS ; n++) {
         real dv = VcDust(n,k,j,i) - VcGas(n,k,j,i);
         UcDust(n,k,j,i) -= dp*dv;
-        if(feedback) UcGas(n,k,j,i) += dp*dv;
+        if(feedback) {
+          UcGas(n,k,j,i) += dp*dv;
+          #if HAVE_ENERGY == 1
+            // We add back the energy dissipated for the dust which is not accounted for
+            // (since there is no energy equation for dust grains)
+
+            // TODO(GL): this should be disabled in the case of a true multifluid system where
+            // both fluids have a proper energy equation
+            UcGas(ENG,k,j,i) += dp*dv*VcDust(n,k,j,i);
+          #endif
+        } // feedback
+      }
+      // Cfl constraint
+      real idt = gamma*VcGas(RHO,k,j,i);
+      if(feedback) idt += gamma*VcDust(RHO,k,j,i);
+      InvDt(k,j,i) += idt;
+    });
+  idfx::popRegion();
+}
+//
+// $ p_g^{(n+1)}=p_g^{(n)}+\sum_i\frac{\rho_g\gamma_i dt}{1+\rho_g \gamma_i dt}p_i^{(n)} $
+// We accumulate in the array "prefactor"
+// $ \sum_i\frac{\rho_i\gamma_i dt}{1+\rho_g\gamma_i dt}
+//
+
+void Drag::AddImplicitBackReaction(const real dt, IdefixArray3D<real> preFactor) {
+  if(!feedback) {
+    // no feedback, no need for anything
+    return;
+  }
+  idfx::pushRegion("AddImplicitFluidMomentum");
+
+  auto UcGas = this->UcGas;
+  auto VcGas = this->VcGas;
+  auto VcDust = this->VcDust;
+  auto UcDust = this->UcDust;
+
+  bool isFirst = this->instanceNumber == 0;
+
+
+  if(!implicit) {
+    IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit");
+  }
+
+  auto gammaDrag = this->gammaDrag;
+  gammaDrag.RefreshUserDrag(data);
+
+  // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg)
+  // Where gamma is computed according to the choice of drag type
+  idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR],
+    KOKKOS_LAMBDA (int k, int j, int i) {
+      real gamma = gammaDrag.GetGamma(k,j,i);  // The drag coefficient
+
+
+      const real factor = UcDust(RHO,k,j,i)*gamma*dt/(1+UcGas(RHO,k,j,i)*gamma*dt);
+      if(isFirst) {
+        preFactor(k,j,i) = factor;
+      } else {
+        preFactor(k,j,i) += factor;
+      }
+
+      for(int n = MX1 ; n < MX1+COMPONENTS ; n++) {
+        UcGas(n,k,j,i) +=  dt * gamma * UcGas(RHO,k,j,i) * UcDust(n,k,j,i) /
+                                                                   (1 + UcGas(RHO,k,j,i)*dt*gamma);
+      }
+    });
+  idfx::popRegion();
+}
+
+void Drag::NormalizeImplicitBackReaction(const real dt) {
+  if(!feedback) {
+    // no feedback, no need for anything
+    return;
+  }
+  idfx::pushRegion("AddImplicitFluidMomentum");
+
+  auto UcGas = this->UcGas;
+  auto preFactor = this->implicitFactor;
+
+  if(!implicit) {
+    IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit");
+  }
+
+  // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg)
+  // Where gamma is computed according to the choice of drag type
+  idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR],
+    KOKKOS_LAMBDA (int k, int j, int i) {
+      const real factor = 1+preFactor(k,j,i);
+      for(int n = MX1 ; n < MX1+COMPONENTS ; n++) {
+        UcGas(n,k,j,i) /= factor;
+      }
+    });
+  idfx::popRegion();
+}
+
+void Drag::AddImplicitFluidMomentum(const real dt) {
+  idfx::pushRegion("AddImplicitFluidMomentum");
+
+  auto UcGas = this->UcGas;
+  auto VcGas = this->VcGas;
+  auto UcDust = this->UcDust;
+  auto VcDust = this->VcDust;
+  auto InvDt = this->InvDt;
+
+  bool feedback = this->feedback;
+
+  if(!implicit) {
+    IDEFIX_ERROR("AddImplicitGasMomentum should not be called when drag is explicit");
+  }
+
+  auto gammaDrag = this->gammaDrag;
+  if(!feedback) gammaDrag.RefreshUserDrag(data);
+
+  // Compute a drag force fd = - gamma*rhod*rhog*(vd-vg)
+  // Where gamma is computed according to the choice of drag type
+  idefix_for("DragForce",0,data->np_tot[KDIR],0,data->np_tot[JDIR],0,data->np_tot[IDIR],
+    KOKKOS_LAMBDA (int k, int j, int i) {
+      real gamma = gammaDrag.GetGamma(k,j,i);  // The drag coefficient
+
+      for(int n = MX1 ; n < MX1+COMPONENTS ; n++) {
+        real oldUc = UcDust(n,k,j,i);
+        UcDust(n,k,j,i) = (oldUc + dt * gamma * UcDust(RHO,k,j,i) * UcGas(n,k,j,i)) /
+                          (1 + UcGas(RHO,k,j,i)*dt*gamma);
+
         #if HAVE_ENERGY == 1
+          real dp = UcDust(n,k,j,i) - oldUc;
+
           // We add back the energy dissipated for the dust which is not accounted for
           // (since there is no energy equation for dust grains)
 
           // TODO(GL): this should be disabled in the case of a true multifluid system where
           // both fluids have a proper energy equation
-          UcGas(ENG,k,j,i) += dp*dv*VcDust(n,k,j,i);
+          if(feedback) UcGas(ENG,k,j,i) -= dp*VcDust(n,k,j,i);
         #endif
       }
-      // Cfl constraint
-      real idt = gamma*VcGas(RHO,k,j,i);
-      if(feedback) idt += gamma*VcDust(RHO,k,j,i);
-      InvDt(k,j,i) += idt;
     });
   idfx::popRegion();
 }
 
 void Drag::ShowConfig() {
   idfx::cout << "Drag: Using ";
-  switch(type) {
-    case Type::Gamma:
+  if(implicit) {
+    idfx::cout << " IMPLICIT ";
+  } else {
+    idfx::cout << " EXPLICIT ";
+  }
+  switch(gammaDrag.type) {
+    case GammaDrag::Type::Gamma:
       idfx::cout << "constant gamma";
       break;
-    case Type::Tau:
+    case GammaDrag::Type::Tau:
       idfx::cout << "constant stopping time";
       break;
-    case Type::Size:
+    case GammaDrag::Type::Size:
       idfx::cout << "constant dust size";
       break;
-    case Type::Userdef:
+    case GammaDrag::Type::Userdef:
       idfx::cout << "user-defined";
       break;
   }
@@ -105,8 +206,60 @@ void Drag::ShowConfig() {
 }
 
 void Drag::EnrollUserDrag(UserDefDragFunc func) {
+  gammaDrag.EnrollUserDrag(func);
+}
+
+////////////////////////////////////////////
+// GammaDrag function definitions
+////////////////////////////////////////////
+
+void GammaDrag::RefreshUserDrag(DataBlock *data) {
+  if(type == Type::Userdef) {
+    if(userDrag != NULL) {
+      idfx::pushRegion("GammaDrag::UserDrag");
+        userDrag(data, instanceNumber, dragCoeff, gammai);
+      idfx::popRegion();
+    } else {
+      IDEFIX_ERROR("No User-defined drag function has been enrolled");
+    }
+  }
+}
+
+void GammaDrag::EnrollUserDrag(UserDefDragFunc func) {
   if(type != Type::Userdef) {
     IDEFIX_ERROR("User-defined drag function requires drag entry to be set to \"userdef\"");
   }
   this->userDrag = func;
 }
+
+GammaDrag::GammaDrag(Input &input, std::string BlockName, int instanceNumber, DataBlock *data) {
+  if(input.CheckEntry(BlockName,"drag")>=0) {
+    std::string dragType = input.Get<std::string>(BlockName,"drag",0);
+    if(dragType.compare("gamma") == 0) {
+      this->type = Type::Gamma;
+    } else if(dragType.compare("tau") == 0) {
+      this->type = Type::Tau;
+      this->VcGas = data->hydro->Vc;
+    } else if(dragType.compare("size") == 0) {
+      this->type = Type::Size;
+      this->eos = *(data->hydro->eos.get());
+      this->VcGas = data->hydro->Vc;
+    } else if(dragType.compare("userdef") == 0) {
+      this->type = Type::Userdef;
+      this->gammai = IdefixArray3D<real>("UserDrag",
+                                  data->np_tot[KDIR],
+                                  data->np_tot[JDIR],
+                                  data->np_tot[IDIR]);
+    } else {
+      std::stringstream msg;
+      msg << "Unknown drag type \"" <<  dragType
+          << "\" in your input file." << std::endl
+          << "Allowed values are: gamma, tau, epstein, stokes, userdef." << std::endl;
+
+      IDEFIX_ERROR(msg);
+    }
+  }
+  // Fetch the drag coefficient for the current specie.
+  this->dragCoeff = input.Get<real>(BlockName,"drag",instanceNumber+1);
+  this->instanceNumber = instanceNumber;
+}