Implement Darcy Physics

[mmegally]

Current situation

See issue #616 on introducing darcy physics.

Release Notes

• Adding new equation type "darcy" to solve for pressure/flows in porous domains
• Test cases added (referenced in issue)

Code of Conduct & Contributing Guidelines

• I agree to follow the Code of Conduct and Contributing Guidelines.

[ktbolt]

@mmegally You've changed 47 files for this PR. It would be helpful to describe at a high level in #616 the implementation, how this is used (XML) and any related issues for data files, linear solver, etc.

[zasexton]

@ktbolt i think the majority of these file changes have to do with adding any new physics into the current procedural framework. Unfortunate but expected at this point in my opinion. @mmegally I'll do my review in a few hours when I have some time. From my first glance id say that you should remove the sbatch slurm scripts from the commits. These shouldn't be in the testing infrastructure. I'll give a more formal review in a bit though

[mmegally]

@ktbolt I went ahead and added a new comment to issue #616 with some details on the file change split, implementation, and XML usage. Let me know if there is anything else I can clarify. To reiterate my comment from there, 31/47 of the file changes are related to the test cases. I will also go ahead and remove the SLURMs from the test folder as zack mentioned right now

[ktbolt]

@mmegally Got it, thanks for adding the Issue comments !

[ktbolt]

@mmegally The darcy_convergence tests are not part of the CI tests so shouldn't they be stored somewhere else ?

Also rename the svFSI.xml files to solver.xml.

[codecov]

Codecov Report

❌ Patch coverage is 50.97561% with 201 lines in your changes missing coverage. Please review.
✅ Project coverage is 72.56%. Comparing base (11e2ee1) to head (5682189).
⚠️ Report is 4 commits behind head on main.

Files with missing lines	Patch %	Lines
Code/Source/solver/darcy.cpp	0.00%	120 Missing ⚠️
Code/Source/solver/set_equation_props.h	86.84%	15 Missing ⚠️
Code/Source/solver/post.cpp	26.31%	14 Missing ⚠️
Code/Source/solver/read_files.cpp	63.33%	11 Missing ⚠️
Code/Source/solver/stokes.cpp	41.66%	7 Missing ⚠️
Code/Source/solver/eq_assem.cpp	0.00%	5 Missing ⚠️
Code/Source/solver/l_elas.cpp	54.54%	5 Missing ⚠️
Code/Source/solver/shells.cpp	66.66%	5 Missing ⚠️
Code/Source/solver/ustruct.cpp	54.54%	5 Missing ⚠️
Code/Source/solver/heats.cpp	50.00%	3 Missing ⚠️
... and 5 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #620      +/-   ##
==========================================
- Coverage   72.86%   72.56%   -0.30%     
==========================================
  Files         258      259       +1     
  Lines       39498    39658     +160     
  Branches     6731     6745      +14     
==========================================
- Hits        28780    28779       -1     
- Misses      10475    10636     +161     
  Partials      243      243              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[zasexton]

Suggested changes for Darcy equation initialization and fiber-flux reconstruction.

[zasexton]

A Darcy-only problem currently sets dFlag=true, causing the predictor to maintain a meaningless Dn state and changing the restart record to include displacement data. The suggested changes below address this by adding the conventional Equation_darcy alias and including Darcy in the first-order/non-displacement classification. With these changes applied, Darcy is handled consistently with the other first-order scalar equations.

Suggested code changes

Add the conventional equation alias in Code/Source/solver/consts.h at line 294:

 constexpr auto Equation_CMM = EquationType::phys_CMM;
 constexpr auto Equation_CEP = EquationType::phys_CEP;
+constexpr auto Equation_darcy = EquationType::phys_darcy;
 constexpr auto Equation_fluid = EquationType::phys_fluid;

Update the non-displacement classification in Code/Source/solver/initialize.cpp at line 443:

