Implement Darcy Physics by mmegally · Pull Request #620 · SimVascular/svMultiPhysics

mmegally · 2026-08-26T19:47:19Z

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 · 2026-08-26T20:10:10Z

@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.

… after fork)

zasexton · 2026-08-26T22:32:41Z

@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 · 2026-08-26T22:36:36Z

@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 · 2026-08-26T22:43:09Z

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

ktbolt · 2026-08-26T23:22:27Z

@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 · 2026-08-27T00:59:20Z

Codecov Report

❌ Patch coverage is 52.01699% with 226 lines in your changes missing coverage. Please review.
✅ Project coverage is 73.11%. Comparing base (00f5c92) to head (c159793).

Files with missing lines	Patch %	Lines
Code/Source/solver/darcy.cpp	0.00%	142 Missing ⚠️
Code/Source/solver/set_equation_props.h	85.34%	17 Missing ⚠️
Code/Source/solver/post.cpp	26.31%	14 Missing ⚠️
Code/Source/solver/read_files.cpp	67.74%	10 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 6 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #620      +/-   ##
==========================================
- Coverage   73.45%   73.11%   -0.35%     
==========================================
  Files         267      268       +1     
  Lines       40041    40241     +200     
  Branches     6706     6729      +23     
==========================================
+ Hits        29413    29423      +10     
- Misses      10385    10575     +190     
  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.

…k's suggested changes)

michelebucelli

Thank you @mmegally! I left a few comments.

michelebucelli · 2026-08-29T16:37:36Z

+void darcy_1d(ComMod& com_mod, const int eNoN, const double w, const Vector<double>& N, const Array<double>& Nx,
+              const Array<double>& al, const Array<double>& yl, Array<double>& lR, Array3<double>& lK)
+{
+  using namespace consts;
+
+  const int cEq = com_mod.cEq;
+  auto& eq = com_mod.eq[cEq];
+  const int cDmn = com_mod.cDmn;
+  auto& dmn = eq.dmn[cDmn];
+  const double dt = com_mod.dt;
+  const int i = eq.s;
+
+  double k = dmn.prop.at(PhysicalProperyType::permeability);
+  double source = dmn.prop.at(PhysicalProperyType::source_term);
+  double beta_0 = dmn.prop.at(PhysicalProperyType::media_compressibility);
+  double rho_0 = dmn.prop.at(PhysicalProperyType::fluid_density);
+  double mu = dmn.prop.at(PhysicalProperyType::darcy_fluid_viscosity);
+
+  double T1 = eq.af * eq.gam * dt;
+  double amd = eq.am / T1;
+  double wl = w * T1;
+
+  double p_dot = 0.0;
+  double Px = 0.0;
+
+  for (int a = 0; a < eNoN; 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 * 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)))));
+    }
+  }
+}


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.

The terms assembled here do not have anything that corresponds to the sink term $\beta_1 (p - p_{sink})$.

The term multiplied by $\beta_0$ is not $\beta_0 (p - p_{source})$, but rather $\beta_0 \rho_0 \frac{\partial p}{\partial t}$

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.

Opting to keep the viscosity and permeability separate at this time since there will likely be some ongoing discussions with the wet-lab experimentalists that will likely want some 1:1 correspondence for their measurements to input parameters. Also keeping the 2d and 3d assembly kernels separate at this time to mirror the other physics files; although I agree that future refactoring should cleanup these to have better concision of code.

Also yes the equation that we should be solving for this model is:

$$\rho\beta\frac{\partial p}{\partial t}-\nabla\cdot\left(\frac{\rho K}{\mu}\nabla p\right)=\rho s$$.

$\beta$ should be the storage/compressibility and $s$ is the prescribed volumetric source. You are correct that we do not implement the $p_{source}$ and $p_{sink}$ exchange terms from the original documentation that was written. That documentation has been updated to reflect the current assembly implementation.

Thanks for clarifying! I assume that you opted to keep the density separate for the same reason as viscosity and compressibility, right?

