Implement Darcy Physics

Code/Source/solver/CMakeLists.txt

@@ -212,6 +212,7 @@ set(CSRCS
   cmm.h cmm.cpp
   consts.h consts.cpp
   contact.h contact.cpp
+  darcy.h darcy.cpp
   distribute.h distribute.cpp
   eq_assem.h eq_assem.cpp
   fluid.h fluid.cpp

Code/Source/solver/ComMod.h

@@ -387,7 +387,7 @@ class dmnType
 
     // General physical properties such as density, elastic modulus...
     // FIX davep double prop[maxNProp] ;
-    std::map<consts::PhysicalProperyType,double> prop;
+    std::map<consts::PhysicalPropertyType,double> prop;
     //double prop[consts::maxNProp];
 
     // Electrophysiology model

Code/Source/solver/Parameters.cpp

@@ -2052,6 +2052,11 @@ DomainParameters::DomainParameters() {
 
   set_parameter("Penalty_parameter", 0.0, !required, penalty_parameter);
   set_parameter("Poisson_ratio", 0.3, !required, poisson_ratio);
+
+  set_parameter("Darcy_permeability", 1e-15, !required, darcy_permeability);
+  set_parameter("Darcy_compressibility", 0.0, !required,
+                darcy_compressibility);
+  set_parameter("Darcy_fluid_viscosity", 1.0, !required, darcy_fluid_viscosity);
 
   set_parameter("Relative_tolerance", 1e-4, !required, relative_tolerance);
   set_parameter("Shell_thickness", 0.0, !required, shell_thickness);
@@ -2060,8 +2065,8 @@ DomainParameters::DomainParameters() {
   set_parameter("Time_step_for_integration", 0.0, !required,
                 time_step_for_integration);
 
-  set_parameter("Inverse_darcy_permeability", 0.0, !required,
-                inverse_darcy_permeability);
+  set_parameter("Brinkman_inverse_permeability", 0.0, !required,
+                brinkman_inverse_permeability);
 
   // Ionic model parameters.
   IonicModelFactory::visit(

Code/Source/solver/Parameters.h

@@ -1661,8 +1661,13 @@ class DomainParameters : public ParameterLists
     Parameter<double> source_term;
     Parameter<double> time_step_for_integration;
 
-    // Inverse of Darcy permeability. Default value of 0.0 for Navier-Stokes and non-zero for Navier-Stokes-Brinkman
-    Parameter<double> inverse_darcy_permeability;
+    Parameter<double> darcy_permeability;
+    Parameter<double> darcy_compressibility;
+    Parameter<double> darcy_fluid_viscosity;
+
+    // Inverse permeability K^{-1} used in the Brinkman drag term
+    // mu K^{-1} u. A value of zero disables Brinkman drag.
+    Parameter<double> brinkman_inverse_permeability;
 };
 
 /// @brief The RemesherParameters class stores parameters for the 
@@ -1783,9 +1788,6 @@ class EquationParameters : public ParameterLists
     // and only then is the mesh equation solved.
     Parameter<bool> explicit_geometric_coupling;
 
-    // Inverse of Darcy permeability. Default value of 0.0 for Navier-Stokes and non-zero for Navier-Stokes-Brinkman
-    Parameter<double> inverse_darcy_permeability;
-
     // Sub-element parameters.
     //
     std::vector<BodyForceParameters*> body_forces;

Code/Source/solver/README.md

@@ -178,8 +178,8 @@ C++ functions are defined within a `namespace` defined for each Fortran file. Fo
              - [<b> fs::get_thood_fs(com_mod, fs, lM, vmsStab, 1) </b>](#)
              - [<b> nn::gnn(fs[1].eNoN, nsd, nsd, Nx, xql, Nqx, Jac, ksix) </b>](#)
              - [<b> nn::gn_nxx(l, fs[0].eNoN, nsd, nsd, Nx, Nxx, xwl, Nwx, Nwxx) </b>](#)
-             - [<b> fluid_3d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability)</b>](#) - If nsd=3
-             - [<b> fluid_2d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability)</b>](#) - If nsd=2
+             - [<b> fluid_3d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, brinkman_inverse_permeability)</b>](#) - If nsd=3
+             - [<b> fluid_2d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, brinkman_inverse_permeability)</b>](#) - If nsd=2
              - [<b> trilinos_doassem_(const_cast<int&>(eNoN), ptr.data(), lK.data(), lR.data())</b>](#) - If using Trilinos
              - [<b> lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR)</b>](#do_assem) - If not using Trilinos
          - [<b> set_bc::set_bc_neu(com_mod, cm_mod, Yg, Dg) </b>](#set_bc_neu)
@@ -1492,9 +1492,9 @@ strongly or weakly.
 
  - `nn::gn_nxx(l, fs[0].eNoN, nsd, nsd, Nx, Nxx, xwl, Nwx, Nwxx)`
 
- - `fluid_3d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability)` - If nsd=3
+ - `fluid_3d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, brinkman_inverse_permeability)` - If nsd=3
 
- - `fluid_2d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability)` - If nsd=2
+ - `fluid_2d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, brinkman_inverse_permeability)` - If nsd=2
 
  - `trilinos_doassem_(const_cast<int&>(eNoN), ptr.data(), lK.data(), lR.data())` - If using Trilinos
 

Code/Source/solver/cmm.cpp

@@ -37,10 +37,10 @@ void cmm_3d(ComMod& com_mod, const int eNoN, const double w, const Vector<double
   const double ctM = 1.0;
   const double ctC = 36.0;
 
-  double rho = dmn.prop.at(PhysicalProperyType::fluid_density);
-  Vector<double> f({dmn.prop.at(PhysicalProperyType::f_x), 
-                    dmn.prop.at(PhysicalProperyType::f_y), 
-                    dmn.prop.at(PhysicalProperyType::f_z)});
+  double rho = dmn.prop.at(PhysicalPropertyType::fluid_density);
+  Vector<double> f({dmn.prop.at(PhysicalPropertyType::f_x), 
+                    dmn.prop.at(PhysicalPropertyType::f_y), 
+                    dmn.prop.at(PhysicalPropertyType::f_z)});
 
   double T1 = eq.af * eq.gam * dt;
   double amd = eq.am/T1;
@@ -429,10 +429,10 @@ void cmm_mass(ComMod& com_mod, const double w, const Vector<double>& N, const Ar
   #endif
 
   Vector<double> f(3);
-  double rho = eq.dmn[cDmn].prop.at(PhysicalProperyType::solid_density);
-  f(0) = eq.dmn[cDmn].prop.at(PhysicalProperyType::f_x);
-  f(1) = eq.dmn[cDmn].prop.at(PhysicalProperyType::f_y);
-  f(2) = eq.dmn[cDmn].prop.at(PhysicalProperyType::f_z);
+  double rho = eq.dmn[cDmn].prop.at(PhysicalPropertyType::solid_density);
+  f(0) = eq.dmn[cDmn].prop.at(PhysicalPropertyType::f_x);
+  f(1) = eq.dmn[cDmn].prop.at(PhysicalPropertyType::f_y);
+  f(2) = eq.dmn[cDmn].prop.at(PhysicalPropertyType::f_z);
   #ifdef debug_cmm_mass
   dmsg << "rho: " << rho ;
   dmsg << "f: " << f ;
@@ -444,7 +444,7 @@ void cmm_mass(ComMod& com_mod, const double w, const Vector<double>& N, const Ar
   if (com_mod.cmmVarWall) { 
     ht = vwp(0);
   } else { 
-    ht = eq.dmn[cDmn].prop.at(PhysicalProperyType::shell_thickness);
+    ht = eq.dmn[cDmn].prop.at(PhysicalPropertyType::shell_thickness);
   }
 
   double wl = w * ht * rho;
@@ -497,7 +497,7 @@ void cmm_stiffness(ComMod& com_mod, const Array<double>& Nxi, const Array<double
   const double dt = com_mod.dt;
   const auto cDmn = com_mod.cDmn;
 
-  double nu = eq.dmn[cDmn].prop.at(PhysicalProperyType::poisson_ratio);
+  double nu = eq.dmn[cDmn].prop.at(PhysicalPropertyType::poisson_ratio);
   double ht, elM; 
 
   if (com_mod.cmmVarWall) {
@@ -507,8 +507,8 @@ void cmm_stiffness(ComMod& com_mod, const Array<double>& Nxi, const Array<double
     // elasticity modulus
     elM = vwp[1];
   } else { 
-    ht = eq.dmn[cDmn].prop.at(PhysicalProperyType::shell_thickness);
-    elM = eq.dmn[cDmn].prop.at(PhysicalProperyType::elasticity_modulus);
+    ht = eq.dmn[cDmn].prop.at(PhysicalPropertyType::shell_thickness);
+    elM = eq.dmn[cDmn].prop.at(PhysicalPropertyType::elasticity_modulus);
   }
 
   double lam = elM /(1.0 - nu*nu);

Code/Source/solver/consts.cpp

@@ -219,6 +219,8 @@ const std::map<std::string,EquationType> equation_name_to_type = {
 
     {"shell", EquationType::phys_shell},
 
+    {"darcy", EquationType::phys_darcy},
+
     {"solid_heat", EquationType::phys_heatS},
     {"heatS", EquationType::phys_heatS},
     {"laplace", EquationType::phys_heatS},

Code/Source/solver/consts.h

@@ -287,9 +287,11 @@ enum class EquationType
   phys_CMM = 209, 
   phys_CEP = 210,
   phys_ustruct = 211,  // Nonlinear elastodynamics using mixed VMS-stabilized formulation 
-  phys_stokes = 212
+  phys_stokes = 212,
+  phys_darcy = 213
 };
 
+constexpr auto Equation_darcy = EquationType::phys_darcy;
 constexpr auto Equation_CMM = EquationType::phys_CMM;
 constexpr auto Equation_CEP = EquationType::phys_CEP;
 constexpr auto Equation_fluid = EquationType::phys_fluid;
@@ -347,6 +349,7 @@ enum class OutputNameType {
   outGrp_activeTensionFibers = 529,
   outGrp_activeTensionSheets = 530,
   outGrp_activeTensionNormal = 531,
+  outGrp_darcyFlux = 532,
 
   out_velocity = 599,
   out_pressure = 598,
@@ -380,7 +383,9 @@ enum class OutputNameType {
   out_fibStretchRate = 570,
   out_activeTensionFibers = 569,
   out_activeTensionSheets = 568,
-  out_activeTensionNormal = 567
+  out_activeTensionNormal = 567,
+  out_darcyPressure = 566,
+  out_darcyFlux = 565
 };
 
 /// @brief Simulation output file types. 
@@ -396,7 +401,7 @@ extern const std::map<std::string,OutputType> output_type_name_to_type;
 
 /// @brief Possible physical properties. Current maxNPror is 20.
 //
-enum class PhysicalProperyType 
+enum class PhysicalPropertyType 
 {
   NA = 0, 
   fluid_density = 1, 
@@ -413,7 +418,10 @@ enum class PhysicalProperyType
   shell_thickness = 12, 
   ctau_M = 13,                 // stabilization coeffs. for USTRUCT (momentum, continuity)
   ctau_C = 14,
-  inverse_darcy_permeability = 15
+  brinkman_inverse_permeability = 15,
+  darcy_permeability = 16,
+  darcy_compressibility = 17,
+  darcy_fluid_viscosity = 18
 };
 
 enum class PreconditionerType