-    if (std::set<EquationType>{Equation_fluid, Equation_heatF, Equation_heatS, Equation_CEP, Equation_stokes}.count(eq.phys) == 0) {
+    if (std::set<EquationType>{
+            Equation_CEP,
+            Equation_darcy,
+            Equation_fluid,
+            Equation_heatF,
+            Equation_heatS,
+            Equation_stokes
+        }.count(eq.phys) == 0) {
       dFlag = true;
     }

Result of the suggested change: after equation initialization, a Darcy-only configuration leaves com_mod.dFlag == false.

[zasexton]

I missed this when originally writing the darcy equations since I was unsure of the dflag purpose.

[zasexton]

The current reconstruction reports MBF_flux along global x for a fiber mesh because nn::gnn only populates Nx(0,a) when insd=1. The suggested changes below address this by reconstructing the physical gradient as (dP/ds) * tangent and then applying the Darcy mobility. This was originally something that I overlooked when initially writing these governing equations.

Suggested code changes

Assuming one-dimensional Darcy fiber support is intentional, compute the unit physical tangent in Code/Source/solver/post.cpp beginning at line 917:

     Array<double> ksix(nsd,nsd);
     double Jac = 0.0;
+    Vector<double> fiber_tangent(nsd);
 
     for (int g = 0; g < lM.nG; g++) {
       if (g == 0 || !lM.lShpF) {
         auto Nx_g = lM.Nx.slice(g);
         nn::gnn(eNoN, nsd, insd, Nx_g, xl, Nx, Jac, ksix);
+
+        if (lM.lFib) {
+          fiber_tangent = 0.0;
+          for (int a = 0; a < eNoN; a++) {
+            for (int j = 0; j < nsd; j++) {
+              fiber_tangent(j) =
+                  fiber_tangent(j) + xl(j,a) * Nx_g(0,a);
+            }
+          }
+
+          const double tangent_norm = utils::norm(fiber_tangent);
+          if (utils::is_zero(tangent_norm)) {
+            throw std::runtime_error(
+                "[post] Cannot compute Darcy flux for a degenerate fiber element.");
+          }
+          fiber_tangent = fiber_tangent / tangent_norm;
+        }
       }

Replace the Darcy flux block beginning at Code/Source/solver/post.cpp:1016:

       } else if (outGrp == OutputNameType::outGrp_mbfFlx) {
-        double kappa = eq.dmn[cDmn].prop[PhysicalProperyType::permeability];
-        int i = eq.s;
-        Vector<double> q(nsd);
-        for (int a = 0; a < eNoN; a++) {
-          for (int j = 0; j < nsd; j++) {
-            q(j) = q(j) + Nx(j, a) * yl(i, a);
+        const double permeability =
+            eq.dmn[cDmn].prop[PhysicalProperyType::permeability];
+        const double viscosity =
+            eq.dmn[cDmn].prop[PhysicalProperyType::darcy_fluid_viscosity];
+        const double mobility = permeability / viscosity;
+        const int i = eq.s;
+
+        Vector<double> grad_p(nsd);
+
+        if (lM.lFib) {
+          double dp_ds = 0.0;
+          for (int a = 0; a < eNoN; a++) {
+            dp_ds = dp_ds + Nx(0,a) * yl(i,a);
+          }
+          for (int j = 0; j < nsd; j++) {
+            grad_p(j) = dp_ds * fiber_tangent(j);
+          }
+        } else {
+          for (int a = 0; a < eNoN; a++) {
+            for (int j = 0; j < nsd; j++) {
+              grad_p(j) = grad_p(j) + Nx(j,a) * yl(i,a);
+            }
           }
         }
+
         for (int j = 0; j < nsd; j++) {
-          lRes(j) = -kappa * q(j);
+          lRes(j) = -mobility * grad_p(j);
         }

The replacement intentionally includes the viscosity correction so the reconstructed flux implements -(K/mu) grad(P) for both fiber and volume meshes.

[zasexton]

Source_term is currently folded into Pd and the complete expression is multiplied by Media_compressibility. A nonzero source therefore has no effect when the default compressibility is zero, and otherwise behaves as an offset to the pressure time derivative. The changes below separate the volumetric source from the storage term in all three kernels. With these changes applied, the residual represents rho * beta * dp/dt - rho * source plus the Darcy diffusion term, so Source_term remains active in the default steady/incompressible configuration.

Suggested code changes

Update the 1-D kernel beginning at Code/Source/solver/darcy.cpp:164:

-  double Pd = -source;
+  double p_dot = 0.0;
   double Px = 0.0;

   for (int a = 0; a < eNoN; a++) {
-    Pd = Pd + N(a) * al(i, a);
+    p_dot = p_dot + N(a) * al(i, a);
     Px = Px + Nx(0, a) * yl(i, a);
   }

   for (int a = 0; a < eNoN; a++) {
-    lR(0, a) = lR(0, a) + w * (rho_0 * beta_0 * N(a) * Pd +
-                               (((k * rho_0) / mu) * (Nx(0, a) * Px)));
+    lR(0, a) = lR(0, a) +
+        w * (rho_0 * N(a) * (beta_0 * p_dot - source) +
+             ((k * rho_0) / mu) * Nx(0, a) * Px);
     for (int b = 0; b < eNoN; b++) {
       lK(0, a, b) = lK(0, a, b) + wl * (rho_0 * beta_0 * N(a) * N(b) * amd +
                                        ((((rho_0 * k) / mu) * (Nx(0, a) * Nx(0, b)))));

Update the 2-D kernel beginning at Code/Source/solver/darcy.cpp:219:

-  double Pd = -source;
+  double p_dot = 0.0;
   Vector<double> Px(nsd);

   for (int a = 0; a < eNoN; a++) {
-    Pd = Pd + N(a)*al(i,a);
+    p_dot = p_dot + N(a)*al(i,a);
     Px(0) = Px(0) + Nx(0,a)*yl(i,a);
     Px(1) = Px(1) + Nx(1,a)*yl(i,a);
   }

   for (int a = 0; a < eNoN; a++) {
-    lR(0,a) = lR(0,a) + w*(rho_0*beta_0*N(a)*Pd + (((k*rho_0)/mu)*(Nx(0,a)*Px(0)
-                                                                   + Nx(1,a)*Px(1))));
+    lR(0,a) = lR(0,a) +
+        w * (rho_0 * N(a) * (beta_0 * p_dot - source) +
+             ((k * rho_0) / mu) *
+                 (Nx(0,a) * Px(0) + Nx(1,a) * Px(1)));
     for (int b = 0; b < eNoN; b++) {
       lK(0,a,b) = lK(0,a,b) + wl*(rho_0*beta_0*N(a)*N(b)*amd +
                                   ((((rho_0*k)/mu)*(Nx(0,a)*Nx(0,b) +

Update the 3-D kernel beginning at Code/Source/solver/darcy.cpp:276:

-  double Pd = -source;
+  double p_dot = 0.0;
   Vector<double> Px(nsd);

   for (int a = 0; a < eNoN; a++) {
-    Pd = Pd + N(a) * al(i,a);
+    p_dot = p_dot + N(a) * al(i,a);
     Px(0) = Px(0) + Nx(0,a) * yl(i,a);
     Px(1) = Px(1) + Nx(1,a) * yl(i,a);
     Px(2) = Px(2) + Nx(2,a) * yl(i,a);
   }

   for (int a = 0; a < eNoN; a++) {
-    lR(0,a) = lR(0, a) + w*(rho_0*beta_0*N(a)*Pd +
-                            ((k*rho_0)/mu)*(Nx(0,a)*Px(0) + Nx(1,a)*Px(1) + Nx(2,a)*Px(2)));
+    lR(0,a) = lR(0, a) +
+        w * (rho_0 * N(a) * (beta_0 * p_dot - source) +
+             ((k * rho_0) / mu) *
+                 (Nx(0,a) * Px(0) + Nx(1,a) * Px(1) +
+                  Nx(2,a) * Px(2)));
     for (int b = 0; b < eNoN; b++) {
       lK(0,a,b) = lK(0,a,b) + wl*(rho_0*beta_0*N(a)*N(b)*amd +

For Media_compressibility == 0, the source contribution at test-function node a becomes -w * rho_0 * N(a) * source instead of zero. The source remains independent of compressibility, while the existing tangent matrix remains correct because Source_term is prescribed data rather than an unknown.

[zasexton]

Darcy registers PhysicalProperyType::fluid_density, but this branch currently reads the generic <Density> value for every equation except CMM. The added Darcy inputs set <Fluid_density>1</Fluid_density>, so their requested value is silently replaced by the default <Density> value of 0.5. The change below routes Darcy through the explicit fluid_density parameter alongside CMM. With it applied, the documented Darcy input populates dmn.prop[PhysicalProperyType::fluid_density] as intended.

Suggested code change

Update Code/Source/solver/read_files.cpp:1552:

           case PhysicalProperyType::fluid_density:
-            if (lEq.phys == EquationType::phys_CMM) {
+            if (lEq.phys == EquationType::phys_CMM ||
+                lEq.phys == EquationType::phys_darcy) {
               rtmp = domain_params->fluid_density.value();
             } else {
               rtmp = domain_params->density.value();
             }
           break;

For the four added Darcy XML inputs, <Fluid_density>1</Fluid_density> produces a stored Darcy fluid density of 1.0 rather than the unrelated default <Density> value of 0.5.

[zasexton]

Darcy currently registers and parses solid_density, porosity, porosity_pressure, fluid_compressibility, and density_pressure, although no Darcy kernel consumes them. The cleanup below removes the unused properties from the Darcy property list and removes the four newly introduced parser controls that have no implementation. These other parameters were holdovers from when I was trying different models; they can be removed.

Suggested code changes

Restrict the Darcy property list in Code/Source/solver/set_equation_props.h:287 to properties consumed by assembly or output:

     propL[0][0] = PhysicalProperyType::permeability;
     propL[1][0] = PhysicalProperyType::source_term;
-    propL[2][0] = PhysicalProperyType::solid_density;
-    propL[3][0] = PhysicalProperyType::porosity;
-    propL[4][0] = PhysicalProperyType::fluid_density;
-    propL[5][0] = PhysicalProperyType::porosity_pressure;
-    propL[6][0] = PhysicalProperyType::media_compressibility;
-    propL[7][0] = PhysicalProperyType::fluid_compressibility;
-    propL[8][0] = PhysicalProperyType::darcy_fluid_viscosity;
-    propL[9][0] = PhysicalProperyType::density_pressure;
+    propL[2][0] = PhysicalProperyType::fluid_density;
+    propL[3][0] = PhysicalProperyType::media_compressibility;
+    propL[4][0] = PhysicalProperyType::darcy_fluid_viscosity;

Remove the newly added, unused XML parameter registrations from Code/Source/solver/Parameters.cpp:2049:

   set_parameter("Permeability", 0.0, !required, permeability);
-  set_parameter("Porosity", 0.0, !required, porosity);
-  set_parameter("Porosity_pressure", 0.0, !required, porosity_pressure);
   set_parameter("Media_compressibility", 0.0, !required, media_compressibility);
-  set_parameter("Fluid_compressibility", 0.0, !required, fluid_compressibility);
   set_parameter("Darcy_fluid_viscosity", 1.0, !required, darcy_fluid_viscosity);
-  set_parameter("Density_pressure", 0.0, !required, density_pressure);

Remove their storage members from Code/Source/solver/Parameters.h:1627:

     Parameter<double> permeability;
-    Parameter<double> porosity;
-    Parameter<double> porosity_pressure;
     Parameter<double> media_compressibility;
-    Parameter<double> fluid_compressibility;
     Parameter<double> darcy_fluid_viscosity;
-    Parameter<double> density_pressure;

Remove the unused physical-property identifiers and retain explicit values for the supported identifiers in Code/Source/solver/consts.h:421:

   inverse_darcy_permeability = 15,
   permeability = 16,
-  porosity = 17,
-  porosity_pressure = 18,
   media_compressibility = 19,
-  fluid_compressibility = 20,
-  darcy_fluid_viscosity = 21,
-  density_pressure = 22
+  darcy_fluid_viscosity = 21
 };

Remove the corresponding no-op cases from Code/Source/solver/read_files.cpp:1588:

           case PhysicalProperyType::permeability:
             rtmp = domain_params->permeability.value();
           break;

-          case PhysicalProperyType::porosity:
-            rtmp = domain_params->porosity.value();
-          break;
-
-          case PhysicalProperyType::porosity_pressure:
-            rtmp = domain_params->porosity_pressure.value();
-          break;
-
           case PhysicalProperyType::media_compressibility:
             rtmp = domain_params->media_compressibility.value();
           break;

-          case PhysicalProperyType::fluid_compressibility:
-            rtmp = domain_params->fluid_compressibility.value();
-          break;
-
           case PhysicalProperyType::darcy_fluid_viscosity:
             rtmp = domain_params->darcy_fluid_viscosity.value();
           break;
-
-          case PhysicalProperyType::density_pressure:
-            rtmp = domain_params->density_pressure.value();
-          break;

The implemented Darcy model retains exactly five material inputs, and each has a runtime consumer:

• Permeability: Darcy assembly and flux output
• Source_term: Darcy residual
• Fluid_density: Darcy residual and tangent
• Media_compressibility: transient storage term
• Darcy_fluid_viscosity: Darcy mobility in assembly and flux output

The removed names are no longer presented as supported Darcy controls.

[michelebucelli]

Thank you @mmegally! I left a few comments.

[michelebucelli]

This change seems to be unrelated to the PR. I suggest reverting to keep all the active stress files close together in the list.

[michelebucelli]

Might be a good opportunity to fix this typo and rename this to PhysicalPropertyType.

[michelebucelli]

Since there's a correspondence between entries of the outGrp category and of the out category, I recommend using the same naming convention for both (that is, either both mbfFlx or mbfFlux). In a similar way, I'd suggest to use consistent capitalization of MBF (i.e. if it is all caps in out_MBF, I think it should be all caps in out_MBFFlx and outGrp_MBFFlx).

Finally, I would suggest avoiding abbreviations when possible, unless they're very established (that is to say: I don't mind mbf, but I think Flx should be Flux).

(These are not just meant as cosmetic changes, but also to prevent inexperienced users from thinking that flx and flux are two different things since they have two different names)

[michelebucelli]

Unrelated to this PR, but it's been bugging me for a while 😅 (also discussed in issue #548)

@ktbolt Do you know whether this enum class could be split into two enum classes, one for the outGrp_* entries and one for the out_* entries? It seems to me that they are not interchangeable, based on the way they are used, so it would seem they should be separate.

[ktbolt]

@michelebucelli The OutputNameType enum class represents constants for output but the outGrp_ and out_ constants really do represent two different things. So go ahead and split OutputNameType up if you would like to do that now; all of the output code will definitely be refactored.

[michelebucelli]

What is the difference between MBF and mbfFlux? Expanding the acronym, the latter reads "myocardial blood flow flux", and I'm not sure about its physical meaning. Perhaps with a bit of clarification we can come up with a more self-documenting name (or, failing that, we can document it).

[michelebucelli]

Is there any specific reason why this formulation is divided by the permeability $K$, as opposed to keeping it in the diffusion term, which I think is the more common approach?

The advantage of the latter would be that it lends more naturally to extensions (e.g. heterogeneous or anisotropic diffusion), which may not be relevant now but might become so in the future. The current formulation instead only works if $K$ is a constant scalar.

[michelebucelli]

I would suggest also defining $p_{source}$ and $p_{sink}$ in this list, and the two coefficients $\beta_0$ and $\beta_1$. I would also suggest clarifying for each what is a problem unknown, what is a quantity derived from the unknown (e.g. $u = -K\nabla p$) and what is a user-specified parameter. For user-specified parameter, it is also a good idea to document the set they belong to (e.g. $K \in (0, +\infty)$ and so on), so that the reader can understand what is spatially varying and what is assumed to be constant.

Some of this is documented by the assumptions above, e.g. "isotropic permeability", but I think it's not a bad idea to include these assumptions in the formal statement of the equation as well.

[michelebucelli]

I know this is present also elsewhere, but: should we consider cleaning this up? As far as I understand it, NURBS are currently unsupported, I'm not sure there's plans to support them in the immediate future.

If we want to keep it as a placeholder, I'd at least suggest to throw a not-implemented exception here, instead of silently doing nothing, so that if the code somehow gets here the user/developer is made very aware of it.

[michelebucelli]

If I am reading this code right, the problem assembled here does not match the one stated in the documentation at the top of this file. The differences I can notice are the following.

1. The terms assembled here do not have anything that corresponds to the sink term $\beta_1 (p - p_{sink})$.
2. The term multiplied by $\beta_0$ is not $\beta_0 (p - p_{source})$, but rather $\beta_0 \rho_0 \frac{\partial p}{\partial t}$
3. There's a couple of coefficients ($\rho_0$ and $\mu$) that rescale the permeability $K$.

It seems to me that the strong form of the problem assembled here is

$$ \rho_0\beta_0\frac{\partial p}{\partial t} - \nabla \cdot\left(\frac{K\rho_0} {\mu}\nabla p\right) = \rho_0 p_{source} $$

Is this right?

If so, I would suggest confirming that this is indeed the problem that is meant to be solved, and update the documentation accordingly.

Additionally, in the above equation $\rho_0$ can be simplified out (assuming it is homogeneous and constant, which seems to be the case), so I suggest to consider whether it makes sense to remove it altogether. Similarly, $\mu$ only acts as a rescaling term on $K$, so just asking the user to provide a rescaled permemability might make for a simpler parameterization.

The same comments apply to darcy_2d and darcy_3d below.

[michelebucelli]

It seems to me that darcy_1d, darcy_2d and darcy_3d are essentially identical functions, but for the way they unroll the loops needed to compute the scalar products (and maybe the way they assemble the pressure gradient Px).

I think it would be rather easy to merge this into a single function, in one of these ways.

1. Templating the darcy function over the spatial dimension, and then using if constexpr statements to switch between the optimized cases for the individual dimensions. For example, the interpolation of the pressure gradient could look something like

   for (int a = 0; a < eNoN; a++) {
     p_dot = p_dot + N(a)*al(i,a);
   
     Px(0) = Px(0) + Nx(0,a)*yl(i,a);     
     if constexpr (nsd > 1)
       Px(1) = Px(1) + Nx(1,a)*yl(i,a);
     if constexpr (nsd > 2)
       Px(2) = Px(2) + Nx(2,a)*yl(i,a);
   }

2. Better than the above, move this dimension-dependent operations (interpolation of the gradient, dot products) to dedicated helper functions, and template those on spatial dimension if necessary (I'm pretty sure that the dot product between vectors is already implemented somewhere). If manual loop unrolling is needed for efficiency, this can be hidden away in the helper functions, which will be good for readability.

Either way, this would ensure more robustly that the 1D, 2D and 3D variants of the equation behave consistently, and it will make it easier to modify the formulation of the equation, should this be needed in the future.

I am aware that the same pattern is present in several other equations, so this is perhaps a more general comment. But it might be a good idea to start applying this to new equations, leaving the improvement of other equations for incoming refactoring.