…sicalProperyType' in other physics modules)

zasexton · 2026-09-02T20:14:05Z

About to push some commits to this fork then we can elevate this draft to a full PR. in the meantime @mmegally can you resolve the active comments that you have addressed along with their appropriate commit tags? you've been doing a great job at fixing up the code so far :)

…nges

Resolve remaining Darcy review comments

zasexton · 2026-09-03T20:35:30Z

@mmegally let's mark this as "Ready to review"

claude

Claude Code Review

This pull request is from a fork — automated review is disabled. A repository maintainer can comment @claude review to run a one-time review.

zasexton

In my opinion we only piece requiring a bit more review is the zero initialization for the 2D cases to make sure that we are not getting heap error allocations and out-of-bounds. I'll fork this branch and double-check now to make any appropriate changes, if necessary.

michelebucelli

Thanks! I did a second pass, I only have a few minor suggestions left.

michelebucelli · 2026-09-10T22:01:08Z

+void darcy_1d(ComMod& com_mod, const int eNoN, const double w, const Vector<double>& N, const Array<double>& Nx,
+              const Array<double>& al, const Array<double>& yl, Array<double>& lR, Array3<double>& lK)
+{
+  using namespace consts;
+
+  const int cEq = com_mod.cEq;
+  auto& eq = com_mod.eq[cEq];
+  const int cDmn = com_mod.cDmn;
+  auto& dmn = eq.dmn[cDmn];
+  const double dt = com_mod.dt;
+  const int i = eq.s;
+
+  double k = dmn.prop.at(PhysicalProperyType::permeability);
+  double source = dmn.prop.at(PhysicalProperyType::source_term);
+  double beta_0 = dmn.prop.at(PhysicalProperyType::media_compressibility);
+  double rho_0 = dmn.prop.at(PhysicalProperyType::fluid_density);
+  double mu = dmn.prop.at(PhysicalProperyType::darcy_fluid_viscosity);
+
+  double T1 = eq.af * eq.gam * dt;
+  double amd = eq.am / T1;
+  double wl = w * T1;
+
+  double p_dot = 0.0;
+  double Px = 0.0;
+
+  for (int a = 0; a < eNoN; 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 * 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)))));
+    }
+  }
+}


Thanks for clarifying! I assume that you opted to keep the density separate for the same reason as viscosity and compressibility, right?

…laces and adding parameter dimensions to documentation

Fix Darcy initial pressure indexing

…g exceptions from svmp namespace

mmegally added 3 commits August 21, 2026 13:35

Add darcy solver implementation and test case

13a64a2

Fix function signatures in header file

49c5bf6

add convergence tests for darcy implementation

e284118

mmegally requested review from michelebucelli and zasexton August 26, 2026 19:47

Add Regazzoni model to validate Ubuntu/macOS integration tests (added…

b1cd373

… after fork)

Remove unintended SLURM files

bae4ba1

Merge branch 'main' into darcy-solver

12e5427

zasexton reviewed Aug 27, 2026

View reviewed changes

Comment thread Code/Source/solver/consts.h

Comment thread Code/Source/solver/post.cpp Outdated

zasexton reviewed Aug 27, 2026

View reviewed changes

Comment thread Code/Source/solver/darcy.cpp Outdated

zasexton reviewed Aug 27, 2026

View reviewed changes

Comment thread Code/Source/solver/set_equation_props.h Outdated

zasexton reviewed Aug 27, 2026

View reviewed changes

Comment thread Code/Source/solver/set_equation_props.h Outdated

Addressing pull request comments (removed tests from CI suite and Zac…

e65e247

…k's suggested changes)

mmegally requested a review from zasexton August 27, 2026 22:26

michelebucelli reviewed Aug 29, 2026

View reviewed changes

mmegally and others added 4 commits September 1, 2026 10:51

Addressing first set of comments from Michele (includes changing 'Phy…

993bbc5

…sicalProperyType' in other physics modules)

fixing unit test PhysicalPropertyType dependency

d063397

Reducing darcy coverage for 2D/3D and changing output names

5682189

Merge branch 'main' into darcy-solver

dc2b9bb

mmegally and others added 4 commits September 3, 2026 11:48

Merge branch 'main' into darcy-solver

294af52

Removing empty NURBS block and fixing recent PhysicalPropertyType cha…

68fe6ae

…nges

docs: clarify Darcy pressure and flux outputs

59f6444

refactor: scope Darcy material parameter names

6189ba6

zasexton added 8 commits September 3, 2026 12:41

fix: validate Darcy material parameters

1d9d562

docs: align the Darcy equation with assembly

99c80d1

docs: state the Darcy single-field formulation

3344986

docs: retain permeability in the diffusion operator

5b094b4

docs: define Darcy model quantities and ranges

7b4ca1c

docs: document the Darcy namespace and references

a785875

fix: reject embedded-line Darcy use

5b7d582

refactor: clarify Darcy flux reconstruction

8eb5ec5

zasexton mentioned this pull request Sep 3, 2026

Resolve remaining Darcy review comments mmegally/svMultiPhysics#1

Merged

Merge pull request #1 from zasexton/review/pr-620-michele

71d99b2

Resolve remaining Darcy review comments

mmegally marked this pull request as ready for review September 3, 2026 20:36

claude Bot reviewed Sep 3, 2026

View reviewed changes

zasexton added the enhancement New feature or request label Sep 3, 2026

mmegally requested review from ktbolt and michelebucelli September 4, 2026 19:30

mmegally mentioned this pull request Sep 8, 2026

Tissue perfusion model #71

Open

Merge branch 'main' into darcy-solver

c83bcec

zasexton reviewed Sep 9, 2026

View reviewed changes

Comment thread Code/Source/solver/initialize.cpp

zasexton reviewed Sep 9, 2026

View reviewed changes

fix: initialize pressure using equation DOFs

a87e729

zasexton self-requested a review September 10, 2026 18:30

michelebucelli reviewed Sep 10, 2026

View reviewed changes

zasexton and others added 7 commits September 10, 2026 15:09

Merge branch 'main' into darcy-solver

3b014e7

First set of comments: dropping element validation in inappropriate p…

3d3d457

…laces and adding parameter dimensions to documentation

Merge pull request #2 from zasexton/review/pr-620-comment-3887103709

e8bda34

Fix Darcy initial pressure indexing

Renaming Darcy_media_compressibility to Darcy_compressibility & addin…

a9bbd48

…g exceptions from svmp namespace

Fixing instance of darcy_compressibility

a68f7c4

fixing incomplete removal of runtime error call from std namespace

623192f

fixing svmp exception calls

c159793

Uh oh!

Conversation

mmegally commented Aug 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Current situation

Release Notes

Code of Conduct & Contributing Guidelines

Uh oh!

ktbolt commented Aug 26, 2026

Uh oh!

zasexton commented Aug 26, 2026

Uh oh!

mmegally commented Aug 26, 2026

Uh oh!

ktbolt commented Aug 26, 2026

Uh oh!

ktbolt commented Aug 26, 2026

Uh oh!

codecov Bot commented Aug 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

zasexton left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

michelebucelli left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

michelebucelli Aug 29, 2026

Choose a reason for hiding this comment

Uh oh!

zasexton Sep 3, 2026

Choose a reason for hiding this comment

Uh oh!

zasexton Sep 9, 2026

Choose a reason for hiding this comment

Uh oh!

michelebucelli Sep 10, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

zasexton commented Sep 2, 2026

Uh oh!

zasexton commented Sep 3, 2026

Uh oh!

claude Bot left a comment

Choose a reason for hiding this comment

Claude Code Review

Uh oh!

Uh oh!

zasexton left a comment

Choose a reason for hiding this comment

Uh oh!

michelebucelli left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

michelebucelli Sep 10, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

mmegally commented Aug 26, 2026 •

edited

Loading

codecov Bot commented Aug 27, 2026 •

edited

Loading