From 13a64a26c102b835c07cd313299205b22c8c9ac9 Mon Sep 17 00:00:00 2001 From: Michael Date: Fri, 21 Aug 2026 13:35:19 -0700 Subject: [PATCH 01/28] Add darcy solver implementation and test case --- Code/Source/solver/CMakeLists.txt | 4 +- Code/Source/solver/Parameters.cpp | 8 + Code/Source/solver/Parameters.h | 7 + Code/Source/solver/consts.cpp | 2 + Code/Source/solver/consts.h | 17 +- Code/Source/solver/darcy.cpp | 297 ++++++++++++++++++ Code/Source/solver/darcy.h | 28 ++ Code/Source/solver/eq_assem.cpp | 9 + Code/Source/solver/load_msh.h | 31 +- Code/Source/solver/petsc_impl.cpp | 3 + Code/Source/solver/post.cpp | 15 + Code/Source/solver/read_files.cpp | 32 +- Code/Source/solver/read_msh.h | 49 ++- Code/Source/solver/set_equation_dof.h | 3 +- Code/Source/solver/set_equation_props.h | 31 ++ Code/Source/solver/set_output_props.h | 4 +- Code/Source/solver/txt.cpp | 1 + Code/Source/solver/vtk_xml.cpp | 3 +- .../darcy_cylinder/darcy_validation_plot.png | 3 + .../fluid/darcy_cylinder/flow_solving_darcy | 42 +++ .../darcy_cylinder/mesh-surfaces/bottom.vtp | 3 + .../mesh-surfaces/inner_wall.vtp | 3 + .../mesh-surfaces/outer_wall.vtp | 3 + .../darcy_cylinder/mesh-surfaces/top.vtp | 3 + .../cases/fluid/darcy_cylinder/result_020.vtu | 3 + tests/cases/fluid/darcy_cylinder/svFSI.xml | 92 ++++++ .../cases/fluid/darcy_cylinder/tissue_vol.vtu | 3 + 27 files changed, 674 insertions(+), 25 deletions(-) create mode 100644 Code/Source/solver/darcy.cpp create mode 100644 Code/Source/solver/darcy.h create mode 100644 tests/cases/fluid/darcy_cylinder/darcy_validation_plot.png create mode 100644 tests/cases/fluid/darcy_cylinder/flow_solving_darcy create mode 100644 tests/cases/fluid/darcy_cylinder/mesh-surfaces/bottom.vtp create mode 100644 tests/cases/fluid/darcy_cylinder/mesh-surfaces/inner_wall.vtp create mode 100644 tests/cases/fluid/darcy_cylinder/mesh-surfaces/outer_wall.vtp create mode 100644 tests/cases/fluid/darcy_cylinder/mesh-surfaces/top.vtp create mode 100644 tests/cases/fluid/darcy_cylinder/result_020.vtu create mode 100644 tests/cases/fluid/darcy_cylinder/svFSI.xml create mode 100644 tests/cases/fluid/darcy_cylinder/tissue_vol.vtu diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index b37e8b985..ddb614b29 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -178,6 +178,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 @@ -234,8 +235,7 @@ set(CSRCS active_stress_uniform_unsteady.cpp active_stress_ode.cpp active_stress_nash_panfilov.cpp - active_stress_regazzoni.cpp - + SPLIT.c svZeroD_interface/LPNSolverInterface.h svZeroD_interface/LPNSolverInterface.cpp diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 8f9eed3d6..0be19fd1c 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2044,6 +2044,14 @@ DomainParameters::DomainParameters() { set_parameter("Penalty_parameter", 0.0, !required, penalty_parameter); set_parameter("Poisson_ratio", 0.3, !required, poisson_ratio); + + 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); set_parameter("Relative_tolerance", 1e-4, !required, relative_tolerance); set_parameter("Shell_thickness", 0.0, !required, shell_thickness); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index b18bd4661..325b0ce08 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1623,6 +1623,13 @@ class DomainParameters : public ParameterLists Parameter source_term; Parameter time_step_for_integration; + Parameter permeability; + Parameter porosity; + Parameter porosity_pressure; + Parameter media_compressibility; + Parameter fluid_compressibility; + Parameter darcy_fluid_viscosity; + Parameter density_pressure; // Inverse of Darcy permeability. Default value of 0.0 for Navier-Stokes and non-zero for Navier-Stokes-Brinkman Parameter inverse_darcy_permeability; }; diff --git a/Code/Source/solver/consts.cpp b/Code/Source/solver/consts.cpp index a7ec61a9d..886d669a9 100644 --- a/Code/Source/solver/consts.cpp +++ b/Code/Source/solver/consts.cpp @@ -219,6 +219,8 @@ const std::map 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}, diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index 05983c4b0..c8817548e 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -287,7 +287,8 @@ 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_CMM = EquationType::phys_CMM; @@ -347,6 +348,7 @@ enum class OutputNameType { outGrp_activeTensionFibers = 529, outGrp_activeTensionSheets = 530, outGrp_activeTensionNormal = 531, + outGrp_mbfFlx = 532, out_velocity = 599, out_pressure = 598, @@ -380,7 +382,9 @@ enum class OutputNameType { out_fibStretchRate = 570, out_activeTensionFibers = 569, out_activeTensionSheets = 568, - out_activeTensionNormal = 567 + out_activeTensionNormal = 567, + out_MBF = 566, + out_mbfFlux = 565 }; /// @brief Simulation output file types. @@ -413,7 +417,14 @@ enum class PhysicalProperyType shell_thickness = 12, ctau_M = 13, // stabilization coeffs. for USTRUCT (momentum, continuity) ctau_C = 14, - inverse_darcy_permeability = 15 + inverse_darcy_permeability = 15, + permeability = 16, + porosity = 17, + porosity_pressure = 18, + media_compressibility = 19, + fluid_compressibility = 20, + darcy_fluid_viscosity = 21, + density_pressure = 22 }; enum class PreconditionerType diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp new file mode 100644 index 000000000..4abf8c5a3 --- /dev/null +++ b/Code/Source/solver/darcy.cpp @@ -0,0 +1,297 @@ +/* + This code implements the Darcy Equation for 2D and 3D + problems in perfusion of porous media. + ------------------------------------------------------------- + Assumptions: + - Homogeneous Permeability + - Homogeneous Density + - Isotropic Permeability + - Assumptions of Stokes Flow + - Steady-State + ------------------------------------------------------------- + Strong form of the Single-Compartment Darcy equation: + u = -K∇(P) + ∇⋅u = β0(P_source - P) - β1(P - P_sink) + where: + u -> Volume flux vector + K -> Permeability tensor + P -> Pressure + ------------------------------------------------------------- + Weak form of the Single-Compartment Darcy equation: + -∫(∇q∇P)dΩ - λ∫qPdΩ = ∫qFdΩ - ∫q∇P⋅nvdΓ + where: + q -> Test function + λ -> (β0 + β1)/K + F -> -(β0(P_source) + β1(P_sink))/K + n -> Normal vector to the boundary +*/ + +#include "darcy.h" + +#include "all_fun.h" +#include "mat_fun.h" +#include "nn.h" +#include "utils.h" + +namespace darcy { + +void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR) +{ + for (int a = 0; a < eNoN; a++) { + lR(0,a) = lR(0,a) + w * N(a) * h; + } +} + +void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& solutions) +{ + const auto& Ag = solutions.intermediate.get_acceleration(); + const auto& Yg = solutions.intermediate.get_velocity(); + #define n_debug_construct_darcy + #ifdef debug_construct_darcy + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + + using namespace consts; + + const int nsd = com_mod.nsd; + const int tDof = com_mod.tDof; + const int dof = com_mod.dof; + const int cEq = com_mod.cEq; + const auto& eq = com_mod.eq[cEq]; + auto& cDmn = com_mod.cDmn; + + int eNoN = lM.eNoN; + int insd = nsd; + if (lM.lFib) { + insd = 1; + } + #ifdef debug_construct_darcy + dmsg << "cEq: " << cEq; + dmsg << "cDmn: " << cDmn; + dmsg << "insd: " << insd; + #endif + + Vector ptr(eNoN); + Vector N(eNoN); + Array xl(nsd, eNoN), al(tDof, eNoN), yl(tDof, eNoN); + Array Nx(insd, eNoN), lR(dof, eNoN); + Array3 lK(dof * dof, eNoN, eNoN); + Array ksix(nsd, nsd); + + for (int e = 0; e < lM.nEl; e++) { + cDmn = all_fun::domain(com_mod, lM, cEq, e); + auto cPhys = eq.dmn[cDmn].phys; + if (cPhys != EquationType::phys_darcy) { + continue; + } + + // Update shape function for NURBS + if (lM.eType == ElementType::NRB) { + //CALL NRMNNX(lm, e) + } + + // Create local copies + for (int a = 0; a < eNoN; a++) { + int Ac = lM.IEN(a, e); + ptr(a) = Ac; + + for (int i = 0; i < nsd; i++) { + xl(i, a) = com_mod.x(i, Ac); + } + + for (int i = 0; i < tDof; i++) { + al(i, a) = Ag(i, Ac); + yl(i, a) = Yg(i, Ac); + } + } + + // Gauss integration + lR = 0.0; + lK = 0.0; + double Jac{0.0}; + + 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 (utils::is_zero(Jac)) { + throw std::runtime_error( + "[construct_darcy] Jacobian for element " + std::to_string(e) + " is < 0."); + } + } + + double w = lM.w(g) * Jac; + N = lM.N.col(g); + + if (insd == 3) { + darcy_3d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); + } else if (insd == 2) { + darcy_2d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); + } else if (insd == 1) { + darcy_1d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); + } else { + throw std::runtime_error("[construct_darcy] insd must be 1, 2 or 3."); + } + } + + eq.linear_algebra->assemble(com_mod, eNoN, ptr, lK, lR); + } +} + +void darcy_1d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, + const Array& al, const Array& yl, Array& lR, Array3& 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 Pd = -source; + double Px = 0.0; + + for (int a = 0; a < eNoN; a++) { + Pd = Pd + 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))); + 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))))); + } + } +} + +void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, + const Array& al, const Array& yl, Array& lR, Array3& lK) +{ + #define n_debug_darcy_2d + #ifdef debug_darcy_2d + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "w: " << w; + #endif + using namespace consts; + + const int nsd = com_mod.nsd; + 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; + + #ifdef debug_darcy_2d + dmsg << "k: " << k; + dmsg << "source: " << source; + dmsg << "T1: " << T1; + dmsg << "i: " << i; + dmsg << "wl: " << wl; + #endif + + double Pd = -source; + Vector Px(nsd); + + for (int a = 0; a < eNoN; a++) { + Pd = Pd + 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)))); + 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) + + Nx(1,a)*Nx(1,b))))); + } + } +} + +void darcy_3d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, + const Array& al, const Array& yl, Array& lR, Array3& lK) +{ + #define n_debug_darcy_3d + #ifdef debug_darcy_3d + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "w: " << w; + #endif + using namespace consts; + + const int nsd = com_mod.nsd; + 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; + + #ifdef debug_darcy_3d + dmsg << "k: " << k; + dmsg << "source: " << source; + dmsg << "T1: " << T1; + dmsg << "i: " << i; + dmsg << "wl: " << wl; + #endif + + double Pd = -source; + Vector Px(nsd); + + for (int a = 0; a < eNoN; a++) { + Pd = Pd + 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))); + 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 + + ((k*rho_0)/mu)*(Nx(0,a)*Nx(0,b) + Nx(1,a)*Nx(1,b) + + Nx(2,a)*Nx(2,b))); + } + } +} + +} \ No newline at end of file diff --git a/Code/Source/solver/darcy.h b/Code/Source/solver/darcy.h new file mode 100644 index 000000000..b80535cff --- /dev/null +++ b/Code/Source/solver/darcy.h @@ -0,0 +1,28 @@ +// +// This code implements the Darcy Equation for 2D and 3D +// problems in perfusion of porus media. +// + +#ifndef DARCY_H +#define DARCY_H + +#include "ComMod.h" +#include "SolutionStates.h" + +namespace darcy { + + void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR); + + void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& solutions); + + void darcy_1d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, + const Array& al, const Array& yl, Array& lR, Array3& lK); + + void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, + const Array& al, const Array& yl, Array& lR, Array3& lK); + + void darcy_3d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, + const Array& al, const Array& yl, Array& lR, Array3& lK); +} + +#endif //DARCY_H \ No newline at end of file diff --git a/Code/Source/solver/eq_assem.cpp b/Code/Source/solver/eq_assem.cpp index 9cea03b30..abc09ff43 100644 --- a/Code/Source/solver/eq_assem.cpp +++ b/Code/Source/solver/eq_assem.cpp @@ -11,6 +11,7 @@ #include "cep.h" #include "cmm.h" +#include "darcy.h" #include "fluid.h" #include "fsi.h" #include "heatf.h" @@ -111,6 +112,10 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& heatf::b_heatf(com_mod, eNoN, w, N, y, h, nV, lR, lK); break; + case EquationType::phys_darcy: + darcy::b_darcy(com_mod, eNoN, w, N, h, lR); + break; + case EquationType::phys_lElas: l_elas::b_l_elas(com_mod, eNoN, w, N, h, nV, lR); break; @@ -408,6 +413,10 @@ void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const heats::construct_heats(com_mod, lM, solutions); break; + case EquationType::phys_darcy: + darcy::construct_darcy(com_mod, lM, solutions); + break; + case EquationType::phys_lElas: l_elas::construct_l_elas(com_mod, lM, solutions); break; diff --git a/Code/Source/solver/load_msh.h b/Code/Source/solver/load_msh.h index 530252710..294807aa2 100644 --- a/Code/Source/solver/load_msh.h +++ b/Code/Source/solver/load_msh.h @@ -1,5 +1,32 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. -// SPDX-License-Identifier: BSD-3-Clause +/* Copyright (c) Stanford University, The Regents of the University of California, and others. + * + * All Rights Reserved. + * + * See Copyright-SimVascular.txt for additional details. + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject + * to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS + * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER + * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ #ifndef LOAD_MSH_H #define LOAD_MSH_H diff --git a/Code/Source/solver/petsc_impl.cpp b/Code/Source/solver/petsc_impl.cpp index eed226d7b..916a55c2e 100644 --- a/Code/Source/solver/petsc_impl.cpp +++ b/Code/Source/solver/petsc_impl.cpp @@ -145,6 +145,9 @@ void petsc_create_linearsolver(const consts::SolverType lsType, const consts::Pr case EquationType::phys_stokes: psol[cEq].pre = "ss_"; break; + case EquationType::phys_darcy: + psol[cEq].pre = "dr_"; + break; default: PetscPrintf(MPI_COMM_WORLD, "ERROR : " "equation type %d is not defined.\n", phys); diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index fc3990173..9daccf7f1 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -1011,6 +1011,21 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S } } + // MBF Flux calculation + // + } else if (outGrp == OutputNameType::outGrp_mbfFlx) { + double kappa = eq.dmn[cDmn].prop[PhysicalProperyType::permeability]; + int i = eq.s; + Vector 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); + } + } + for (int j = 0; j < nsd; j++) { + lRes(j) = -kappa * q(j); + } + // Strain tensor invariants calculation // } else if (outGrp == OutputNameType::outGrp_stInv) { diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index befd3b9e4..f1686b493 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -206,8 +206,8 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, if (effective_direction.size() != 0) { if (effective_direction.size() != com_mod.nsd) { - auto effective_size = (std::stringstream() << "(" << effective_direction.size() << ")").str(); - auto space_dim = (std::stringstream() << "(" << com_mod.nsd << ")").str(); + auto effective_size = "(" + std::to_string(effective_direction.size()) + ")"; + auto space_dim = "(" + std::to_string(com_mod.nsd) + ")"; svmp::raise("The size of the effective direction " + effective_size + " does not equal the number of space dimensions " + space_dim); } @@ -1580,6 +1580,34 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& case PhysicalProperyType::inverse_darcy_permeability: rtmp = domain_params->inverse_darcy_permeability.value(); break; + + 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; } lEq.dmn[iDmn].prop[prop] = rtmp; diff --git a/Code/Source/solver/read_msh.h b/Code/Source/solver/read_msh.h index c8d336eed..0712da0d3 100644 --- a/Code/Source/solver/read_msh.h +++ b/Code/Source/solver/read_msh.h @@ -1,11 +1,37 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. -// SPDX-License-Identifier: BSD-3-Clause +/* Copyright (c) Stanford University, The Regents of the University of California, and others. + * + * All Rights Reserved. + * + * See Copyright-SimVascular.txt for additional details. + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject + * to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS + * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER + * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ #ifndef READ_MSH_H #define READ_MSH_H #include "ComMod.h" -#include "SolutionStates.h" #include "Simulation.h" #include "Vector.h" @@ -22,11 +48,11 @@ namespace read_msh_ns { Vector gN; }; - void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const SolutionStates& solutions); - void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const SolutionStates& solutions); - void calc_elem_skew(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const SolutionStates& solutions); + void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag); + void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag); + void calc_elem_skew(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag); - void calc_mesh_props(ComMod& com_mod, const CmMod& cm_mod, const int nMesh, std::vector& mesh, const SolutionStates& solutions); + void calc_mesh_props(ComMod& com_mod, const CmMod& cm_mod, const int nMesh, std::vector& mesh); void calc_nbc(mshType& mesh, faceType& face); @@ -46,18 +72,15 @@ namespace read_msh_ns { void load_var_ini(Simulation* simulation, const ComMod& com_mod); void match_faces(const ComMod& com_mod, const faceType& face1, const faceType& face2, const double tol, utils::stackType& lPrj); - void match_nodes(const ComMod& com_mod, const faceType& lFa, const faceType& pFa, - const double ptol, const int nNds, Array& map); void read_fib_nff(Simulation* simulation, mshType& mesh, const std::string& fName, const std::string& kwrd, const int idx); void read_msh(Simulation* simulation); void set_dmn_id_ff(Simulation* simulation, mshType& mesh, const std::string& file_name); - void set_dmn_id_vtk(Simulation* simulation, mshType& lM, const std::string& file_name, const std::string& kwrd); + void set_dmn_id_vtk(Simulation* simulation, mshType& mesh, const std::string& file_name, const std::string& kwrd); void set_projector(Simulation* simulation, utils::stackType& avNds); - void set_ris_projector(Simulation* simulation); - void set_uris_meshes(Simulation* simulation); - + + }; #endif diff --git a/Code/Source/solver/set_equation_dof.h b/Code/Source/solver/set_equation_dof.h index eb01a31e5..096a5b5cf 100644 --- a/Code/Source/solver/set_equation_dof.h +++ b/Code/Source/solver/set_equation_dof.h @@ -21,6 +21,7 @@ std::map equation_dof_map = {EquationType::phys_FSI, std::make_tuple(nsd+1, "FS") }, {EquationType::phys_mesh, std::make_tuple(nsd, "MS") }, {EquationType::phys_CEP, std::make_tuple(1, "EP") }, - {EquationType::phys_stokes, std::make_tuple(nsd+1, "SS") } + {EquationType::phys_stokes, std::make_tuple(nsd+1, "SS") }, + {EquationType::phys_darcy, std::make_tuple(1, "DR") } }; diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index 1cd236b7f..aa492d44b 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -274,6 +274,37 @@ SetEquationPropertiesMapType set_equation_props = { } }, +//---------------------------// +// phys_darcy // +//---------------------------// +{consts::EquationType::phys_darcy, [](Simulation* simulation, EquationParameters* eq_params, eqType& lEq, EquationProps& propL, + EquationOutputs& outPuts, EquationNdop& nDOP) -> void +{ + using namespace consts; + auto& com_mod = simulation->get_com_mod(); + lEq.phys = consts::EquationType::phys_darcy; + + 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; + + read_domain(simulation, eq_params, lEq, propL); + + nDOP = {2,1,1,0}; + outPuts = {OutputNameType::out_MBF, OutputNameType::out_mbfFlux}; + + // Set solver parameters. + read_ls(simulation, eq_params, SolverType::lSolver_CG, lEq); +} }, + + //---------------------------// // phys_FSI // //---------------------------// diff --git a/Code/Source/solver/set_output_props.h b/Code/Source/solver/set_output_props.h index 35cfbd8c4..493d8684f 100644 --- a/Code/Source/solver/set_output_props.h +++ b/Code/Source/solver/set_output_props.h @@ -59,5 +59,7 @@ std::map output_props_map = {OutputNameType::out_voltage, std::make_tuple(OutputNameType::outGrp_Y, 0, 1, "Membrane_potential") }, {OutputNameType::out_vortex, std::make_tuple(OutputNameType::outGrp_vortex, 0, 1, "Vortex") }, {OutputNameType::out_vorticity, std::make_tuple(OutputNameType::outGrp_vort, 0, maxNSD, "Vorticity") }, - {OutputNameType::out_WSS, std::make_tuple(OutputNameType::outGrp_WSS, 0, maxNSD, "WSS") } + {OutputNameType::out_WSS, std::make_tuple(OutputNameType::outGrp_WSS, 0, maxNSD, "WSS") }, + {OutputNameType::out_MBF, std::make_tuple(OutputNameType::outGrp_Y, 0, 1, "MBF")}, + {OutputNameType::out_mbfFlux, std::make_tuple(OutputNameType::outGrp_mbfFlx, 0, nsd, "MBF_flux")} }; diff --git a/Code/Source/solver/txt.cpp b/Code/Source/solver/txt.cpp index 89f03c907..287e56c6d 100644 --- a/Code/Source/solver/txt.cpp +++ b/Code/Source/solver/txt.cpp @@ -321,6 +321,7 @@ void txt(Simulation* simulation, const bool init_write, const SolutionStates& so case OutputNameType::outGrp_divV: case OutputNameType::outGrp_J: case OutputNameType::outGrp_mises: + case OutputNameType::outGrp_mbfFlx: post::all_post(simulation, tmpV, solutions, oGrp, iEq); break; diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index c9c94f476..915c71565 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1163,7 +1163,8 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b case OutputNameType::outGrp_hFlx: case OutputNameType::outGrp_stInv: case OutputNameType::outGrp_vortex: - case OutputNameType::outGrp_Visc: + case OutputNameType::outGrp_Visc: + case OutputNameType::outGrp_mbfFlx: post::post(simulation, msh, tmpV, solutions, oGrp, iEq); for (int a = 0; a < msh.nNo; a++) { int Ac = msh.gN(a); diff --git a/tests/cases/fluid/darcy_cylinder/darcy_validation_plot.png b/tests/cases/fluid/darcy_cylinder/darcy_validation_plot.png new file mode 100644 index 000000000..b5fd098f6 --- /dev/null +++ b/tests/cases/fluid/darcy_cylinder/darcy_validation_plot.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:55c9f391196e0a0ba7d5d09cdade95fc6531e0ad741f4581dab9c2d93c598bee +size 148969 diff --git a/tests/cases/fluid/darcy_cylinder/flow_solving_darcy b/tests/cases/fluid/darcy_cylinder/flow_solving_darcy new file mode 100644 index 000000000..6766e7429 --- /dev/null +++ b/tests/cases/fluid/darcy_cylinder/flow_solving_darcy @@ -0,0 +1,42 @@ +#!/bin/bash + +# Name of your job +#SBATCH --job-name=TEST +#SBATCH --partition=amarsden + +#SBATCH --output=TEST.o%j + +# Specifiy the name of the error file. The %j specifies the job ID +#SBATCH --error=TEST_ERR.e%j + +# The walltime you require for your job +#SBATCH --time=24:00:00 + +# Job priority. Usually leave normal. +#SBATCH --qos=normal + +# Number of nodes you are requesting for your job. You can have 24 processors per node +#SBATCH --nodes=1 + +# Amount of memory you require per node. The default is 4000 MB per node. +#SBATCH --mem=6000 + +# Number of processors per node +#SBATCH --ntasks-per-node=16 + +# Send an email to this address when the job starts and finishes +#SBATCH --mail-user=mmegally@stanford.edu +#SBATCH --mail-type=begin +#SBATCH --mail-type=end + +# Run normal batch commands +module load system +module load gcc +module load openmpi/3.1.2 +module load openblas +module load llvm/4.0.0 +module load x11 +module load mesa + +srun /home/users/mmegally/Perfusion/svMultiPhysics/build_darcy_2/svMultiPhysics-build/bin/svmultiphysics svFSI.xml + diff --git a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/bottom.vtp b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/bottom.vtp new file mode 100644 index 000000000..c98f09ae0 --- /dev/null +++ b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/bottom.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8c47e55f124d935488fb0a6da47cf417e5225fb6ede73b2edbcd1239a2b4f498 +size 83044 diff --git a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/inner_wall.vtp b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/inner_wall.vtp new file mode 100644 index 000000000..c558c1ee7 --- /dev/null +++ b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/inner_wall.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:78b09a14424cd64e35a4be1d67831bfa52ce92d444737649b99eadfefd188fdd +size 92490 diff --git a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/outer_wall.vtp b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/outer_wall.vtp new file mode 100644 index 000000000..0725d76a2 --- /dev/null +++ b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/outer_wall.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b749fd60aa89473e6d31efa576170bc94fb54aabccdabf71ca19b3540c15e189 +size 92130 diff --git a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/top.vtp b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/top.vtp new file mode 100644 index 000000000..5258f8abe --- /dev/null +++ b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:eb329a65a000898e8a6cc49d306cf34dade035e9a37fa7c0f42de683ac571dac +size 83072 diff --git a/tests/cases/fluid/darcy_cylinder/result_020.vtu b/tests/cases/fluid/darcy_cylinder/result_020.vtu new file mode 100644 index 000000000..64bb9c721 --- /dev/null +++ b/tests/cases/fluid/darcy_cylinder/result_020.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6c787705b921939503b9968d3139ffcad41e2d17d94254be9f2487a5e7f830e6 +size 1745017 diff --git a/tests/cases/fluid/darcy_cylinder/svFSI.xml b/tests/cases/fluid/darcy_cylinder/svFSI.xml new file mode 100644 index 000000000..a5974a32d --- /dev/null +++ b/tests/cases/fluid/darcy_cylinder/svFSI.xml @@ -0,0 +1,92 @@ + + + + + 0 + 3 + 20 + 0.001 + 0.50 + STOP_SIM + + 1 + result + 1 + 0 + + 5 + 1 + + 1 + 0 + 1 + + + + + tissue_vol.vtu + + + mesh-surfaces/inner_wall.vtp + + + + mesh-surfaces/outer_wall.vtp + + + + mesh-surfaces/bottom.vtp + + + + mesh-surfaces/top.vtp + + + + + + + false + 2 + 5 + 1e-6 + + 1 + 1 + 1e-11 + 0.0 + + + true + true + + + + + fsils + + 1e-6 + 100 + 50 + + + + Dirichlet + Steady + 1.0 + false + + + + Dirichlet + Steady + 0.0 + false + + + + + + + + diff --git a/tests/cases/fluid/darcy_cylinder/tissue_vol.vtu b/tests/cases/fluid/darcy_cylinder/tissue_vol.vtu new file mode 100644 index 000000000..09862e9e4 --- /dev/null +++ b/tests/cases/fluid/darcy_cylinder/tissue_vol.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:518e93edb3a819af183d35aa3ccb35a7a41d10e2546b3abeb869464d5015946e +size 1448656 From 49c5bf69b4ad323ea748ecac8f143d1fe3cc0aee Mon Sep 17 00:00:00 2001 From: Michael Date: Fri, 21 Aug 2026 16:41:48 -0700 Subject: [PATCH 02/28] Fix function signatures in header file --- Code/Source/solver/load_msh.h | 31 ++-------------------- Code/Source/solver/read_msh.h | 49 ++++++++++------------------------- 2 files changed, 15 insertions(+), 65 deletions(-) diff --git a/Code/Source/solver/load_msh.h b/Code/Source/solver/load_msh.h index 294807aa2..530252710 100644 --- a/Code/Source/solver/load_msh.h +++ b/Code/Source/solver/load_msh.h @@ -1,32 +1,5 @@ -/* Copyright (c) Stanford University, The Regents of the University of California, and others. - * - * All Rights Reserved. - * - * See Copyright-SimVascular.txt for additional details. - * - * Permission is hereby granted, free of charge, to any person obtaining - * a copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject - * to the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS - * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED - * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A - * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER - * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF - * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause #ifndef LOAD_MSH_H #define LOAD_MSH_H diff --git a/Code/Source/solver/read_msh.h b/Code/Source/solver/read_msh.h index 0712da0d3..c8d336eed 100644 --- a/Code/Source/solver/read_msh.h +++ b/Code/Source/solver/read_msh.h @@ -1,37 +1,11 @@ -/* Copyright (c) Stanford University, The Regents of the University of California, and others. - * - * All Rights Reserved. - * - * See Copyright-SimVascular.txt for additional details. - * - * Permission is hereby granted, free of charge, to any person obtaining - * a copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject - * to the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS - * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED - * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A - * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER - * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF - * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause #ifndef READ_MSH_H #define READ_MSH_H #include "ComMod.h" +#include "SolutionStates.h" #include "Simulation.h" #include "Vector.h" @@ -48,11 +22,11 @@ namespace read_msh_ns { Vector gN; }; - void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag); - void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag); - void calc_elem_skew(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag); + void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const SolutionStates& solutions); + void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const SolutionStates& solutions); + void calc_elem_skew(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const SolutionStates& solutions); - void calc_mesh_props(ComMod& com_mod, const CmMod& cm_mod, const int nMesh, std::vector& mesh); + void calc_mesh_props(ComMod& com_mod, const CmMod& cm_mod, const int nMesh, std::vector& mesh, const SolutionStates& solutions); void calc_nbc(mshType& mesh, faceType& face); @@ -72,15 +46,18 @@ namespace read_msh_ns { void load_var_ini(Simulation* simulation, const ComMod& com_mod); void match_faces(const ComMod& com_mod, const faceType& face1, const faceType& face2, const double tol, utils::stackType& lPrj); + void match_nodes(const ComMod& com_mod, const faceType& lFa, const faceType& pFa, + const double ptol, const int nNds, Array& map); void read_fib_nff(Simulation* simulation, mshType& mesh, const std::string& fName, const std::string& kwrd, const int idx); void read_msh(Simulation* simulation); void set_dmn_id_ff(Simulation* simulation, mshType& mesh, const std::string& file_name); - void set_dmn_id_vtk(Simulation* simulation, mshType& mesh, const std::string& file_name, const std::string& kwrd); + void set_dmn_id_vtk(Simulation* simulation, mshType& lM, const std::string& file_name, const std::string& kwrd); void set_projector(Simulation* simulation, utils::stackType& avNds); - - + void set_ris_projector(Simulation* simulation); + void set_uris_meshes(Simulation* simulation); + }; #endif From e2841188ee6007aab2bbf02cd7c817a822d5ba93 Mon Sep 17 00:00:00 2001 From: Michael Date: Tue, 25 Aug 2026 11:49:16 -0700 Subject: [PATCH 03/28] add convergence tests for darcy implementation --- .../coarse-mesh/flow_solving_darcy | 36 ++++++ .../coarse-mesh/mesh-surfaces/bottom.vtp | 3 + .../coarse-mesh/mesh-surfaces/inner_wall.vtp | 3 + .../coarse-mesh/mesh-surfaces/outer_wall.vtp | 3 + .../coarse-mesh/mesh-surfaces/top.vtp | 3 + .../darcy_convergence/coarse-mesh/svFSI.xml | 92 +++++++++++++++ .../coarse-mesh/tissue_vol.vtu | 3 + .../fluid/darcy_convergence/conv_analysis.py | 110 ++++++++++++++++++ .../fine-mesh/flow_solving_darcy | 36 ++++++ .../fine-mesh/mesh-surfaces/bottom.vtp | 3 + .../fine-mesh/mesh-surfaces/inner_wall.vtp | 3 + .../fine-mesh/mesh-surfaces/outer_wall.vtp | 3 + .../fine-mesh/mesh-surfaces/top.vtp | 3 + .../darcy_convergence/fine-mesh/svFSI.xml | 92 +++++++++++++++ .../fine-mesh/tissue_vol.vtu | 3 + .../med-mesh/flow_solving_darcy | 36 ++++++ .../med-mesh/mesh-surfaces/bottom.vtp | 3 + .../med-mesh/mesh-surfaces/inner_wall.vtp | 3 + .../med-mesh/mesh-surfaces/outer_wall.vtp | 3 + .../med-mesh/mesh-surfaces/top.vtp | 3 + .../darcy_convergence/med-mesh/svFSI.xml | 92 +++++++++++++++ .../darcy_convergence/med-mesh/tissue_vol.vtu | 3 + 22 files changed, 539 insertions(+) create mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/flow_solving_darcy create mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/bottom.vtp create mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/inner_wall.vtp create mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/outer_wall.vtp create mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/top.vtp create mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/svFSI.xml create mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/tissue_vol.vtu create mode 100644 tests/cases/fluid/darcy_convergence/conv_analysis.py create mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/flow_solving_darcy create mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/bottom.vtp create mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/inner_wall.vtp create mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/outer_wall.vtp create mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/top.vtp create mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/svFSI.xml create mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/tissue_vol.vtu create mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/flow_solving_darcy create mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/bottom.vtp create mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/inner_wall.vtp create mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/outer_wall.vtp create mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/top.vtp create mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/svFSI.xml create mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/tissue_vol.vtu diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/flow_solving_darcy b/tests/cases/fluid/darcy_convergence/coarse-mesh/flow_solving_darcy new file mode 100644 index 000000000..1c7331020 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/coarse-mesh/flow_solving_darcy @@ -0,0 +1,36 @@ +#!/bin/bash + +# Name of your job +#SBATCH --job-name=TEST +#SBATCH --partition=amarsden + +#SBATCH --output=TEST.o%j + +# Specifiy the name of the error file. The %j specifies the job ID +#SBATCH --error=TEST_ERR.e%j + +# The walltime you require for your job +#SBATCH --time=24:00:00 + +# Job priority. Usually leave normal. +#SBATCH --qos=normal + +# Number of nodes you are requesting for your job. You can have 24 processors per node +#SBATCH --nodes=1 + +# Amount of memory you require per node. The default is 4000 MB per node. +#SBATCH --mem=6000 + +# Number of processors per node +#SBATCH --ntasks-per-node=16 + +# Send an email to this address when the job starts and finishes +#SBATCH --mail-user=mmegally@stanford.edu +#SBATCH --mail-type=begin +#SBATCH --mail-type=end + +# Run normal batch commands +module load system viz gcc/10.1.0 openmpi/4.1.2 vtk/9.4.1 openblas/0.3.10 + +srun /home/users/mmegally/darcy_build/build/svMultiPhysics-build/bin/svmultiphysics svFSI.xml + diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/bottom.vtp b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/bottom.vtp new file mode 100644 index 000000000..1f4334702 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/bottom.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0bf1e60b12de03cf56316d42515c6b26cded32268b31cfc63cf896b05a9f9c1f +size 13638 diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/inner_wall.vtp b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/inner_wall.vtp new file mode 100644 index 000000000..17bc528a9 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/inner_wall.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ac60179c80b16b98eb4c59b705c28b0a90363f7ace3a1e9a6f31c02470d71b82 +size 14993 diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/outer_wall.vtp b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/outer_wall.vtp new file mode 100644 index 000000000..724e6d591 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/outer_wall.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bd551041810e4d97f715f0a17f1515d61e43a2bb182dd21d3b91b7c350a334d4 +size 14984 diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/top.vtp b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/top.vtp new file mode 100644 index 000000000..0151f335d --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e254917ca313d3a2130a0a14ab735d6cca682f62c9b76d28901a43ebd8ef8bb0 +size 13640 diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/svFSI.xml b/tests/cases/fluid/darcy_convergence/coarse-mesh/svFSI.xml new file mode 100644 index 000000000..a5974a32d --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/coarse-mesh/svFSI.xml @@ -0,0 +1,92 @@ + + + + + 0 + 3 + 20 + 0.001 + 0.50 + STOP_SIM + + 1 + result + 1 + 0 + + 5 + 1 + + 1 + 0 + 1 + + + + + tissue_vol.vtu + + + mesh-surfaces/inner_wall.vtp + + + + mesh-surfaces/outer_wall.vtp + + + + mesh-surfaces/bottom.vtp + + + + mesh-surfaces/top.vtp + + + + + + + false + 2 + 5 + 1e-6 + + 1 + 1 + 1e-11 + 0.0 + + + true + true + + + + + fsils + + 1e-6 + 100 + 50 + + + + Dirichlet + Steady + 1.0 + false + + + + Dirichlet + Steady + 0.0 + false + + + + + + + + diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/tissue_vol.vtu b/tests/cases/fluid/darcy_convergence/coarse-mesh/tissue_vol.vtu new file mode 100644 index 000000000..cd5a6a4b5 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/coarse-mesh/tissue_vol.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4475a6e60f24b8966911571d14af9dc145bf66355ff5345905a6f4feb80deebb +size 108738 diff --git a/tests/cases/fluid/darcy_convergence/conv_analysis.py b/tests/cases/fluid/darcy_convergence/conv_analysis.py new file mode 100644 index 000000000..343337a92 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/conv_analysis.py @@ -0,0 +1,110 @@ +import pyvista as pv +import numpy as np +import matplotlib.pyplot as plt + + +# Analytical Solutions: +def exact_pressure(x, y, z): + r = np.sqrt(x ** 2 + y ** 2) + return 1.0 - (np.log(r) / np.log(2.0)) + + +def exact_velocity(x, y, z): + K = 1e-11 # Permeability + mu = 1.0 # Viscosity + + r = np.sqrt(x ** 2 + y ** 2) + + # Apply the Darcy material scale to the pressure gradient + vel_r = (K / mu) * (1.0 / (r * np.log(2.0))) + + theta = np.arctan2(y, x) + u_x = vel_r * np.cos(theta) + u_y = vel_r * np.sin(theta) + u_z = np.zeros_like(r) + return np.column_stack((u_x, u_y, u_z)) + + + +# Compute L2 Norms: +def compute_l2_errors(result_file, pressure_array="MBF", velocity_array="MBF_flux"): + mesh = pv.read(result_file) + mesh = mesh.compute_cell_sizes() + volumes = np.abs(mesh.cell_data["Volume"]) + + centers = mesh.cell_centers().points + x, y, z = centers[:, 0], centers[:, 1], centers[:, 2] + + if pressure_array in mesh.point_data: + mesh = mesh.ptc() + + num_p = mesh.cell_data[pressure_array] + num_v = mesh.cell_data[velocity_array] + + ana_p = exact_pressure(x, y, z) + ana_v = exact_velocity(x, y, z) + + p_err_sq = (num_p - ana_p) ** 2 + L2_p = np.sqrt(np.sum(p_err_sq * volumes)) + + v_err_sq = np.sum((num_v - ana_v) ** 2, axis=1) + L2_v = np.sqrt(np.sum(v_err_sq * volumes)) + + total_volume = np.sum(volumes) + h = (total_volume / mesh.n_cells) ** (1 / 3.0) + + return h, L2_p, L2_v + + +if __name__ == "__main__": + base_dir = "." + meshes = [ + f"{base_dir}/coarse-mesh/16-procs/result_020.vtu", + f"{base_dir}/med-mesh/16-procs/result_020.vtu", + f"{base_dir}/fine-mesh/16-procs/result_020.vtu" + ] + + h_vals, p_errs, v_errs = [], [], [] + + print(f"{'Mesh Level':<15} | {'h (Elem Size)':<15} | {'L2 Pressure':<15} | {'L2 Velocity':<15}") + print("-" * 65) + + for mesh_file in meshes: + mesh_name = mesh_file.split('/')[1] + h, L2_p, L2_v = compute_l2_errors(mesh_file) + h_vals.append(h) + p_errs.append(L2_p) + v_errs.append(L2_v) + print(f"{mesh_name:<15} | {h:<15.6e} | {L2_p:<15.6e} | {L2_v:<15.6e}") + + + # Plot the Convergence Rates + h_vals = np.array(h_vals) + p_errs = np.array(p_errs) + v_errs = np.array(v_errs) + + # Create theoretical reference lines starting from the coarse mesh error + p_ref = p_errs[0] * (h_vals / h_vals[0]) ** 2 # O(h^2) slope + v_ref = v_errs[0] * (h_vals / h_vals[0]) ** 1 # O(h^1) slope + + plt.figure(figsize=(9, 7)) + + # Plot actual calculated errors + plt.loglog(h_vals, p_errs, 'o-', linewidth=2, markersize=8, label='Pressure $L_2$ Error', color='blue') + plt.loglog(h_vals, v_errs, 's-', linewidth=2, markersize=8, label='Velocity $L_2$ Error', color='red') + + # Plot reference slopes (dashed) + plt.loglog(h_vals, p_ref, '--', color='lightblue', label='Expected $O(h^2)$') + plt.loglog(h_vals, v_ref, '--', color='lightcoral', label='Expected $O(h)$') + + plt.xlabel('Element Size ($h$)', fontsize=12) + plt.ylabel('$L_2$ Norm Error', fontsize=12) + plt.title('Mesh Convergence: Darcy Flow', fontsize=14) + plt.grid(True, which="both", ls="--", alpha=0.5) + plt.legend(fontsize=11) + + # Invert X-axis so finer meshes (smaller h) are on the right + plt.gca().invert_xaxis() + + plt.tight_layout() + plt.show() diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/flow_solving_darcy b/tests/cases/fluid/darcy_convergence/fine-mesh/flow_solving_darcy new file mode 100644 index 000000000..1c7331020 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/fine-mesh/flow_solving_darcy @@ -0,0 +1,36 @@ +#!/bin/bash + +# Name of your job +#SBATCH --job-name=TEST +#SBATCH --partition=amarsden + +#SBATCH --output=TEST.o%j + +# Specifiy the name of the error file. The %j specifies the job ID +#SBATCH --error=TEST_ERR.e%j + +# The walltime you require for your job +#SBATCH --time=24:00:00 + +# Job priority. Usually leave normal. +#SBATCH --qos=normal + +# Number of nodes you are requesting for your job. You can have 24 processors per node +#SBATCH --nodes=1 + +# Amount of memory you require per node. The default is 4000 MB per node. +#SBATCH --mem=6000 + +# Number of processors per node +#SBATCH --ntasks-per-node=16 + +# Send an email to this address when the job starts and finishes +#SBATCH --mail-user=mmegally@stanford.edu +#SBATCH --mail-type=begin +#SBATCH --mail-type=end + +# Run normal batch commands +module load system viz gcc/10.1.0 openmpi/4.1.2 vtk/9.4.1 openblas/0.3.10 + +srun /home/users/mmegally/darcy_build/build/svMultiPhysics-build/bin/svmultiphysics svFSI.xml + diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/bottom.vtp b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/bottom.vtp new file mode 100644 index 000000000..ddda62b1c --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/bottom.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f23c63eb33a679d161b2588d901529b68c845725717d9547d77695bf61c02139 +size 382399 diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/inner_wall.vtp b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/inner_wall.vtp new file mode 100644 index 000000000..5cb89de59 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/inner_wall.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b9aeeeb1e1e6905d094524c86961e0be1a18bb19f705ba49a826a75656642fe6 +size 389434 diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/outer_wall.vtp b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/outer_wall.vtp new file mode 100644 index 000000000..6ec31850c --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/outer_wall.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f3ca33b843158a517a47ef3d8abc38a86b9356993beff11d696a6d24f19f8e9e +size 388018 diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/top.vtp b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/top.vtp new file mode 100644 index 000000000..567d3b5c5 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e0bf4150e7937e82039433eceb8a772241077639f2446292e1e6143528a9bf26 +size 381992 diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/svFSI.xml b/tests/cases/fluid/darcy_convergence/fine-mesh/svFSI.xml new file mode 100644 index 000000000..a5974a32d --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/fine-mesh/svFSI.xml @@ -0,0 +1,92 @@ + + + + + 0 + 3 + 20 + 0.001 + 0.50 + STOP_SIM + + 1 + result + 1 + 0 + + 5 + 1 + + 1 + 0 + 1 + + + + + tissue_vol.vtu + + + mesh-surfaces/inner_wall.vtp + + + + mesh-surfaces/outer_wall.vtp + + + + mesh-surfaces/bottom.vtp + + + + mesh-surfaces/top.vtp + + + + + + + false + 2 + 5 + 1e-6 + + 1 + 1 + 1e-11 + 0.0 + + + true + true + + + + + fsils + + 1e-6 + 100 + 50 + + + + Dirichlet + Steady + 1.0 + false + + + + Dirichlet + Steady + 0.0 + false + + + + + + + + diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/tissue_vol.vtu b/tests/cases/fluid/darcy_convergence/fine-mesh/tissue_vol.vtu new file mode 100644 index 000000000..8f03b4a00 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/fine-mesh/tissue_vol.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:08c08fc4a71a6170e5129d82d0939a94781d4f39a7d8987077af53eaa7fd6edb +size 9690906 diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/flow_solving_darcy b/tests/cases/fluid/darcy_convergence/med-mesh/flow_solving_darcy new file mode 100644 index 000000000..1c7331020 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/med-mesh/flow_solving_darcy @@ -0,0 +1,36 @@ +#!/bin/bash + +# Name of your job +#SBATCH --job-name=TEST +#SBATCH --partition=amarsden + +#SBATCH --output=TEST.o%j + +# Specifiy the name of the error file. The %j specifies the job ID +#SBATCH --error=TEST_ERR.e%j + +# The walltime you require for your job +#SBATCH --time=24:00:00 + +# Job priority. Usually leave normal. +#SBATCH --qos=normal + +# Number of nodes you are requesting for your job. You can have 24 processors per node +#SBATCH --nodes=1 + +# Amount of memory you require per node. The default is 4000 MB per node. +#SBATCH --mem=6000 + +# Number of processors per node +#SBATCH --ntasks-per-node=16 + +# Send an email to this address when the job starts and finishes +#SBATCH --mail-user=mmegally@stanford.edu +#SBATCH --mail-type=begin +#SBATCH --mail-type=end + +# Run normal batch commands +module load system viz gcc/10.1.0 openmpi/4.1.2 vtk/9.4.1 openblas/0.3.10 + +srun /home/users/mmegally/darcy_build/build/svMultiPhysics-build/bin/svmultiphysics svFSI.xml + diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/bottom.vtp b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/bottom.vtp new file mode 100644 index 000000000..3d1ff6632 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/bottom.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ca7129083bd1a7d2ba8e56fecb0b4e0c83e03c4b90dc71770a01135d000ad3e8 +size 61132 diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/inner_wall.vtp b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/inner_wall.vtp new file mode 100644 index 000000000..57cb3ac98 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/inner_wall.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:29535ba2be5b2ce0956eb9847e8f306ecd7bb34a8cb856c71ee22afedf924b6d +size 65378 diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/outer_wall.vtp b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/outer_wall.vtp new file mode 100644 index 000000000..38abe1060 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/outer_wall.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:69bfbca17619c9288260303cfa6639c35163bd0e73897d78d7398611152dfbb9 +size 65194 diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/top.vtp b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/top.vtp new file mode 100644 index 000000000..2da7351b5 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e7e31b71e8e44475ad02880a745b82aee444868e031fd29aa939f6962f2a468a +size 61108 diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/svFSI.xml b/tests/cases/fluid/darcy_convergence/med-mesh/svFSI.xml new file mode 100644 index 000000000..a5974a32d --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/med-mesh/svFSI.xml @@ -0,0 +1,92 @@ + + + + + 0 + 3 + 20 + 0.001 + 0.50 + STOP_SIM + + 1 + result + 1 + 0 + + 5 + 1 + + 1 + 0 + 1 + + + + + tissue_vol.vtu + + + mesh-surfaces/inner_wall.vtp + + + + mesh-surfaces/outer_wall.vtp + + + + mesh-surfaces/bottom.vtp + + + + mesh-surfaces/top.vtp + + + + + + + false + 2 + 5 + 1e-6 + + 1 + 1 + 1e-11 + 0.0 + + + true + true + + + + + fsils + + 1e-6 + 100 + 50 + + + + Dirichlet + Steady + 1.0 + false + + + + Dirichlet + Steady + 0.0 + false + + + + + + + + diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/tissue_vol.vtu b/tests/cases/fluid/darcy_convergence/med-mesh/tissue_vol.vtu new file mode 100644 index 000000000..98ab74186 --- /dev/null +++ b/tests/cases/fluid/darcy_convergence/med-mesh/tissue_vol.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:de64d090aed9534ff5bd50a438e25d7e90d2d1f89558568996c8378f0cb8c3bf +size 1032350 From b1cd37315e23bb2852ad9dd0c8cfc6f5bc11138e Mon Sep 17 00:00:00 2001 From: Michael Date: Wed, 26 Aug 2026 15:12:02 -0700 Subject: [PATCH 04/28] Add Regazzoni model to validate Ubuntu/macOS integration tests (added after fork) --- Code/Source/solver/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index ddb614b29..58a41d667 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -170,6 +170,7 @@ set(CSRCS SimulationLogger.h VtkData.h VtkData.cpp + active_stress_regazzoni.h active_stress_regazzoni.cpp all_fun.h all_fun.cpp baf_ini.h baf_ini.cpp bf.h bf.cpp From bae4ba1eeaf5b211c12b9ffeb453a3a9d090e5f6 Mon Sep 17 00:00:00 2001 From: Michael Date: Wed, 26 Aug 2026 16:03:38 -0700 Subject: [PATCH 05/28] Remove unintended SLURM files --- .../coarse-mesh/flow_solving_darcy | 36 ---------------- .../fine-mesh/flow_solving_darcy | 36 ---------------- .../med-mesh/flow_solving_darcy | 36 ---------------- .../fluid/darcy_cylinder/flow_solving_darcy | 42 ------------------- 4 files changed, 150 deletions(-) delete mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/flow_solving_darcy delete mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/flow_solving_darcy delete mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/flow_solving_darcy delete mode 100644 tests/cases/fluid/darcy_cylinder/flow_solving_darcy diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/flow_solving_darcy b/tests/cases/fluid/darcy_convergence/coarse-mesh/flow_solving_darcy deleted file mode 100644 index 1c7331020..000000000 --- a/tests/cases/fluid/darcy_convergence/coarse-mesh/flow_solving_darcy +++ /dev/null @@ -1,36 +0,0 @@ -#!/bin/bash - -# Name of your job -#SBATCH --job-name=TEST -#SBATCH --partition=amarsden - -#SBATCH --output=TEST.o%j - -# Specifiy the name of the error file. The %j specifies the job ID -#SBATCH --error=TEST_ERR.e%j - -# The walltime you require for your job -#SBATCH --time=24:00:00 - -# Job priority. Usually leave normal. -#SBATCH --qos=normal - -# Number of nodes you are requesting for your job. You can have 24 processors per node -#SBATCH --nodes=1 - -# Amount of memory you require per node. The default is 4000 MB per node. -#SBATCH --mem=6000 - -# Number of processors per node -#SBATCH --ntasks-per-node=16 - -# Send an email to this address when the job starts and finishes -#SBATCH --mail-user=mmegally@stanford.edu -#SBATCH --mail-type=begin -#SBATCH --mail-type=end - -# Run normal batch commands -module load system viz gcc/10.1.0 openmpi/4.1.2 vtk/9.4.1 openblas/0.3.10 - -srun /home/users/mmegally/darcy_build/build/svMultiPhysics-build/bin/svmultiphysics svFSI.xml - diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/flow_solving_darcy b/tests/cases/fluid/darcy_convergence/fine-mesh/flow_solving_darcy deleted file mode 100644 index 1c7331020..000000000 --- a/tests/cases/fluid/darcy_convergence/fine-mesh/flow_solving_darcy +++ /dev/null @@ -1,36 +0,0 @@ -#!/bin/bash - -# Name of your job -#SBATCH --job-name=TEST -#SBATCH --partition=amarsden - -#SBATCH --output=TEST.o%j - -# Specifiy the name of the error file. The %j specifies the job ID -#SBATCH --error=TEST_ERR.e%j - -# The walltime you require for your job -#SBATCH --time=24:00:00 - -# Job priority. Usually leave normal. -#SBATCH --qos=normal - -# Number of nodes you are requesting for your job. You can have 24 processors per node -#SBATCH --nodes=1 - -# Amount of memory you require per node. The default is 4000 MB per node. -#SBATCH --mem=6000 - -# Number of processors per node -#SBATCH --ntasks-per-node=16 - -# Send an email to this address when the job starts and finishes -#SBATCH --mail-user=mmegally@stanford.edu -#SBATCH --mail-type=begin -#SBATCH --mail-type=end - -# Run normal batch commands -module load system viz gcc/10.1.0 openmpi/4.1.2 vtk/9.4.1 openblas/0.3.10 - -srun /home/users/mmegally/darcy_build/build/svMultiPhysics-build/bin/svmultiphysics svFSI.xml - diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/flow_solving_darcy b/tests/cases/fluid/darcy_convergence/med-mesh/flow_solving_darcy deleted file mode 100644 index 1c7331020..000000000 --- a/tests/cases/fluid/darcy_convergence/med-mesh/flow_solving_darcy +++ /dev/null @@ -1,36 +0,0 @@ -#!/bin/bash - -# Name of your job -#SBATCH --job-name=TEST -#SBATCH --partition=amarsden - -#SBATCH --output=TEST.o%j - -# Specifiy the name of the error file. The %j specifies the job ID -#SBATCH --error=TEST_ERR.e%j - -# The walltime you require for your job -#SBATCH --time=24:00:00 - -# Job priority. Usually leave normal. -#SBATCH --qos=normal - -# Number of nodes you are requesting for your job. You can have 24 processors per node -#SBATCH --nodes=1 - -# Amount of memory you require per node. The default is 4000 MB per node. -#SBATCH --mem=6000 - -# Number of processors per node -#SBATCH --ntasks-per-node=16 - -# Send an email to this address when the job starts and finishes -#SBATCH --mail-user=mmegally@stanford.edu -#SBATCH --mail-type=begin -#SBATCH --mail-type=end - -# Run normal batch commands -module load system viz gcc/10.1.0 openmpi/4.1.2 vtk/9.4.1 openblas/0.3.10 - -srun /home/users/mmegally/darcy_build/build/svMultiPhysics-build/bin/svmultiphysics svFSI.xml - diff --git a/tests/cases/fluid/darcy_cylinder/flow_solving_darcy b/tests/cases/fluid/darcy_cylinder/flow_solving_darcy deleted file mode 100644 index 6766e7429..000000000 --- a/tests/cases/fluid/darcy_cylinder/flow_solving_darcy +++ /dev/null @@ -1,42 +0,0 @@ -#!/bin/bash - -# Name of your job -#SBATCH --job-name=TEST -#SBATCH --partition=amarsden - -#SBATCH --output=TEST.o%j - -# Specifiy the name of the error file. The %j specifies the job ID -#SBATCH --error=TEST_ERR.e%j - -# The walltime you require for your job -#SBATCH --time=24:00:00 - -# Job priority. Usually leave normal. -#SBATCH --qos=normal - -# Number of nodes you are requesting for your job. You can have 24 processors per node -#SBATCH --nodes=1 - -# Amount of memory you require per node. The default is 4000 MB per node. -#SBATCH --mem=6000 - -# Number of processors per node -#SBATCH --ntasks-per-node=16 - -# Send an email to this address when the job starts and finishes -#SBATCH --mail-user=mmegally@stanford.edu -#SBATCH --mail-type=begin -#SBATCH --mail-type=end - -# Run normal batch commands -module load system -module load gcc -module load openmpi/3.1.2 -module load openblas -module load llvm/4.0.0 -module load x11 -module load mesa - -srun /home/users/mmegally/Perfusion/svMultiPhysics/build_darcy_2/svMultiPhysics-build/bin/svmultiphysics svFSI.xml - From e65e247e8167eddd30bf1177c5d29ea91e7f7cd0 Mon Sep 17 00:00:00 2001 From: Michael Date: Thu, 27 Aug 2026 13:36:17 -0700 Subject: [PATCH 06/28] Addressing pull request comments (removed tests from CI suite and Zack's suggested changes) --- Code/Source/solver/Parameters.cpp | 4 - Code/Source/solver/Parameters.h | 4 - Code/Source/solver/consts.h | 9 +- Code/Source/solver/darcy.cpp | 30 +++-- Code/Source/solver/initialize.cpp | 6 +- Code/Source/solver/post.cpp | 53 +++++++-- Code/Source/solver/read_files.cpp | 18 +-- Code/Source/solver/set_equation_props.h | 11 +- .../coarse-mesh/mesh-surfaces/bottom.vtp | 3 - .../coarse-mesh/mesh-surfaces/inner_wall.vtp | 3 - .../coarse-mesh/mesh-surfaces/outer_wall.vtp | 3 - .../coarse-mesh/mesh-surfaces/top.vtp | 3 - .../darcy_convergence/coarse-mesh/svFSI.xml | 92 --------------- .../coarse-mesh/tissue_vol.vtu | 3 - .../fluid/darcy_convergence/conv_analysis.py | 110 ------------------ .../fine-mesh/mesh-surfaces/bottom.vtp | 3 - .../fine-mesh/mesh-surfaces/inner_wall.vtp | 3 - .../fine-mesh/mesh-surfaces/outer_wall.vtp | 3 - .../fine-mesh/mesh-surfaces/top.vtp | 3 - .../darcy_convergence/fine-mesh/svFSI.xml | 92 --------------- .../fine-mesh/tissue_vol.vtu | 3 - .../med-mesh/mesh-surfaces/bottom.vtp | 3 - .../med-mesh/mesh-surfaces/inner_wall.vtp | 3 - .../med-mesh/mesh-surfaces/outer_wall.vtp | 3 - .../med-mesh/mesh-surfaces/top.vtp | 3 - .../darcy_convergence/med-mesh/svFSI.xml | 92 --------------- .../darcy_convergence/med-mesh/tissue_vol.vtu | 3 - .../darcy_cylinder/darcy_validation_plot.png | 3 - .../darcy_cylinder/mesh-surfaces/bottom.vtp | 3 - .../mesh-surfaces/inner_wall.vtp | 3 - .../mesh-surfaces/outer_wall.vtp | 3 - .../darcy_cylinder/mesh-surfaces/top.vtp | 3 - .../cases/fluid/darcy_cylinder/result_020.vtu | 3 - tests/cases/fluid/darcy_cylinder/svFSI.xml | 92 --------------- .../cases/fluid/darcy_cylinder/tissue_vol.vtu | 3 - 35 files changed, 72 insertions(+), 607 deletions(-) delete mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/bottom.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/inner_wall.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/outer_wall.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/top.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/svFSI.xml delete mode 100644 tests/cases/fluid/darcy_convergence/coarse-mesh/tissue_vol.vtu delete mode 100644 tests/cases/fluid/darcy_convergence/conv_analysis.py delete mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/bottom.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/inner_wall.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/outer_wall.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/top.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/svFSI.xml delete mode 100644 tests/cases/fluid/darcy_convergence/fine-mesh/tissue_vol.vtu delete mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/bottom.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/inner_wall.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/outer_wall.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/top.vtp delete mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/svFSI.xml delete mode 100644 tests/cases/fluid/darcy_convergence/med-mesh/tissue_vol.vtu delete mode 100644 tests/cases/fluid/darcy_cylinder/darcy_validation_plot.png delete mode 100644 tests/cases/fluid/darcy_cylinder/mesh-surfaces/bottom.vtp delete mode 100644 tests/cases/fluid/darcy_cylinder/mesh-surfaces/inner_wall.vtp delete mode 100644 tests/cases/fluid/darcy_cylinder/mesh-surfaces/outer_wall.vtp delete mode 100644 tests/cases/fluid/darcy_cylinder/mesh-surfaces/top.vtp delete mode 100644 tests/cases/fluid/darcy_cylinder/result_020.vtu delete mode 100644 tests/cases/fluid/darcy_cylinder/svFSI.xml delete mode 100644 tests/cases/fluid/darcy_cylinder/tissue_vol.vtu diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 0be19fd1c..fdcfadc23 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2046,12 +2046,8 @@ DomainParameters::DomainParameters() { set_parameter("Poisson_ratio", 0.3, !required, poisson_ratio); 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); set_parameter("Relative_tolerance", 1e-4, !required, relative_tolerance); set_parameter("Shell_thickness", 0.0, !required, shell_thickness); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 325b0ce08..29445c859 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1624,12 +1624,8 @@ class DomainParameters : public ParameterLists Parameter time_step_for_integration; Parameter permeability; - Parameter porosity; - Parameter porosity_pressure; Parameter media_compressibility; - Parameter fluid_compressibility; Parameter darcy_fluid_viscosity; - Parameter density_pressure; // Inverse of Darcy permeability. Default value of 0.0 for Navier-Stokes and non-zero for Navier-Stokes-Brinkman Parameter inverse_darcy_permeability; }; diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index c8817548e..4f042557a 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -291,6 +291,7 @@ enum class EquationType 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; @@ -419,12 +420,8 @@ enum class PhysicalProperyType ctau_C = 14, inverse_darcy_permeability = 15, permeability = 16, - porosity = 17, - porosity_pressure = 18, - media_compressibility = 19, - fluid_compressibility = 20, - darcy_fluid_viscosity = 21, - density_pressure = 22 + media_compressibility = 17, + darcy_fluid_viscosity = 18 }; enum class PreconditionerType diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 4abf8c5a3..01fbfec91 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -161,17 +161,18 @@ void darcy_1d(ComMod& com_mod, const int eNoN, const double w, const Vector 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) + @@ -273,19 +276,22 @@ void darcy_3d(ComMod& com_mod, const int eNoN, const double w, const Vector 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 + ((k*rho_0)/mu)*(Nx(0,a)*Nx(0,b) + Nx(1,a)*Nx(1,b) + diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 7baf4f330..f59513013 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -440,9 +440,9 @@ void initialize(Simulation* simulation, Vector& timeP) // std::tie(eq.dof, eq.sym) = equation_dof_map.at(eq.phys); - if (std::set{Equation_fluid, Equation_heatF, Equation_heatS, Equation_CEP, Equation_stokes}.count(eq.phys) == 0) { - dFlag = true; - } + if (std::set{Equation_CEP, Equation_darcy, Equation_fluid, Equation_heatF, Equation_heatS, Equation_stokes}.count(eq.phys) == 0) { + dFlag = true; + } // For second order eqs. if (std::set{Equation_lElas, Equation_struct, Equation_shell, Equation_mesh}.count(eq.phys) != 0) { diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 9daccf7f1..14d04eeac 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -916,11 +916,29 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S Array ksix(nsd,nsd); double Jac = 0.0; + Vector 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; + } } double w = lM.w(g) * Jac; @@ -1014,17 +1032,34 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S // MBF Flux calculation // } else if (outGrp == OutputNameType::outGrp_mbfFlx) { - double kappa = eq.dmn[cDmn].prop[PhysicalProperyType::permeability]; - int i = eq.s; - Vector q(nsd); - for (int a = 0; a < eNoN; 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 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++) { - q(j) = q(j) + Nx(j, a) * yl(i, a); + grad_p(j) = dp_ds * fiber_tangent(j); } - } - for (int j = 0; j < nsd; j++) { - lRes(j) = -kappa * q(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) = -mobility * grad_p(j); + } // Strain tensor invariants calculation // diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index f1686b493..ff069516a 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1550,7 +1550,7 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& break; 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(); @@ -1585,29 +1585,13 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& 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; } lEq.dmn[iDmn].prop[prop] = rtmp; diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index aa492d44b..27e8d003d 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -286,14 +286,9 @@ SetEquationPropertiesMapType set_equation_props = { 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; read_domain(simulation, eq_params, lEq, propL); diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/bottom.vtp b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/bottom.vtp deleted file mode 100644 index 1f4334702..000000000 --- a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/bottom.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:0bf1e60b12de03cf56316d42515c6b26cded32268b31cfc63cf896b05a9f9c1f -size 13638 diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/inner_wall.vtp b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/inner_wall.vtp deleted file mode 100644 index 17bc528a9..000000000 --- a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/inner_wall.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:ac60179c80b16b98eb4c59b705c28b0a90363f7ace3a1e9a6f31c02470d71b82 -size 14993 diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/outer_wall.vtp b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/outer_wall.vtp deleted file mode 100644 index 724e6d591..000000000 --- a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/outer_wall.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:bd551041810e4d97f715f0a17f1515d61e43a2bb182dd21d3b91b7c350a334d4 -size 14984 diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/top.vtp b/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/top.vtp deleted file mode 100644 index 0151f335d..000000000 --- a/tests/cases/fluid/darcy_convergence/coarse-mesh/mesh-surfaces/top.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:e254917ca313d3a2130a0a14ab735d6cca682f62c9b76d28901a43ebd8ef8bb0 -size 13640 diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/svFSI.xml b/tests/cases/fluid/darcy_convergence/coarse-mesh/svFSI.xml deleted file mode 100644 index a5974a32d..000000000 --- a/tests/cases/fluid/darcy_convergence/coarse-mesh/svFSI.xml +++ /dev/null @@ -1,92 +0,0 @@ - - - - - 0 - 3 - 20 - 0.001 - 0.50 - STOP_SIM - - 1 - result - 1 - 0 - - 5 - 1 - - 1 - 0 - 1 - - - - - tissue_vol.vtu - - - mesh-surfaces/inner_wall.vtp - - - - mesh-surfaces/outer_wall.vtp - - - - mesh-surfaces/bottom.vtp - - - - mesh-surfaces/top.vtp - - - - - - - false - 2 - 5 - 1e-6 - - 1 - 1 - 1e-11 - 0.0 - - - true - true - - - - - fsils - - 1e-6 - 100 - 50 - - - - Dirichlet - Steady - 1.0 - false - - - - Dirichlet - Steady - 0.0 - false - - - - - - - - diff --git a/tests/cases/fluid/darcy_convergence/coarse-mesh/tissue_vol.vtu b/tests/cases/fluid/darcy_convergence/coarse-mesh/tissue_vol.vtu deleted file mode 100644 index cd5a6a4b5..000000000 --- a/tests/cases/fluid/darcy_convergence/coarse-mesh/tissue_vol.vtu +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:4475a6e60f24b8966911571d14af9dc145bf66355ff5345905a6f4feb80deebb -size 108738 diff --git a/tests/cases/fluid/darcy_convergence/conv_analysis.py b/tests/cases/fluid/darcy_convergence/conv_analysis.py deleted file mode 100644 index 343337a92..000000000 --- a/tests/cases/fluid/darcy_convergence/conv_analysis.py +++ /dev/null @@ -1,110 +0,0 @@ -import pyvista as pv -import numpy as np -import matplotlib.pyplot as plt - - -# Analytical Solutions: -def exact_pressure(x, y, z): - r = np.sqrt(x ** 2 + y ** 2) - return 1.0 - (np.log(r) / np.log(2.0)) - - -def exact_velocity(x, y, z): - K = 1e-11 # Permeability - mu = 1.0 # Viscosity - - r = np.sqrt(x ** 2 + y ** 2) - - # Apply the Darcy material scale to the pressure gradient - vel_r = (K / mu) * (1.0 / (r * np.log(2.0))) - - theta = np.arctan2(y, x) - u_x = vel_r * np.cos(theta) - u_y = vel_r * np.sin(theta) - u_z = np.zeros_like(r) - return np.column_stack((u_x, u_y, u_z)) - - - -# Compute L2 Norms: -def compute_l2_errors(result_file, pressure_array="MBF", velocity_array="MBF_flux"): - mesh = pv.read(result_file) - mesh = mesh.compute_cell_sizes() - volumes = np.abs(mesh.cell_data["Volume"]) - - centers = mesh.cell_centers().points - x, y, z = centers[:, 0], centers[:, 1], centers[:, 2] - - if pressure_array in mesh.point_data: - mesh = mesh.ptc() - - num_p = mesh.cell_data[pressure_array] - num_v = mesh.cell_data[velocity_array] - - ana_p = exact_pressure(x, y, z) - ana_v = exact_velocity(x, y, z) - - p_err_sq = (num_p - ana_p) ** 2 - L2_p = np.sqrt(np.sum(p_err_sq * volumes)) - - v_err_sq = np.sum((num_v - ana_v) ** 2, axis=1) - L2_v = np.sqrt(np.sum(v_err_sq * volumes)) - - total_volume = np.sum(volumes) - h = (total_volume / mesh.n_cells) ** (1 / 3.0) - - return h, L2_p, L2_v - - -if __name__ == "__main__": - base_dir = "." - meshes = [ - f"{base_dir}/coarse-mesh/16-procs/result_020.vtu", - f"{base_dir}/med-mesh/16-procs/result_020.vtu", - f"{base_dir}/fine-mesh/16-procs/result_020.vtu" - ] - - h_vals, p_errs, v_errs = [], [], [] - - print(f"{'Mesh Level':<15} | {'h (Elem Size)':<15} | {'L2 Pressure':<15} | {'L2 Velocity':<15}") - print("-" * 65) - - for mesh_file in meshes: - mesh_name = mesh_file.split('/')[1] - h, L2_p, L2_v = compute_l2_errors(mesh_file) - h_vals.append(h) - p_errs.append(L2_p) - v_errs.append(L2_v) - print(f"{mesh_name:<15} | {h:<15.6e} | {L2_p:<15.6e} | {L2_v:<15.6e}") - - - # Plot the Convergence Rates - h_vals = np.array(h_vals) - p_errs = np.array(p_errs) - v_errs = np.array(v_errs) - - # Create theoretical reference lines starting from the coarse mesh error - p_ref = p_errs[0] * (h_vals / h_vals[0]) ** 2 # O(h^2) slope - v_ref = v_errs[0] * (h_vals / h_vals[0]) ** 1 # O(h^1) slope - - plt.figure(figsize=(9, 7)) - - # Plot actual calculated errors - plt.loglog(h_vals, p_errs, 'o-', linewidth=2, markersize=8, label='Pressure $L_2$ Error', color='blue') - plt.loglog(h_vals, v_errs, 's-', linewidth=2, markersize=8, label='Velocity $L_2$ Error', color='red') - - # Plot reference slopes (dashed) - plt.loglog(h_vals, p_ref, '--', color='lightblue', label='Expected $O(h^2)$') - plt.loglog(h_vals, v_ref, '--', color='lightcoral', label='Expected $O(h)$') - - plt.xlabel('Element Size ($h$)', fontsize=12) - plt.ylabel('$L_2$ Norm Error', fontsize=12) - plt.title('Mesh Convergence: Darcy Flow', fontsize=14) - plt.grid(True, which="both", ls="--", alpha=0.5) - plt.legend(fontsize=11) - - # Invert X-axis so finer meshes (smaller h) are on the right - plt.gca().invert_xaxis() - - plt.tight_layout() - plt.show() diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/bottom.vtp b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/bottom.vtp deleted file mode 100644 index ddda62b1c..000000000 --- a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/bottom.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:f23c63eb33a679d161b2588d901529b68c845725717d9547d77695bf61c02139 -size 382399 diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/inner_wall.vtp b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/inner_wall.vtp deleted file mode 100644 index 5cb89de59..000000000 --- a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/inner_wall.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:b9aeeeb1e1e6905d094524c86961e0be1a18bb19f705ba49a826a75656642fe6 -size 389434 diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/outer_wall.vtp b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/outer_wall.vtp deleted file mode 100644 index 6ec31850c..000000000 --- a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/outer_wall.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:f3ca33b843158a517a47ef3d8abc38a86b9356993beff11d696a6d24f19f8e9e -size 388018 diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/top.vtp b/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/top.vtp deleted file mode 100644 index 567d3b5c5..000000000 --- a/tests/cases/fluid/darcy_convergence/fine-mesh/mesh-surfaces/top.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:e0bf4150e7937e82039433eceb8a772241077639f2446292e1e6143528a9bf26 -size 381992 diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/svFSI.xml b/tests/cases/fluid/darcy_convergence/fine-mesh/svFSI.xml deleted file mode 100644 index a5974a32d..000000000 --- a/tests/cases/fluid/darcy_convergence/fine-mesh/svFSI.xml +++ /dev/null @@ -1,92 +0,0 @@ - - - - - 0 - 3 - 20 - 0.001 - 0.50 - STOP_SIM - - 1 - result - 1 - 0 - - 5 - 1 - - 1 - 0 - 1 - - - - - tissue_vol.vtu - - - mesh-surfaces/inner_wall.vtp - - - - mesh-surfaces/outer_wall.vtp - - - - mesh-surfaces/bottom.vtp - - - - mesh-surfaces/top.vtp - - - - - - - false - 2 - 5 - 1e-6 - - 1 - 1 - 1e-11 - 0.0 - - - true - true - - - - - fsils - - 1e-6 - 100 - 50 - - - - Dirichlet - Steady - 1.0 - false - - - - Dirichlet - Steady - 0.0 - false - - - - - - - - diff --git a/tests/cases/fluid/darcy_convergence/fine-mesh/tissue_vol.vtu b/tests/cases/fluid/darcy_convergence/fine-mesh/tissue_vol.vtu deleted file mode 100644 index 8f03b4a00..000000000 --- a/tests/cases/fluid/darcy_convergence/fine-mesh/tissue_vol.vtu +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:08c08fc4a71a6170e5129d82d0939a94781d4f39a7d8987077af53eaa7fd6edb -size 9690906 diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/bottom.vtp b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/bottom.vtp deleted file mode 100644 index 3d1ff6632..000000000 --- a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/bottom.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:ca7129083bd1a7d2ba8e56fecb0b4e0c83e03c4b90dc71770a01135d000ad3e8 -size 61132 diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/inner_wall.vtp b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/inner_wall.vtp deleted file mode 100644 index 57cb3ac98..000000000 --- a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/inner_wall.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:29535ba2be5b2ce0956eb9847e8f306ecd7bb34a8cb856c71ee22afedf924b6d -size 65378 diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/outer_wall.vtp b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/outer_wall.vtp deleted file mode 100644 index 38abe1060..000000000 --- a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/outer_wall.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:69bfbca17619c9288260303cfa6639c35163bd0e73897d78d7398611152dfbb9 -size 65194 diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/top.vtp b/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/top.vtp deleted file mode 100644 index 2da7351b5..000000000 --- a/tests/cases/fluid/darcy_convergence/med-mesh/mesh-surfaces/top.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:e7e31b71e8e44475ad02880a745b82aee444868e031fd29aa939f6962f2a468a -size 61108 diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/svFSI.xml b/tests/cases/fluid/darcy_convergence/med-mesh/svFSI.xml deleted file mode 100644 index a5974a32d..000000000 --- a/tests/cases/fluid/darcy_convergence/med-mesh/svFSI.xml +++ /dev/null @@ -1,92 +0,0 @@ - - - - - 0 - 3 - 20 - 0.001 - 0.50 - STOP_SIM - - 1 - result - 1 - 0 - - 5 - 1 - - 1 - 0 - 1 - - - - - tissue_vol.vtu - - - mesh-surfaces/inner_wall.vtp - - - - mesh-surfaces/outer_wall.vtp - - - - mesh-surfaces/bottom.vtp - - - - mesh-surfaces/top.vtp - - - - - - - false - 2 - 5 - 1e-6 - - 1 - 1 - 1e-11 - 0.0 - - - true - true - - - - - fsils - - 1e-6 - 100 - 50 - - - - Dirichlet - Steady - 1.0 - false - - - - Dirichlet - Steady - 0.0 - false - - - - - - - - diff --git a/tests/cases/fluid/darcy_convergence/med-mesh/tissue_vol.vtu b/tests/cases/fluid/darcy_convergence/med-mesh/tissue_vol.vtu deleted file mode 100644 index 98ab74186..000000000 --- a/tests/cases/fluid/darcy_convergence/med-mesh/tissue_vol.vtu +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:de64d090aed9534ff5bd50a438e25d7e90d2d1f89558568996c8378f0cb8c3bf -size 1032350 diff --git a/tests/cases/fluid/darcy_cylinder/darcy_validation_plot.png b/tests/cases/fluid/darcy_cylinder/darcy_validation_plot.png deleted file mode 100644 index b5fd098f6..000000000 --- a/tests/cases/fluid/darcy_cylinder/darcy_validation_plot.png +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:55c9f391196e0a0ba7d5d09cdade95fc6531e0ad741f4581dab9c2d93c598bee -size 148969 diff --git a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/bottom.vtp b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/bottom.vtp deleted file mode 100644 index c98f09ae0..000000000 --- a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/bottom.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:8c47e55f124d935488fb0a6da47cf417e5225fb6ede73b2edbcd1239a2b4f498 -size 83044 diff --git a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/inner_wall.vtp b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/inner_wall.vtp deleted file mode 100644 index c558c1ee7..000000000 --- a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/inner_wall.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:78b09a14424cd64e35a4be1d67831bfa52ce92d444737649b99eadfefd188fdd -size 92490 diff --git a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/outer_wall.vtp b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/outer_wall.vtp deleted file mode 100644 index 0725d76a2..000000000 --- a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/outer_wall.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:b749fd60aa89473e6d31efa576170bc94fb54aabccdabf71ca19b3540c15e189 -size 92130 diff --git a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/top.vtp b/tests/cases/fluid/darcy_cylinder/mesh-surfaces/top.vtp deleted file mode 100644 index 5258f8abe..000000000 --- a/tests/cases/fluid/darcy_cylinder/mesh-surfaces/top.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:eb329a65a000898e8a6cc49d306cf34dade035e9a37fa7c0f42de683ac571dac -size 83072 diff --git a/tests/cases/fluid/darcy_cylinder/result_020.vtu b/tests/cases/fluid/darcy_cylinder/result_020.vtu deleted file mode 100644 index 64bb9c721..000000000 --- a/tests/cases/fluid/darcy_cylinder/result_020.vtu +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:6c787705b921939503b9968d3139ffcad41e2d17d94254be9f2487a5e7f830e6 -size 1745017 diff --git a/tests/cases/fluid/darcy_cylinder/svFSI.xml b/tests/cases/fluid/darcy_cylinder/svFSI.xml deleted file mode 100644 index a5974a32d..000000000 --- a/tests/cases/fluid/darcy_cylinder/svFSI.xml +++ /dev/null @@ -1,92 +0,0 @@ - - - - - 0 - 3 - 20 - 0.001 - 0.50 - STOP_SIM - - 1 - result - 1 - 0 - - 5 - 1 - - 1 - 0 - 1 - - - - - tissue_vol.vtu - - - mesh-surfaces/inner_wall.vtp - - - - mesh-surfaces/outer_wall.vtp - - - - mesh-surfaces/bottom.vtp - - - - mesh-surfaces/top.vtp - - - - - - - false - 2 - 5 - 1e-6 - - 1 - 1 - 1e-11 - 0.0 - - - true - true - - - - - fsils - - 1e-6 - 100 - 50 - - - - Dirichlet - Steady - 1.0 - false - - - - Dirichlet - Steady - 0.0 - false - - - - - - - - diff --git a/tests/cases/fluid/darcy_cylinder/tissue_vol.vtu b/tests/cases/fluid/darcy_cylinder/tissue_vol.vtu deleted file mode 100644 index 09862e9e4..000000000 --- a/tests/cases/fluid/darcy_cylinder/tissue_vol.vtu +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:518e93edb3a819af183d35aa3ccb35a7a41d10e2546b3abeb869464d5015946e -size 1448656 From 993bbc5f27b90570cc187ff3ac686d91867467eb Mon Sep 17 00:00:00 2001 From: Michael Date: Tue, 1 Sep 2026 10:51:32 -0700 Subject: [PATCH 07/28] Addressing first set of comments from Michele (includes changing 'PhysicalProperyType' in other physics modules) --- Code/Source/solver/CMakeLists.txt | 2 +- Code/Source/solver/ComMod.h | 2 +- Code/Source/solver/Parameters.cpp | 4 +- Code/Source/solver/Parameters.h | 4 +- Code/Source/solver/cmm.cpp | 24 +-- Code/Source/solver/consts.h | 8 +- Code/Source/solver/darcy.cpp | 91 +++++----- Code/Source/solver/fluid.cpp | 36 ++-- Code/Source/solver/heatf.cpp | 8 +- Code/Source/solver/heats.cpp | 12 +- Code/Source/solver/l_elas.cpp | 22 +-- Code/Source/solver/mat_models.cpp | 12 +- Code/Source/solver/post.cpp | 30 ++-- Code/Source/solver/read_files.cpp | 46 +++--- Code/Source/solver/read_files.h | 2 +- Code/Source/solver/set_bc.cpp | 4 +- Code/Source/solver/set_equation_props.h | 210 ++++++++++++------------ Code/Source/solver/set_output_props.h | 2 +- Code/Source/solver/shells.cpp | 30 ++-- Code/Source/solver/stokes.cpp | 24 +-- Code/Source/solver/sv_struct.cpp | 16 +- Code/Source/solver/txt.cpp | 2 +- Code/Source/solver/ustruct.cpp | 22 +-- Code/Source/solver/vtk_xml.cpp | 2 +- 24 files changed, 314 insertions(+), 301 deletions(-) diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 97a0b3660..266031484 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -170,7 +170,6 @@ set(CSRCS SimulationLogger.h VtkData.h VtkData.cpp - active_stress_regazzoni.h active_stress_regazzoni.cpp all_fun.h all_fun.cpp baf_ini.h baf_ini.cpp bf.h bf.cpp @@ -236,6 +235,7 @@ set(CSRCS active_stress_uniform_unsteady.cpp active_stress_ode.cpp active_stress_nash_panfilov.cpp + active_stress_regazzoni.cpp SPLIT.c diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 2e9f60068..3d2e389be 100644 --- a/Code/Source/solver/ComMod.h +++ b/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 prop; + std::map prop; //double prop[consts::maxNProp]; // Electrophysiology model diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index fdcfadc23..c60291bbd 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2045,8 +2045,8 @@ DomainParameters::DomainParameters() { set_parameter("Penalty_parameter", 0.0, !required, penalty_parameter); set_parameter("Poisson_ratio", 0.3, !required, poisson_ratio); - set_parameter("Permeability", 0.0, !required, permeability); - set_parameter("Media_compressibility", 0.0, !required, media_compressibility); + set_parameter("Permeability", 1e-15, !required, darcy_permeability); + set_parameter("Media_compressibility", 0.0, !required, darcy_media_compressibility); set_parameter("Darcy_fluid_viscosity", 1.0, !required, darcy_fluid_viscosity); set_parameter("Relative_tolerance", 1e-4, !required, relative_tolerance); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 29445c859..f2423ec0e 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1623,8 +1623,8 @@ class DomainParameters : public ParameterLists Parameter source_term; Parameter time_step_for_integration; - Parameter permeability; - Parameter media_compressibility; + Parameter darcy_permeability; + Parameter darcy_media_compressibility; Parameter darcy_fluid_viscosity; // Inverse of Darcy permeability. Default value of 0.0 for Navier-Stokes and non-zero for Navier-Stokes-Brinkman Parameter inverse_darcy_permeability; diff --git a/Code/Source/solver/cmm.cpp b/Code/Source/solver/cmm.cpp index 3d78b5b96..2a62cb660 100644 --- a/Code/Source/solver/cmm.cpp +++ b/Code/Source/solver/cmm.cpp @@ -37,10 +37,10 @@ void cmm_3d(ComMod& com_mod, const int eNoN, const double w, const Vector 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 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; @@ -430,10 +430,10 @@ void cmm_mass(ComMod& com_mod, const double w, const Vector& N, const Ar #endif Vector 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 ; @@ -445,7 +445,7 @@ void cmm_mass(ComMod& com_mod, const double w, const Vector& 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; @@ -498,7 +498,7 @@ void cmm_stiffness(ComMod& com_mod, const Array& Nxi, const Array& Nxi, const Array output_type_name_to_type; /// @brief Possible physical properties. Current maxNPror is 20. // -enum class PhysicalProperyType +enum class PhysicalPropertyType { NA = 0, fluid_density = 1, @@ -419,8 +419,8 @@ enum class PhysicalProperyType ctau_M = 13, // stabilization coeffs. for USTRUCT (momentum, continuity) ctau_C = 14, inverse_darcy_permeability = 15, - permeability = 16, - media_compressibility = 17, + darcy_permeability = 16, + darcy_media_compressibility = 17, darcy_fluid_viscosity = 18 }; diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 01fbfec91..989e5537c 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -1,3 +1,11 @@ +#include "darcy.h" + +#include "all_fun.h" +#include "mat_fun.h" +#include "nn.h" +#include "utils.h" + +namespace darcy { /* This code implements the Darcy Equation for 2D and 3D problems in perfusion of porous media. @@ -9,32 +17,37 @@ - Assumptions of Stokes Flow - Steady-State ------------------------------------------------------------- - Strong form of the Single-Compartment Darcy equation: - u = -K∇(P) - ∇⋅u = β0(P_source - P) - β1(P - P_sink) - where: - u -> Volume flux vector - K -> Permeability tensor - P -> Pressure - ------------------------------------------------------------- - Weak form of the Single-Compartment Darcy equation: - -∫(∇q∇P)dΩ - λ∫qPdΩ = ∫qFdΩ - ∫q∇P⋅nvdΓ - where: - q -> Test function - λ -> (β0 + β1)/K - F -> -(β0(P_source) + β1(P_sink))/K - n -> Normal vector to the boundary + * Strong form of the Single-Compartment Darcy equation: + * \f[ u = -K\nabla P \f] + * \f[ \nabla \cdot u = \beta_0(P_{source} - P) - \beta_1(P - P_{sink}) \f] + * Note: See equations 8(a)/(b) in https://doi.org/10.1007/s10439-020-02681-z + * + * Combined Strong form: + * \f[ -\nabla \cdot (K \nabla P) - \beta_0(P_{source} - P) + \beta_1(P - P_{sink}) = 0 \f] + * + * Where: + * - \f$ u \f$ : Darcy flux + * - \f$ K \f$ : Permeability tensor + * - \f$ P \f$ : Pressure + * - \f$ P_{source} \f$ : Source pressure (e.g., arterial pressure) + * - \f$ P_{sink} \f$ : Sink pressure (e.g., venous/extraction pressure) + * - \f$ \beta_0 \f$ : Source coupling term (describes conductance of flow entering myocardium) + * - \f$ \beta_1 \f$ : Sink coupling term (describes conductance of flow exiting myocardium) + * + * ------------------------------------------------------------- + * + * Weak form of the Single-Compartment Darcy equation: + * \f[ -\int_{\Omega} (\nabla q \cdot \nabla P) d\Omega - \lambda \int_{\Omega} q P d\Omega = \int_{\Omega} q F d\Omega - \int_{\Gamma} q (\nabla P \cdot n) d\Gamma \f] + * + * Where: + * - \f$ q \f$ : Test function + * - \f$ \lambda \f$ : \f$ \frac{\beta_0 + \beta_1}{K} \f$ + * - \f$ F \f$ : \f$ -\frac{\beta_0 P_{source} + \beta_1 P_{sink}}{K} \f$ + * - \f$ n \f$ : Normal vector to the boundary + * - \f$ \Omega \f$ : Computational domain + * - \f$ \Gamma \f$ : Domain boundary */ -#include "darcy.h" - -#include "all_fun.h" -#include "mat_fun.h" -#include "nn.h" -#include "utils.h" - -namespace darcy { - void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR) { for (int a = 0; a < eNoN; a++) { @@ -151,11 +164,11 @@ void darcy_1d(ComMod& com_mod, const int eNoN, const double w, const Vector f(2); // f_x is internal force in x-direction; what is internal force? - f[0] = dmn.prop[PhysicalProperyType::f_x]; + f[0] = dmn.prop[PhysicalPropertyType::f_x]; - f[1] = dmn.prop[PhysicalProperyType::f_y]; + f[1] = dmn.prop[PhysicalPropertyType::f_y]; double T1 = eq.af * eq.gam * dt; double amd = eq.am / T1; @@ -1108,13 +1108,13 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e double ctM = 1.0; double ctC = 36.0; - double rho = dmn.prop[PhysicalProperyType::fluid_density]; + double rho = dmn.prop[PhysicalPropertyType::fluid_density]; Vector f(2); // f_x is internal force in x-direction; what is internal force? - f[0] = dmn.prop[PhysicalProperyType::f_x]; + f[0] = dmn.prop[PhysicalPropertyType::f_x]; - f[1] = dmn.prop[PhysicalProperyType::f_y]; + f[1] = dmn.prop[PhysicalPropertyType::f_y]; double T1 = eq.af * eq.gam * dt; double amd = eq.am / T1; @@ -1470,11 +1470,11 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e const double ctM = 1.0; const double ctC = 36.0; - double rho = dmn.prop[PhysicalProperyType::fluid_density]; + double rho = dmn.prop[PhysicalPropertyType::fluid_density]; double f[3]; - f[0] = dmn.prop[PhysicalProperyType::f_x]; - f[1] = dmn.prop[PhysicalProperyType::f_y]; - f[2] = dmn.prop[PhysicalProperyType::f_z]; + f[0] = dmn.prop[PhysicalPropertyType::f_x]; + f[1] = dmn.prop[PhysicalPropertyType::f_y]; + f[2] = dmn.prop[PhysicalPropertyType::f_z]; double T1 = eq.af * eq.gam * dt; double amd = eq.am / T1; @@ -1796,11 +1796,11 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e double ctM = 1.0; double ctC = 36.0; - double rho = dmn.prop[PhysicalProperyType::fluid_density]; + double rho = dmn.prop[PhysicalPropertyType::fluid_density]; std::array f; - f[0] = dmn.prop[PhysicalProperyType::f_x]; - f[1] = dmn.prop[PhysicalProperyType::f_y]; - f[2] = dmn.prop[PhysicalProperyType::f_z]; + f[0] = dmn.prop[PhysicalPropertyType::f_x]; + f[1] = dmn.prop[PhysicalPropertyType::f_y]; + f[2] = dmn.prop[PhysicalPropertyType::f_z]; double T1 = eq.af * eq.gam * dt; double amd = eq.am / T1; diff --git a/Code/Source/solver/heatf.cpp b/Code/Source/solver/heatf.cpp index e5be7f06f..10ab9e11d 100644 --- a/Code/Source/solver/heatf.cpp +++ b/Code/Source/solver/heatf.cpp @@ -158,8 +158,8 @@ void heatf_2d(ComMod& com_mod, const int eNoN, const double w, const Vector f({dmn.prop.at(PhysicalProperyType::f_x), - dmn.prop.at(PhysicalProperyType::f_y),}); + Vector f({dmn.prop.at(PhysicalPropertyType::f_x), + dmn.prop.at(PhysicalPropertyType::f_y),}); int i = eq.s; int j = i + 1; @@ -269,13 +269,13 @@ void l_elas_3d(ComMod& com_mod, const int eNoN, const double w, const Vector f({dmn.prop.at(PhysicalProperyType::f_x), - dmn.prop.at(PhysicalProperyType::f_y), - dmn.prop.at(PhysicalProperyType::f_z)}); + Vector f({dmn.prop.at(PhysicalPropertyType::f_x), + dmn.prop.at(PhysicalPropertyType::f_y), + dmn.prop.at(PhysicalPropertyType::f_z)}); int i = eq.s; int j = i + 1; diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index caef69cea..b23120176 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -1472,11 +1472,11 @@ void compute_tau(const ComMod& com_mod, const dmnType& lDmn, const double detF, using namespace consts; double he = 0.50 * pow(Je,1.0/static_cast(com_mod.nsd)); - double rho0 = lDmn.prop.at(PhysicalProperyType::solid_density); - double Em = lDmn.prop.at(PhysicalProperyType::elasticity_modulus); - double nu = lDmn.prop.at(PhysicalProperyType::poisson_ratio); - double ctM = lDmn.prop.at(PhysicalProperyType::ctau_M); - double ctC = lDmn.prop.at(PhysicalProperyType::ctau_C); + double rho0 = lDmn.prop.at(PhysicalPropertyType::solid_density); + double Em = lDmn.prop.at(PhysicalPropertyType::elasticity_modulus); + double nu = lDmn.prop.at(PhysicalPropertyType::poisson_ratio); + double ctM = lDmn.prop.at(PhysicalPropertyType::ctau_M); + double ctC = lDmn.prop.at(PhysicalPropertyType::ctau_C); double mu = 0.50*Em / (1.0 + nu); double c = 0.0; @@ -1513,7 +1513,7 @@ void g_vol_pen(const ComMod& com_mod, const dmnType& lDmn, const double p, { using namespace consts; - ro = lDmn.prop.at(PhysicalProperyType::solid_density) / Ja; + ro = lDmn.prop.at(PhysicalPropertyType::solid_density) / Ja; bt = 0.0; dbt = 0.0; dro = 0.0; diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 14d04eeac..344ea1b9f 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -853,7 +853,7 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S int nsd = com_mod.nsd; if ((outGrp == OutputNameType::outGrp_eFlx) && (com_mod.dmnId.size() == 0)) { - double rho = eq.dmn[0].prop[PhysicalProperyType::fluid_density]; + double rho = eq.dmn[0].prop[PhysicalPropertyType::fluid_density]; for (int a = 0; a < lM.nNo; a++) { int Ac = lM.gN(a); double p = lY(nsd,Ac); @@ -978,7 +978,7 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S // Energy flux calculation // } else if (outGrp == OutputNameType::outGrp_eFlx) { - double rho = eq.dmn[cDmn].prop[PhysicalProperyType::fluid_density]; + double rho = eq.dmn[cDmn].prop[PhysicalPropertyType::fluid_density]; double p = 0.0; Vector u(nsd); Vector lRes(maxNSD); @@ -998,7 +998,7 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S // Heat flux calculation // } else if (outGrp == OutputNameType::outGrp_hFlx) { - double kappa = eq.dmn[cDmn].prop[PhysicalProperyType::conductivity]; + double kappa = eq.dmn[cDmn].prop[PhysicalPropertyType::conductivity]; int i = eq.s; if (eq.phys == EquationType::phys_heatF) { @@ -1031,20 +1031,20 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S // MBF Flux calculation // - } else if (outGrp == OutputNameType::outGrp_mbfFlx) { + } else if (outGrp == OutputNameType::outGrp_MBFFlux) { const double permeability = - eq.dmn[cDmn].prop[PhysicalProperyType::permeability]; + eq.dmn[cDmn].prop[PhysicalPropertyType::darcy_permeability]; const double viscosity = - eq.dmn[cDmn].prop[PhysicalProperyType::darcy_fluid_viscosity]; + eq.dmn[cDmn].prop[PhysicalPropertyType::darcy_fluid_viscosity]; const double mobility = permeability / viscosity; - const int i = eq.s; + const int equation_index = eq.s; Vector 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); + dp_ds = dp_ds + Nx(0,a) * yl(equation_index,a); } for (int j = 0; j < nsd; j++) { grad_p(j) = dp_ds * fiber_tangent(j); @@ -1052,14 +1052,14 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S } 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); + grad_p(j) = grad_p(j) + Nx(j,a) * yl(equation_index,a); } } } - for (int j = 0; j < nsd; j++) { + for (int j = 0; j < nsd; j++) { lRes(j) = -mobility * grad_p(j); - } + } // Strain tensor invariants calculation // @@ -1310,8 +1310,8 @@ void shl_post(Simulation* simulation, const mshType& lM, const int m, Array double w = 0.0; if (cPhys == EquationType::phys_lElas) { - elM = eq.dmn[cDmn].prop[PhysicalProperyType::elasticity_modulus]; - nu = eq.dmn[cDmn].prop[PhysicalProperyType::poisson_ratio]; + elM = eq.dmn[cDmn].prop[PhysicalPropertyType::elasticity_modulus]; + nu = eq.dmn[cDmn].prop[PhysicalPropertyType::poisson_ratio]; lambda = elM*nu / (1.0 + nu) / (1.0 - 2.0*nu); mu = 0.5*elM / (1.0 + nu); } diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index ff069516a..6417a51c2 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1513,43 +1513,43 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& auto prop = propList[iProp][iPhys]; switch (prop) { - case PhysicalProperyType::backflow_stab: + case PhysicalPropertyType::backflow_stab: rtmp = domain_params->backflow_stabilization_coefficient.value(); break; - case PhysicalProperyType::conductivity: + case PhysicalPropertyType::conductivity: rtmp = domain_params->conductivity.value(); break; - case PhysicalProperyType::ctau_C: + case PhysicalPropertyType::ctau_C: rtmp = domain_params->continuity_stabilization_coefficient.value(); break; - case PhysicalProperyType::ctau_M: + case PhysicalPropertyType::ctau_M: rtmp = domain_params->momentum_stabilization_coefficient.value(); break; - case PhysicalProperyType::damping: + case PhysicalPropertyType::damping: rtmp = domain_params->mass_damping.value(); break; - case PhysicalProperyType::elasticity_modulus: + case PhysicalPropertyType::elasticity_modulus: rtmp = domain_params->elasticity_modulus.value(); break; - case PhysicalProperyType::f_x: + case PhysicalPropertyType::f_x: rtmp = domain_params->force_x.value(); break; - case PhysicalProperyType::f_y: + case PhysicalPropertyType::f_y: rtmp = domain_params->force_y.value(); break; - case PhysicalProperyType::f_z: + case PhysicalPropertyType::f_z: rtmp = domain_params->force_z.value(); break; - case PhysicalProperyType::fluid_density: + case PhysicalPropertyType::fluid_density: if (lEq.phys == EquationType::phys_CMM || lEq.phys == EquationType::phys_darcy) { rtmp = domain_params->fluid_density.value(); } else { @@ -1557,15 +1557,15 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& } break; - case PhysicalProperyType::poisson_ratio: + case PhysicalPropertyType::poisson_ratio: rtmp = domain_params->poisson_ratio.value(); break; - case PhysicalProperyType::shell_thickness: + case PhysicalPropertyType::shell_thickness: rtmp = domain_params->shell_thickness.value(); break; - case PhysicalProperyType::solid_density: + case PhysicalPropertyType::solid_density: if (lEq.phys == EquationType::phys_CMM) { rtmp = domain_params->solid_density.value(); } else { @@ -1573,23 +1573,23 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& } break; - case PhysicalProperyType::source_term: + case PhysicalPropertyType::source_term: rtmp = domain_params->source_term.value(); break; - case PhysicalProperyType::inverse_darcy_permeability: + case PhysicalPropertyType::inverse_darcy_permeability: rtmp = domain_params->inverse_darcy_permeability.value(); break; - case PhysicalProperyType::permeability: - rtmp = domain_params->permeability.value(); + case PhysicalPropertyType::darcy_permeability: + rtmp = domain_params->darcy_permeability.value(); break; - case PhysicalProperyType::media_compressibility: - rtmp = domain_params->media_compressibility.value(); + case PhysicalPropertyType::darcy_media_compressibility: + rtmp = domain_params->darcy_media_compressibility.value(); break; - case PhysicalProperyType::darcy_fluid_viscosity: + case PhysicalPropertyType::darcy_fluid_viscosity: rtmp = domain_params->darcy_fluid_viscosity.value(); break; } @@ -1734,7 +1734,7 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq) if (eq_params->use_taylor_hood_type_basis.defined()) { THflag = eq_params->use_taylor_hood_type_basis.value(); } - EquationProps propL{consts::PhysicalProperyType::NA}; + EquationProps propL{consts::PhysicalPropertyType::NA}; EquationOutputs outPuts; EquationNdop nDOP; @@ -2246,8 +2246,8 @@ void read_mat_model(Simulation* simulation, EquationParameters* eq_params, Domai using namespace consts; // Domain properties: elasticity modulus, poisson ratio - double E = lDmn.prop[PhysicalProperyType::elasticity_modulus]; - double nu = lDmn.prop[PhysicalProperyType::poisson_ratio]; + double E = lDmn.prop[PhysicalPropertyType::elasticity_modulus]; + double nu = lDmn.prop[PhysicalPropertyType::poisson_ratio]; // Shear modulus double mu = 0.5 * E / (1.0 + nu); diff --git a/Code/Source/solver/read_files.h b/Code/Source/solver/read_files.h index 197f83fdb..001e1a6d2 100644 --- a/Code/Source/solver/read_files.h +++ b/Code/Source/solver/read_files.h @@ -19,7 +19,7 @@ namespace read_files_ns { using EquationNdop = std::array; using EquationOutputs = std::array; using EquationPhys = std::vector; - using EquationProps = std::array, 20>; + using EquationProps = std::array, 20>; void face_match(ComMod& com_mod, faceType& lFa, faceType& gFa, Vector& ptr); diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 65f744727..a235bd286 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -1561,9 +1561,9 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const int iM = lFa.iM; int cDmn_local = all_fun::domain(com_mod, com_mod.msh[iM], cEq, lFa.gE(0)); double rho = eq.dmn[cDmn_local].prop.at( - consts::PhysicalProperyType::fluid_density); + consts::PhysicalPropertyType::fluid_density); double beta = eq.dmn[cDmn_local].prop.at( - consts::PhysicalProperyType::backflow_stab); + consts::PhysicalPropertyType::backflow_stab); double A = lFa.area; if (A > 0.0) { double u_n = Q_3D / A; // face-averaged normal velocity (< 0) diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index 27e8d003d..3862bbd8f 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -26,13 +26,13 @@ SetEquationPropertiesMapType set_equation_props = { auto& cep_mod = simulation->get_cep_mod(); lEq.phys = consts::EquationType::phys_CEP; - propL[0][0] = PhysicalProperyType::fluid_density; - propL[1][0] = PhysicalProperyType::backflow_stab; - propL[2][0] = PhysicalProperyType::f_x; - propL[3][0] = PhysicalProperyType::f_y; + propL[0][0] = PhysicalPropertyType::fluid_density; + propL[1][0] = PhysicalPropertyType::backflow_stab; + propL[2][0] = PhysicalPropertyType::f_x; + propL[3][0] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[4][0] = PhysicalProperyType::f_z; + propL[4][0] = PhysicalPropertyType::f_z; } cep_mod.cepEq = true; @@ -110,21 +110,21 @@ SetEquationPropertiesMapType set_equation_props = { } if (!com_mod.cmmInit) { - propL[0][0] = PhysicalProperyType::fluid_density; - propL[1][0] = PhysicalProperyType::backflow_stab; - propL[2][0] = PhysicalProperyType::solid_density; - propL[3][0] = PhysicalProperyType::poisson_ratio; - propL[4][0] = PhysicalProperyType::damping; + propL[0][0] = PhysicalPropertyType::fluid_density; + propL[1][0] = PhysicalPropertyType::backflow_stab; + propL[2][0] = PhysicalPropertyType::solid_density; + propL[3][0] = PhysicalPropertyType::poisson_ratio; + propL[4][0] = PhysicalPropertyType::damping; if (!com_mod.cmmVarWall) { - propL[5][0] = PhysicalProperyType::shell_thickness; - propL[6][0] = PhysicalProperyType::elasticity_modulus; + propL[5][0] = PhysicalPropertyType::shell_thickness; + propL[6][0] = PhysicalPropertyType::elasticity_modulus; } - propL[7][0] = PhysicalProperyType::f_x; - propL[8][0] = PhysicalProperyType::f_y; + propL[7][0] = PhysicalPropertyType::f_x; + propL[8][0] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[9][0] = PhysicalProperyType::f_z; + propL[9][0] = PhysicalPropertyType::f_z; } nDOP = {12, 4, 3, 0}; @@ -144,16 +144,16 @@ SetEquationPropertiesMapType set_equation_props = { }; } else { - propL[0][0] = PhysicalProperyType::poisson_ratio; + propL[0][0] = PhysicalPropertyType::poisson_ratio; if (!com_mod.cmmVarWall) { - propL[1][0] = PhysicalProperyType::shell_thickness; - propL[2][0] = PhysicalProperyType::elasticity_modulus; + propL[1][0] = PhysicalPropertyType::shell_thickness; + propL[2][0] = PhysicalPropertyType::elasticity_modulus; } - propL[7][0] = PhysicalProperyType::f_x; - propL[8][0] = PhysicalProperyType::f_y; + propL[7][0] = PhysicalPropertyType::f_x; + propL[8][0] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[9][0] = PhysicalProperyType::f_z; + propL[9][0] = PhysicalPropertyType::f_z; } if (pstEq) { @@ -169,7 +169,7 @@ SetEquationPropertiesMapType set_equation_props = { if (com_mod.cmmInit) { for (auto& domain : lEq.dmn) { - domain.prop[PhysicalProperyType::solid_density] = 0.0; + domain.prop[PhysicalPropertyType::solid_density] = 0.0; } } @@ -189,14 +189,14 @@ SetEquationPropertiesMapType set_equation_props = { auto& com_mod = simulation->get_com_mod(); lEq.phys = consts::EquationType::phys_fluid; - propL[0][0] = PhysicalProperyType::fluid_density; - propL[1][0] = PhysicalProperyType::backflow_stab; - propL[2][0] = PhysicalProperyType::inverse_darcy_permeability; - propL[3][0] = PhysicalProperyType::f_x; - propL[4][0] = PhysicalProperyType::f_y; + propL[0][0] = PhysicalPropertyType::fluid_density; + propL[1][0] = PhysicalPropertyType::backflow_stab; + propL[2][0] = PhysicalPropertyType::inverse_darcy_permeability; + propL[3][0] = PhysicalPropertyType::f_x; + propL[4][0] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[5][0] = PhysicalProperyType::f_z; + propL[5][0] = PhysicalPropertyType::f_z; } // Set fluid domain properties. @@ -234,8 +234,8 @@ SetEquationPropertiesMapType set_equation_props = { auto& com_mod = simulation->get_com_mod(); lEq.phys = consts::EquationType::phys_heatF; - propL[0][0] = PhysicalProperyType::conductivity; - propL[1][0] = PhysicalProperyType::source_term; + propL[0][0] = PhysicalPropertyType::conductivity; + propL[1][0] = PhysicalPropertyType::source_term; read_domain(simulation, eq_params, lEq, propL); @@ -260,9 +260,9 @@ SetEquationPropertiesMapType set_equation_props = { auto& com_mod = simulation->get_com_mod(); lEq.phys = consts::EquationType::phys_heatS; - propL[0][0] = PhysicalProperyType::conductivity; - propL[1][0] = PhysicalProperyType::source_term; - propL[2][0] = PhysicalProperyType::solid_density; + propL[0][0] = PhysicalPropertyType::conductivity; + propL[1][0] = PhysicalPropertyType::source_term; + propL[2][0] = PhysicalPropertyType::solid_density; read_domain(simulation, eq_params, lEq, propL); @@ -284,11 +284,11 @@ SetEquationPropertiesMapType set_equation_props = { auto& com_mod = simulation->get_com_mod(); lEq.phys = consts::EquationType::phys_darcy; - propL[0][0] = PhysicalProperyType::permeability; - propL[1][0] = PhysicalProperyType::source_term; - propL[2][0] = PhysicalProperyType::fluid_density; - propL[3][0] = PhysicalProperyType::media_compressibility; - propL[4][0] = PhysicalProperyType::darcy_fluid_viscosity; + propL[0][0] = PhysicalPropertyType::darcy_permeability; + propL[1][0] = PhysicalPropertyType::source_term; + propL[2][0] = PhysicalPropertyType::fluid_density; + propL[3][0] = PhysicalPropertyType::darcy_media_compressibility; + propL[4][0] = PhysicalPropertyType::darcy_fluid_viscosity; read_domain(simulation, eq_params, lEq, propL); @@ -317,48 +317,48 @@ SetEquationPropertiesMapType set_equation_props = { // Set fluid properties. int n = 0; - propL[0][n] = PhysicalProperyType::fluid_density; - propL[1][n] = PhysicalProperyType::backflow_stab; - propL[2][n] = PhysicalProperyType::f_x; - propL[3][n] = PhysicalProperyType::f_y; + propL[0][n] = PhysicalPropertyType::fluid_density; + propL[1][n] = PhysicalPropertyType::backflow_stab; + propL[2][n] = PhysicalPropertyType::f_x; + propL[3][n] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[4][n] = PhysicalProperyType::f_z; + propL[4][n] = PhysicalPropertyType::f_z; } // Set struct properties. n += 1; - propL[0][n] = PhysicalProperyType::solid_density; - propL[1][n] = PhysicalProperyType::elasticity_modulus; - propL[2][n] = PhysicalProperyType::poisson_ratio; - propL[3][n] = PhysicalProperyType::damping; - propL[4][n] = PhysicalProperyType::f_x; - propL[5][n] = PhysicalProperyType::f_y; + propL[0][n] = PhysicalPropertyType::solid_density; + propL[1][n] = PhysicalPropertyType::elasticity_modulus; + propL[2][n] = PhysicalPropertyType::poisson_ratio; + propL[3][n] = PhysicalPropertyType::damping; + propL[4][n] = PhysicalPropertyType::f_x; + propL[5][n] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[6][n] = PhysicalProperyType::f_z; + propL[6][n] = PhysicalPropertyType::f_z; } // Set ustruct properties. n += 1; - propL[0][n] = PhysicalProperyType::solid_density; - propL[1][n] = PhysicalProperyType::elasticity_modulus; - propL[2][n] = PhysicalProperyType::poisson_ratio; - propL[3][n] = PhysicalProperyType::ctau_M; - propL[4][n] = PhysicalProperyType::ctau_C; - propL[5][n] = PhysicalProperyType::f_x; - propL[6][n] = PhysicalProperyType::f_y; + propL[0][n] = PhysicalPropertyType::solid_density; + propL[1][n] = PhysicalPropertyType::elasticity_modulus; + propL[2][n] = PhysicalPropertyType::poisson_ratio; + propL[3][n] = PhysicalPropertyType::ctau_M; + propL[4][n] = PhysicalPropertyType::ctau_C; + propL[5][n] = PhysicalPropertyType::f_x; + propL[6][n] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[7][n] = PhysicalProperyType::f_z; + propL[7][n] = PhysicalPropertyType::f_z; } // Set lElas properties. n += 1; - propL[0][n] = PhysicalProperyType::solid_density; - propL[1][n] = PhysicalProperyType::elasticity_modulus; - propL[2][n] = PhysicalProperyType::poisson_ratio; - propL[3][n] = PhysicalProperyType::f_x; - propL[4][n] = PhysicalProperyType::f_y; + propL[0][n] = PhysicalPropertyType::solid_density; + propL[1][n] = PhysicalPropertyType::elasticity_modulus; + propL[2][n] = PhysicalPropertyType::poisson_ratio; + propL[3][n] = PhysicalPropertyType::f_x; + propL[4][n] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[5][n] = PhysicalProperyType::f_z; + propL[5][n] = PhysicalPropertyType::f_z; } // Set lEq properties. @@ -412,13 +412,13 @@ SetEquationPropertiesMapType set_equation_props = { auto& com_mod = simulation->get_com_mod(); lEq.phys = consts::EquationType::phys_lElas; - propL[0][0] = PhysicalProperyType::solid_density; - propL[1][0] = PhysicalProperyType::elasticity_modulus; - propL[2][0] = PhysicalProperyType::poisson_ratio; - propL[3][0] = PhysicalProperyType::f_x; - propL[4][0] = PhysicalProperyType::f_y; + propL[0][0] = PhysicalPropertyType::solid_density; + propL[1][0] = PhysicalPropertyType::elasticity_modulus; + propL[2][0] = PhysicalPropertyType::poisson_ratio; + propL[3][0] = PhysicalPropertyType::f_x; + propL[4][0] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[5][0] = PhysicalProperyType::f_z; + propL[5][0] = PhysicalPropertyType::f_z; } read_domain(simulation, eq_params, lEq, propL); @@ -451,20 +451,20 @@ SetEquationPropertiesMapType set_equation_props = { auto& com_mod = simulation->get_com_mod(); lEq.phys = consts::EquationType::phys_mesh; - propL[0][0] = PhysicalProperyType::solid_density; - propL[1][0] = PhysicalProperyType::elasticity_modulus; - propL[2][0] = PhysicalProperyType::poisson_ratio; - propL[3][0] = PhysicalProperyType::f_x; - propL[4][0] = PhysicalProperyType::f_y; + propL[0][0] = PhysicalPropertyType::solid_density; + propL[1][0] = PhysicalPropertyType::elasticity_modulus; + propL[2][0] = PhysicalPropertyType::poisson_ratio; + propL[3][0] = PhysicalPropertyType::f_x; + propL[4][0] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[5][0] = PhysicalProperyType::f_z; + propL[5][0] = PhysicalPropertyType::f_z; } read_domain(simulation, eq_params, lEq, propL); for (auto& domain : lEq.dmn) { - domain.prop[PhysicalProperyType::solid_density] = 0.0; - domain.prop[PhysicalProperyType::elasticity_modulus] = 1.0; + domain.prop[PhysicalPropertyType::solid_density] = 0.0; + domain.prop[PhysicalPropertyType::elasticity_modulus] = 1.0; } nDOP = {3, 1, 0, 0}; @@ -489,14 +489,14 @@ SetEquationPropertiesMapType set_equation_props = { lEq.phys = consts::EquationType::phys_shell; com_mod.shlEq = true; - propL[0][0] = PhysicalProperyType::solid_density; - propL[1][0] = PhysicalProperyType::damping; - propL[2][0] = PhysicalProperyType::elasticity_modulus; - propL[3][0] = PhysicalProperyType::poisson_ratio; - propL[4][0] = PhysicalProperyType::shell_thickness; - propL[5][0] = PhysicalProperyType::f_x; - propL[6][0] = PhysicalProperyType::f_y; - propL[7][0] = PhysicalProperyType::f_z; + propL[0][0] = PhysicalPropertyType::solid_density; + propL[1][0] = PhysicalPropertyType::damping; + propL[2][0] = PhysicalPropertyType::elasticity_modulus; + propL[3][0] = PhysicalPropertyType::poisson_ratio; + propL[4][0] = PhysicalPropertyType::shell_thickness; + propL[5][0] = PhysicalPropertyType::f_x; + propL[6][0] = PhysicalPropertyType::f_y; + propL[7][0] = PhysicalPropertyType::f_z; read_domain(simulation, eq_params, lEq, propL); @@ -529,11 +529,11 @@ SetEquationPropertiesMapType set_equation_props = { auto& com_mod = simulation->get_com_mod(); lEq.phys = consts::EquationType::phys_stokes; - propL[0][0] = PhysicalProperyType::ctau_M; - propL[1][0] = PhysicalProperyType::f_x; - propL[2][0] = PhysicalProperyType::f_y; + propL[0][0] = PhysicalPropertyType::ctau_M; + propL[1][0] = PhysicalPropertyType::f_x; + propL[2][0] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[3][0] = PhysicalProperyType::f_z; + propL[3][0] = PhysicalPropertyType::f_z; } read_domain(simulation, eq_params, lEq, propL); @@ -565,14 +565,14 @@ SetEquationPropertiesMapType set_equation_props = { auto& com_mod = simulation->get_com_mod(); lEq.phys = consts::EquationType::phys_struct; - propL[0][0] = PhysicalProperyType::solid_density; - propL[1][0] = PhysicalProperyType::damping; - propL[2][0] = PhysicalProperyType::elasticity_modulus; - propL[3][0] = PhysicalProperyType::poisson_ratio; - propL[4][0] = PhysicalProperyType::f_x; - propL[5][0] = PhysicalProperyType::f_y; + propL[0][0] = PhysicalPropertyType::solid_density; + propL[1][0] = PhysicalPropertyType::damping; + propL[2][0] = PhysicalPropertyType::elasticity_modulus; + propL[3][0] = PhysicalPropertyType::poisson_ratio; + propL[4][0] = PhysicalPropertyType::f_x; + propL[5][0] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[6][0] = PhysicalProperyType::f_z; + propL[6][0] = PhysicalPropertyType::f_z; } read_domain(simulation, eq_params, lEq, propL); @@ -620,15 +620,15 @@ SetEquationPropertiesMapType set_equation_props = { lEq.phys = consts::EquationType::phys_ustruct; com_mod.sstEq = true; - propL[0][0] = PhysicalProperyType::solid_density; - propL[1][0] = PhysicalProperyType::elasticity_modulus; - propL[2][0] = PhysicalProperyType::poisson_ratio; - propL[3][0] = PhysicalProperyType::ctau_M; - propL[4][0] = PhysicalProperyType::ctau_C; - propL[5][0] = PhysicalProperyType::f_x; - propL[6][0] = PhysicalProperyType::f_y; + propL[0][0] = PhysicalPropertyType::solid_density; + propL[1][0] = PhysicalPropertyType::elasticity_modulus; + propL[2][0] = PhysicalPropertyType::poisson_ratio; + propL[3][0] = PhysicalPropertyType::ctau_M; + propL[4][0] = PhysicalPropertyType::ctau_C; + propL[5][0] = PhysicalPropertyType::f_x; + propL[6][0] = PhysicalPropertyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[7][0] = PhysicalProperyType::f_z; + propL[7][0] = PhysicalPropertyType::f_z; } read_domain(simulation, eq_params, lEq, propL); diff --git a/Code/Source/solver/set_output_props.h b/Code/Source/solver/set_output_props.h index 493d8684f..384441912 100644 --- a/Code/Source/solver/set_output_props.h +++ b/Code/Source/solver/set_output_props.h @@ -61,5 +61,5 @@ std::map output_props_map = {OutputNameType::out_vorticity, std::make_tuple(OutputNameType::outGrp_vort, 0, maxNSD, "Vorticity") }, {OutputNameType::out_WSS, std::make_tuple(OutputNameType::outGrp_WSS, 0, maxNSD, "WSS") }, {OutputNameType::out_MBF, std::make_tuple(OutputNameType::outGrp_Y, 0, 1, "MBF")}, - {OutputNameType::out_mbfFlux, std::make_tuple(OutputNameType::outGrp_mbfFlx, 0, nsd, "MBF_flux")} + {OutputNameType::out_mbfFlux, std::make_tuple(OutputNameType::outGrp_MBFFlux, 0, nsd, "MBF_flux")} }; diff --git a/Code/Source/solver/shells.cpp b/Code/Source/solver/shells.cpp index f17064463..96932a2df 100644 --- a/Code/Source/solver/shells.cpp +++ b/Code/Source/solver/shells.cpp @@ -188,11 +188,11 @@ void shell_3d(ComMod& com_mod, const mshType& lM, const int g, const int eNoN, const double dt = com_mod.dt; // Define parameters - double rho = eq.dmn[cDmn].prop.at(PhysicalProperyType::solid_density); - double dmp = dmn.prop.at(PhysicalProperyType::damping); - double ht = eq.dmn[cDmn].prop.at(PhysicalProperyType::shell_thickness); - Vector fb({dmn.prop.at(PhysicalProperyType::f_x), dmn.prop.at(PhysicalProperyType::f_y), - dmn.prop.at(PhysicalProperyType::f_z)}); + double rho = eq.dmn[cDmn].prop.at(PhysicalPropertyType::solid_density); + double dmp = dmn.prop.at(PhysicalPropertyType::damping); + double ht = eq.dmn[cDmn].prop.at(PhysicalPropertyType::shell_thickness); + Vector fb({dmn.prop.at(PhysicalPropertyType::f_x), dmn.prop.at(PhysicalPropertyType::f_y), + dmn.prop.at(PhysicalPropertyType::f_z)}); double amd = eq.am * rho + eq.af * eq.gam * dt * dmp; double afl = eq.af * eq.beta * dt * dt; @@ -561,8 +561,8 @@ void shell_bend_cst(ComMod& com_mod, const mshType& lM, const int e, const Vecto auto& dmn = eq.dmn[cDmn]; // Define parameters - double rho = eq.dmn[cDmn].prop.at(PhysicalProperyType::solid_density); - double dmp = dmn.prop.at(PhysicalProperyType::damping); + double rho = eq.dmn[cDmn].prop.at(PhysicalPropertyType::solid_density); + double dmp = dmn.prop.at(PhysicalPropertyType::damping); int nsd = com_mod.nsd; int eNoN = 2 * lM.eNoN; @@ -1206,12 +1206,12 @@ void shell_cst(ComMod& com_mod, const mshType& lM, const int e, const int eNoN, #endif // Define parameters - double rho = eq.dmn[cDmn].prop.at(PhysicalProperyType::solid_density); - double dmp = dmn.prop.at(PhysicalProperyType::damping); - double ht = eq.dmn[cDmn].prop.at(PhysicalProperyType::shell_thickness); - Vector fb({dmn.prop.at(PhysicalProperyType::f_x), - dmn.prop.at(PhysicalProperyType::f_y), - dmn.prop.at(PhysicalProperyType::f_z)}); + double rho = eq.dmn[cDmn].prop.at(PhysicalPropertyType::solid_density); + double dmp = dmn.prop.at(PhysicalPropertyType::damping); + double ht = eq.dmn[cDmn].prop.at(PhysicalPropertyType::shell_thickness); + Vector fb({dmn.prop.at(PhysicalPropertyType::f_x), + dmn.prop.at(PhysicalPropertyType::f_y), + dmn.prop.at(PhysicalPropertyType::f_z)}); double amd = eq.am * rho + eq.af * eq.gam * dt * dmp; double afl = eq.af * eq.beta * dt * dt; @@ -1721,8 +1721,8 @@ void shl_strs_res(const ComMod& com_mod, const dmnType& lDmn, const int nFn, con #endif // Set shell thickness - double ht = lDmn.prop.at(PhysicalProperyType::shell_thickness); - double nu = lDmn.prop.at(PhysicalProperyType::poisson_ratio); + double ht = lDmn.prop.at(PhysicalPropertyType::shell_thickness); + double nu = lDmn.prop.at(PhysicalPropertyType::poisson_ratio); // Check for incompressibility bool flag = false; diff --git a/Code/Source/solver/stokes.cpp b/Code/Source/solver/stokes.cpp index c0f1525a8..f888f6061 100644 --- a/Code/Source/solver/stokes.cpp +++ b/Code/Source/solver/stokes.cpp @@ -231,11 +231,11 @@ void stokes_2d_c(ComMod& com_mod, const int lStab, const int eNoNw, const int eN auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; double mu = dmn.fluid_visc.mu_i; - double ctM = dmn.prop[PhysicalProperyType::ctau_M]; + double ctM = dmn.prop[PhysicalPropertyType::ctau_M]; Vector fb(2); - fb[0] = dmn.prop[PhysicalProperyType::f_x]; - fb[1] = dmn.prop[PhysicalProperyType::f_y]; + fb[0] = dmn.prop[PhysicalPropertyType::f_x]; + fb[1] = dmn.prop[PhysicalPropertyType::f_y]; double wm = w * eq.am; double wf = w * eq.af * eq.gam * dt; @@ -355,8 +355,8 @@ void stokes_2d_m(ComMod& com_mod, const int eNoNw, const int eNoNq, const double double mu = dmn.fluid_visc.mu_i; Vector fb(2); - fb[0] = dmn.prop[PhysicalProperyType::f_x]; - fb[1] = dmn.prop[PhysicalProperyType::f_y]; + fb[0] = dmn.prop[PhysicalPropertyType::f_x]; + fb[1] = dmn.prop[PhysicalPropertyType::f_y]; double af = eq.af * eq.gam * dt; double wf = w * af; @@ -464,12 +464,12 @@ void stokes_3d_c(ComMod& com_mod, const int lStab, const int eNoNw, const int eN auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; double mu = dmn.fluid_visc.mu_i; - double ctM = dmn.prop[PhysicalProperyType::ctau_M]; + double ctM = dmn.prop[PhysicalPropertyType::ctau_M]; Vector fb(3); - fb[0] = dmn.prop[PhysicalProperyType::f_x]; - fb[1] = dmn.prop[PhysicalProperyType::f_y]; - fb[2] = dmn.prop[PhysicalProperyType::f_z]; + fb[0] = dmn.prop[PhysicalPropertyType::f_x]; + fb[1] = dmn.prop[PhysicalPropertyType::f_y]; + fb[2] = dmn.prop[PhysicalPropertyType::f_z]; double wm = w * eq.am; double wf = w * eq.af * eq.gam * dt; @@ -581,9 +581,9 @@ void stokes_3d_m(ComMod& com_mod, const int eNoNw, const int eNoNq, const double double mu = dmn.fluid_visc.mu_i; Vector fb(3); - fb[0] = dmn.prop[PhysicalProperyType::f_x]; - fb[1] = dmn.prop[PhysicalProperyType::f_y]; - fb[2] = dmn.prop[PhysicalProperyType::f_z]; + fb[0] = dmn.prop[PhysicalPropertyType::f_x]; + fb[1] = dmn.prop[PhysicalPropertyType::f_y]; + fb[2] = dmn.prop[PhysicalPropertyType::f_z]; double af = eq.af * eq.gam * dt; double wf = w * af; diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index 561ee73b5..f350c8c6c 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -368,9 +368,9 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, // Set parameters // - double rho = dmn.prop.at(PhysicalProperyType::solid_density); - double dmp = dmn.prop.at(PhysicalProperyType::damping); - Vector fb({dmn.prop.at(PhysicalProperyType::f_x), dmn.prop.at(PhysicalProperyType::f_y)}); + double rho = dmn.prop.at(PhysicalPropertyType::solid_density); + double dmp = dmn.prop.at(PhysicalPropertyType::damping); + Vector fb({dmn.prop.at(PhysicalPropertyType::f_x), dmn.prop.at(PhysicalPropertyType::f_y)}); double afu = eq.af * eq.beta*dt*dt; double afv = eq.af * eq.gam*dt; double amd = eq.am * rho + eq.af * eq.gam * dt * dmp; @@ -566,11 +566,11 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, // Set parameters // - double rho = dmn.prop.at(PhysicalProperyType::solid_density); - double dmp = dmn.prop.at(PhysicalProperyType::damping); - Vector fb({dmn.prop.at(PhysicalProperyType::f_x), - dmn.prop.at(PhysicalProperyType::f_y), - dmn.prop.at(PhysicalProperyType::f_z)}); + double rho = dmn.prop.at(PhysicalPropertyType::solid_density); + double dmp = dmn.prop.at(PhysicalPropertyType::damping); + Vector fb({dmn.prop.at(PhysicalPropertyType::f_x), + dmn.prop.at(PhysicalPropertyType::f_y), + dmn.prop.at(PhysicalPropertyType::f_z)}); double afu = eq.af * eq.beta*dt*dt; double afv = eq.af * eq.gam*dt; diff --git a/Code/Source/solver/txt.cpp b/Code/Source/solver/txt.cpp index 287e56c6d..e84ff556f 100644 --- a/Code/Source/solver/txt.cpp +++ b/Code/Source/solver/txt.cpp @@ -321,7 +321,7 @@ void txt(Simulation* simulation, const bool init_write, const SolutionStates& so case OutputNameType::outGrp_divV: case OutputNameType::outGrp_J: case OutputNameType::outGrp_mises: - case OutputNameType::outGrp_mbfFlx: + case OutputNameType::outGrp_MBFFlux: post::all_post(simulation, tmpV, solutions, oGrp, iEq); break; diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index a698050e7..87f5efa4d 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -449,9 +449,9 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in const double dt = com_mod.dt; Vector fb(2); - fb[0] = dmn.prop[PhysicalProperyType::f_x]; - fb[1] = dmn.prop[PhysicalProperyType::f_y]; - fb[2] = dmn.prop[PhysicalProperyType::f_z]; + fb[0] = dmn.prop[PhysicalPropertyType::f_x]; + fb[1] = dmn.prop[PhysicalPropertyType::f_y]; + fb[2] = dmn.prop[PhysicalPropertyType::f_z]; double am = eq.am; double af = eq.af * eq.gam * dt; @@ -653,9 +653,9 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in const double dt = com_mod.dt; Vector fb(3); - fb[0] = dmn.prop[PhysicalProperyType::f_x]; - fb[1] = dmn.prop[PhysicalProperyType::f_y]; - fb[2] = dmn.prop[PhysicalProperyType::f_z]; + fb[0] = dmn.prop[PhysicalPropertyType::f_x]; + fb[1] = dmn.prop[PhysicalPropertyType::f_y]; + fb[2] = dmn.prop[PhysicalPropertyType::f_z]; double am = eq.am; double af = eq.af * eq.gam * dt; @@ -902,8 +902,8 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // Define parameters // Vector fb(2); - fb[0] = dmn.prop[PhysicalProperyType::f_x]; - fb[1] = dmn.prop[PhysicalProperyType::f_y]; + fb[0] = dmn.prop[PhysicalPropertyType::f_x]; + fb[1] = dmn.prop[PhysicalPropertyType::f_y]; double am = eq.am; double af = eq.af * eq.gam * dt; @@ -1192,9 +1192,9 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // Define parameters Vector fb(3); - fb[0] = dmn.prop[PhysicalProperyType::f_x]; - fb[1] = dmn.prop[PhysicalProperyType::f_y]; - fb[2] = dmn.prop[PhysicalProperyType::f_z]; + fb[0] = dmn.prop[PhysicalPropertyType::f_x]; + fb[1] = dmn.prop[PhysicalPropertyType::f_y]; + fb[2] = dmn.prop[PhysicalPropertyType::f_z]; double am = eq.am; double af = eq.af * eq.gam * dt; diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 915c71565..5247dc871 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1164,7 +1164,7 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b case OutputNameType::outGrp_stInv: case OutputNameType::outGrp_vortex: case OutputNameType::outGrp_Visc: - case OutputNameType::outGrp_mbfFlx: + case OutputNameType::outGrp_MBFFlux: post::post(simulation, msh, tmpV, solutions, oGrp, iEq); for (int a = 0; a < msh.nNo; a++) { int Ac = msh.gN(a); From d0633979ec6b66f78d66bc3f2a517a27baa38be8 Mon Sep 17 00:00:00 2001 From: Michael Date: Tue, 1 Sep 2026 11:17:50 -0700 Subject: [PATCH 08/28] fixing unit test PhysicalPropertyType dependency --- tests/unitTests/material_model_tests/test_material_common.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unitTests/material_model_tests/test_material_common.h b/tests/unitTests/material_model_tests/test_material_common.h index 42f74498e..7aaa81104 100644 --- a/tests/unitTests/material_model_tests/test_material_common.h +++ b/tests/unitTests/material_model_tests/test_material_common.h @@ -380,7 +380,7 @@ class TestMaterialModel : public TestBase { */ void g_vol_pen(const double p, const double rho0, double &rho, double &beta, double &drho, double &dbeta, const double Ja) { auto &dmn = com_mod.mockEq.mockDmn; - dmn.prop[consts::PhysicalProperyType::solid_density] = rho0; // Set initial solid density + dmn.prop[consts::PhysicalPropertyType::solid_density] = rho0; // Set initial solid density mat_models::g_vol_pen(com_mod, dmn, p, rho, beta, drho, dbeta, Ja); } From 5682189b02579ad08f93cc10715fd759524a9256 Mon Sep 17 00:00:00 2001 From: Michael Date: Tue, 1 Sep 2026 16:27:51 -0700 Subject: [PATCH 09/28] Reducing darcy coverage for 2D/3D and changing output names --- Code/Source/solver/consts.h | 6 ++-- Code/Source/solver/darcy.cpp | 47 ++----------------------- Code/Source/solver/darcy.h | 3 -- Code/Source/solver/post.cpp | 40 ++++----------------- Code/Source/solver/set_equation_props.h | 2 +- Code/Source/solver/set_output_props.h | 4 +-- Code/Source/solver/txt.cpp | 2 +- Code/Source/solver/vtk_xml.cpp | 2 +- 8 files changed, 16 insertions(+), 90 deletions(-) diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index 49fbc0ff9..636daa9d4 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -349,7 +349,7 @@ enum class OutputNameType { outGrp_activeTensionFibers = 529, outGrp_activeTensionSheets = 530, outGrp_activeTensionNormal = 531, - outGrp_MBFFlux = 532, + outGrp_darcyFlux = 532, out_velocity = 599, out_pressure = 598, @@ -384,8 +384,8 @@ enum class OutputNameType { out_activeTensionFibers = 569, out_activeTensionSheets = 568, out_activeTensionNormal = 567, - out_MBF = 566, - out_mbfFlux = 565 + out_darcyPressure = 566, + out_darcyFlux = 565 }; /// @brief Simulation output file types. diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 989e5537c..2655640f1 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -29,7 +29,7 @@ namespace darcy { * - \f$ u \f$ : Darcy flux * - \f$ K \f$ : Permeability tensor * - \f$ P \f$ : Pressure - * - \f$ P_{source} \f$ : Source pressure (e.g., arterial pressure) + * - \f$ P_{source} \f$ : Source pressure (e.g., arterial pressure) * - \f$ P_{sink} \f$ : Sink pressure (e.g., venous/extraction pressure) * - \f$ \beta_0 \f$ : Source coupling term (describes conductance of flow entering myocardium) * - \f$ \beta_1 \f$ : Sink coupling term (describes conductance of flow exiting myocardium) @@ -141,10 +141,8 @@ void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& s darcy_3d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); } else if (insd == 2) { darcy_2d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); - } else if (insd == 1) { - darcy_1d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); } else { - throw std::runtime_error("[construct_darcy] insd must be 1, 2 or 3."); + throw std::runtime_error("[construct_darcy] insd must be 2 or 3."); } } @@ -152,47 +150,6 @@ void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& s } } -void darcy_1d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, - const Array& al, const Array& yl, Array& lR, Array3& 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(PhysicalPropertyType::darcy_permeability); - double source = dmn.prop.at(PhysicalPropertyType::source_term); - double beta_0 = dmn.prop.at(PhysicalPropertyType::darcy_media_compressibility); - double rho_0 = dmn.prop.at(PhysicalPropertyType::fluid_density); - double mu = dmn.prop.at(PhysicalPropertyType::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))))); - } - } -} - void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, const Array& al, const Array& yl, Array& lR, Array3& lK) { diff --git a/Code/Source/solver/darcy.h b/Code/Source/solver/darcy.h index b80535cff..d54073571 100644 --- a/Code/Source/solver/darcy.h +++ b/Code/Source/solver/darcy.h @@ -15,9 +15,6 @@ namespace darcy { void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& solutions); - void darcy_1d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, - const Array& al, const Array& yl, Array& lR, Array3& lK); - void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, const Array& al, const Array& yl, Array& lR, Array3& lK); diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 344ea1b9f..fb403b050 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -916,29 +916,11 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S Array ksix(nsd,nsd); double Jac = 0.0; - Vector 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; - } } double w = lM.w(g) * Jac; @@ -1029,9 +1011,9 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S } } - // MBF Flux calculation + // Darcy Flux calculation // - } else if (outGrp == OutputNameType::outGrp_MBFFlux) { + } else if (outGrp == OutputNameType::outGrp_darcyFlux) { const double permeability = eq.dmn[cDmn].prop[PhysicalPropertyType::darcy_permeability]; const double viscosity = @@ -1041,22 +1023,12 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S Vector 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(equation_index,a); - } + for (int a = 0; a < eNoN; a++) { for (int j = 0; j < nsd; j++) { - grad_p(j) = dp_ds * fiber_tangent(j); + grad_p(j) = grad_p(j) + Nx(j,a) * yl(equation_index,a); } - } 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(equation_index,a); - } - } - } - + } + for (int j = 0; j < nsd; j++) { lRes(j) = -mobility * grad_p(j); } diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index 3862bbd8f..f68a5569b 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -293,7 +293,7 @@ SetEquationPropertiesMapType set_equation_props = { read_domain(simulation, eq_params, lEq, propL); nDOP = {2,1,1,0}; - outPuts = {OutputNameType::out_MBF, OutputNameType::out_mbfFlux}; + outPuts = {OutputNameType::out_darcyPressure, OutputNameType::out_darcyFlux}; // Set solver parameters. read_ls(simulation, eq_params, SolverType::lSolver_CG, lEq); diff --git a/Code/Source/solver/set_output_props.h b/Code/Source/solver/set_output_props.h index 384441912..b5f6a8b41 100644 --- a/Code/Source/solver/set_output_props.h +++ b/Code/Source/solver/set_output_props.h @@ -60,6 +60,6 @@ std::map output_props_map = {OutputNameType::out_vortex, std::make_tuple(OutputNameType::outGrp_vortex, 0, 1, "Vortex") }, {OutputNameType::out_vorticity, std::make_tuple(OutputNameType::outGrp_vort, 0, maxNSD, "Vorticity") }, {OutputNameType::out_WSS, std::make_tuple(OutputNameType::outGrp_WSS, 0, maxNSD, "WSS") }, - {OutputNameType::out_MBF, std::make_tuple(OutputNameType::outGrp_Y, 0, 1, "MBF")}, - {OutputNameType::out_mbfFlux, std::make_tuple(OutputNameType::outGrp_MBFFlux, 0, nsd, "MBF_flux")} + {OutputNameType::out_darcyPressure, std::make_tuple(OutputNameType::outGrp_Y, 0, 1, "Darcy_pressure")}, + {OutputNameType::out_darcyFlux, std::make_tuple(OutputNameType::outGrp_darcyFlux, 0, nsd, "Darcy_flux")} }; diff --git a/Code/Source/solver/txt.cpp b/Code/Source/solver/txt.cpp index e84ff556f..009913bf7 100644 --- a/Code/Source/solver/txt.cpp +++ b/Code/Source/solver/txt.cpp @@ -321,7 +321,7 @@ void txt(Simulation* simulation, const bool init_write, const SolutionStates& so case OutputNameType::outGrp_divV: case OutputNameType::outGrp_J: case OutputNameType::outGrp_mises: - case OutputNameType::outGrp_MBFFlux: + case OutputNameType::outGrp_darcyFlux: post::all_post(simulation, tmpV, solutions, oGrp, iEq); break; diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 5247dc871..5a3ea1203 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1164,7 +1164,7 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b case OutputNameType::outGrp_stInv: case OutputNameType::outGrp_vortex: case OutputNameType::outGrp_Visc: - case OutputNameType::outGrp_MBFFlux: + case OutputNameType::outGrp_darcyFlux: post::post(simulation, msh, tmpV, solutions, oGrp, iEq); for (int a = 0; a < msh.nNo; a++) { int Ac = msh.gN(a); From 68fe6ae58887a445ce564727abc7370a36bb2042 Mon Sep 17 00:00:00 2001 From: Michael Date: Thu, 3 Sep 2026 12:08:48 -0700 Subject: [PATCH 10/28] Removing empty NURBS block and fixing recent PhysicalPropertyType changes --- Code/Source/solver/darcy.cpp | 5 ----- Code/Source/solver/set_bc.cpp | 4 ++-- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 2655640f1..30319b5c9 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -99,11 +99,6 @@ void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& s continue; } - // Update shape function for NURBS - if (lM.eType == ElementType::NRB) { - //CALL NRMNNX(lm, e) - } - // Create local copies for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a, e); diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 2e32b85d5..9a4ef9355 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -1592,9 +1592,9 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const int cDmn_local = all_fun::domain(com_mod, com_mod.msh[iM], cEq, lFa.gE(0)); double rho = eq.dmn[cDmn_local].prop.at( - consts::PhysicalProperyType::fluid_density); + consts::PhysicalPropertyType::fluid_density); double beta = eq.dmn[cDmn_local].prop.at( - consts::PhysicalProperyType::backflow_stab); + consts::PhysicalPropertyType::backflow_stab); double A = lFa.area; if (A > 0.0) { double u_n = Q_3D / A; // face-averaged normal velocity (< 0) From 59f6444a54daa60a5e203cefe5a86176f44486ce Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 1 Sep 2026 15:03:39 -0700 Subject: [PATCH 11/28] docs: clarify Darcy pressure and flux outputs --- Code/Source/solver/post.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index fb403b050..1c39aae8b 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -1011,8 +1011,8 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S } } - // Darcy Flux calculation - // + // Darcy flux, derived from pressure: + // q = -(K/mu) grad(p). } else if (outGrp == OutputNameType::outGrp_darcyFlux) { const double permeability = eq.dmn[cDmn].prop[PhysicalPropertyType::darcy_permeability]; From 6189ba623efe167979be58591eab983b577958ac Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 1 Sep 2026 15:38:32 -0700 Subject: [PATCH 12/28] refactor: scope Darcy material parameter names --- Code/Source/solver/Parameters.cpp | 9 +- Code/Source/solver/Parameters.h | 9 +- Code/Source/solver/README.md | 8 +- Code/Source/solver/consts.h | 2 +- Code/Source/solver/fluid.cpp | 96 +++++++++---------- Code/Source/solver/fluid.h | 9 +- Code/Source/solver/read_files.cpp | 4 +- Code/Source/solver/set_equation_props.h | 2 +- .../fluid/driven_cavity_2d_porous/solver.xml | 6 +- 9 files changed, 72 insertions(+), 73 deletions(-) diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 730dea29d..ce2fdba61 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2053,8 +2053,9 @@ DomainParameters::DomainParameters() { set_parameter("Penalty_parameter", 0.0, !required, penalty_parameter); set_parameter("Poisson_ratio", 0.3, !required, poisson_ratio); - set_parameter("Permeability", 1e-15, !required, darcy_permeability); - set_parameter("Media_compressibility", 0.0, !required, darcy_media_compressibility); + set_parameter("Darcy_permeability", 1e-15, !required, darcy_permeability); + set_parameter("Darcy_media_compressibility", 0.0, !required, + darcy_media_compressibility); set_parameter("Darcy_fluid_viscosity", 1.0, !required, darcy_fluid_viscosity); set_parameter("Relative_tolerance", 1e-4, !required, relative_tolerance); @@ -2064,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( diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index f7ce9582d..c6895e5ed 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1640,8 +1640,10 @@ class DomainParameters : public ParameterLists Parameter darcy_permeability; Parameter darcy_media_compressibility; Parameter darcy_fluid_viscosity; - // Inverse of Darcy permeability. Default value of 0.0 for Navier-Stokes and non-zero for Navier-Stokes-Brinkman - Parameter inverse_darcy_permeability; + + // Inverse permeability K^{-1} used in the Brinkman drag term + // mu K^{-1} u. A value of zero disables Brinkman drag. + Parameter brinkman_inverse_permeability; }; /// @brief The RemesherParameters class stores parameters for the @@ -1762,9 +1764,6 @@ class EquationParameters : public ParameterLists // and only then is the mesh equation solved. Parameter explicit_geometric_coupling; - // Inverse of Darcy permeability. Default value of 0.0 for Navier-Stokes and non-zero for Navier-Stokes-Brinkman - Parameter inverse_darcy_permeability; - // Sub-element parameters. // std::vector body_forces; diff --git a/Code/Source/solver/README.md b/Code/Source/solver/README.md index d11378e35..89cdfde10 100644 --- a/Code/Source/solver/README.md +++ b/Code/Source/solver/README.md @@ -178,8 +178,8 @@ C++ functions are defined within a `namespace` defined for each Fortran file. Fo - [ fs::get_thood_fs(com_mod, fs, lM, vmsStab, 1) ](#) - [ nn::gnn(fs[1].eNoN, nsd, nsd, Nx, xql, Nqx, Jac, ksix) ](#) - [ 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_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_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, brinkman_inverse_permeability)](#) - If nsd=2 - [ trilinos_doassem_(const_cast(eNoN), ptr.data(), lK.data(), lR.data())](#) - If using Trilinos - [ lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR)](#do_assem) - If not using Trilinos - [ set_bc::set_bc_neu(com_mod, cm_mod, Yg, Dg) ](#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(eNoN), ptr.data(), lK.data(), lR.data())` - If using Trilinos diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index 636daa9d4..d08145f84 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -418,7 +418,7 @@ enum class PhysicalPropertyType 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_media_compressibility = 17, darcy_fluid_viscosity = 18 diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index 907ee5a85..c0ef6eb65 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -553,7 +553,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s continue; } - double K_inverse_darcy_permeability = eq.dmn[cDmn].prop.at(PhysicalPropertyType::inverse_darcy_permeability); + double brinkman_inverse_permeability = eq.dmn[cDmn].prop.at(PhysicalPropertyType::brinkman_inverse_permeability); // Update shape functions for NURBS if (lM.eType == ElementType::NRB) { @@ -668,14 +668,14 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s urisValveVelTermTotal = urisValveVelTermTotalEl.rcol(g); } 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, + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, brinkman_inverse_permeability, urisFactorTotal, urisValveVelTermTotal); } else if (nsd == 2) { auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); 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); + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, brinkman_inverse_permeability); } } // g: loop @@ -734,14 +734,14 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s urisValveVelTermTotal = urisValveVelTermTotalEl.rcol(g); } fluid_3d_c(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, + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, brinkman_inverse_permeability, urisFactorTotal, urisValveVelTermTotal); } else if (nsd == 2) { auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); fluid_2d_c(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); + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, brinkman_inverse_permeability); } } // g: loop @@ -766,7 +766,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s void fluid_2d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability) + const Array& bfl, Array& lR, Array3& lK, double brinkman_inverse_permeability) { using namespace consts; @@ -991,7 +991,7 @@ void fluid_2d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e double kT = 4.0 * pow(ctM/dt,2.0); // If we consider the NSB model, we need to add an extra term inside the computation for the stab parameter - kT = kT + pow(K_inverse_darcy_permeability*mu/rho, 2.0); + kT = kT + pow(brinkman_inverse_permeability*mu/rho, 2.0); double kU = u(0)*u(0)*Kxi(0,0) + u(1)*u(0)*Kxi(1,0) + u(0)*u(1)*Kxi(0,1) + u(1)*u(1)*Kxi(1,1); double kS = Kxi(0,0)*Kxi(0,0) + Kxi(1,0)*Kxi(1,0) + Kxi(0,1)*Kxi(0,1) + Kxi(1,1)*Kxi(1,1); @@ -1009,12 +1009,12 @@ void fluid_2d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e rS(0) = mu_x(0)*es(0,0) + mu_x(1)*es(1,0) + mu*d2u2(0); rS(1) = mu_x(0)*es(0,1) + mu_x(1)*es(1,1) + mu*d2u2(1); - up(0) = -tauM*(rho*rV(0) + px(0) - rS(0) + mu*K_inverse_darcy_permeability*u(0)); - up(1) = -tauM*(rho*rV(1) + px(1) - rS(1) + mu*K_inverse_darcy_permeability*u(1)); + up(0) = -tauM*(rho*rV(0) + px(0) - rS(0) + mu*brinkman_inverse_permeability*u(0)); + up(1) = -tauM*(rho*rV(1) + px(1) - rS(1) + mu*brinkman_inverse_permeability*u(1)); for (int a = 0; a < eNoNw; a++) { double uNx = u(0)*Nwx(0,a) + u(1)*Nwx(1,a); - T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a)) + mu_x(0)*Nwx(0,a) + mu_x(1)*Nwx(1,a) - mu*K_inverse_darcy_permeability*Nw(a); + T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a)) + mu_x(0)*Nwx(0,a) + mu_x(1)*Nwx(1,a) - mu*brinkman_inverse_permeability*Nw(a); updu(0,0,a) = mu_x(0)*Nwx(0,a) + d2u2(0)*mu_g*esNx(0,a) + T1; updu(1,0,a) = mu_x(1)*Nwx(0,a) + d2u2(1)*mu_g*esNx(0,a); @@ -1086,7 +1086,7 @@ void fluid_2d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability) + const Array& bfl, Array& lR, Array3& lK, double brinkman_inverse_permeability) { using namespace consts; @@ -1279,7 +1279,7 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e double kT = 4.0 * pow(ctM/dt,2.0); // If we consider the NSB model, we need to add an extra term inside the computation for the stab parameter - kT = kT + pow(K_inverse_darcy_permeability*mu/rho, 2.0); + kT = kT + pow(brinkman_inverse_permeability*mu/rho, 2.0); double kU = u(0)*u(0)*Kxi(0,0) + u(1)*u(0)*Kxi(1,0) + u(0)*u(1)*Kxi(0,1) + u(1)*u(1)*Kxi(1,1); double kS = Kxi(0,0)*Kxi(0,0) + Kxi(1,0)*Kxi(1,0) + Kxi(0,1)*Kxi(0,1) + Kxi(1,1)*Kxi(1,1); @@ -1301,8 +1301,8 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // up[i] = ith component of u_prime (where u_prime = fine-scale velocity in VMS) = -tau_M / rho * ith component of momentum PDE residual (not weak form residual) Vector up(2); - up(0) = -tauM*(rho*rV(0) + px(0) - rS(0) + mu*K_inverse_darcy_permeability*u(0)); - up(1) = -tauM*(rho*rV(1) + px(1) - rS(1) + mu*K_inverse_darcy_permeability*u(1)); + up(0) = -tauM*(rho*rV(0) + px(0) - rS(0) + mu*brinkman_inverse_permeability*u(0)); + up(1) = -tauM*(rho*rV(1) + px(1) - rS(1) + mu*brinkman_inverse_permeability*u(1)); // tauC = rho * tau_C; tauB = rho * tau_bar; pa = pressure - rho * tau_C * divergence of velocity double tauC, tauB, pa; @@ -1367,7 +1367,7 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e uaNx(a) = uNx(a); } - T1 = -rho*uNx(a) + mu*(Nwxx(0,a) + Nwxx(1,a)) + mu_x(0)*Nwx(0,a) + mu_x(1)*Nwx(1,a) - mu*K_inverse_darcy_permeability*Nw(a); + T1 = -rho*uNx(a) + mu*(Nwxx(0,a) + Nwxx(1,a)) + mu_x(0)*Nwx(0,a) + mu_x(1)*Nwx(1,a) - mu*brinkman_inverse_permeability*Nw(a); updu(0,0,a) = mu_x(0)*Nwx(0,a) + d2u2(0)*mu_g*esNx(0,a) + T1; updu(1,0,a) = mu_x(1)*Nwx(0,a) + d2u2(1)*mu_g*esNx(0,a); @@ -1392,7 +1392,7 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // dRm_a1/du_b1 // derivative of x-component of momentum (weak form) residual with respect to the x-component of (the acceleration at the next time step) - lK(0,a,b) = lK(0,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); + lK(0,a,b) = lK(0,a,b) + mu*brinkman_inverse_permeability*wl*Nw(b)*Nw(a); T2 = mu*rM(1,0) + tauC*rM(0,1) + esNx(0,a)*mu_g*esNx(1,b) - rho*tauM*uaNx(a)*updu(1,0,b); @@ -1411,7 +1411,7 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // dRm_a2/du_b2 // derivative of y-component of momentum (weak form) residual with respect to the y-component of (the acceleration at the next time step) lK(4,a,b) = lK(4,a,b) + wl*(T2 + T1); - lK(4,a,b) = lK(4,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); + lK(4,a,b) = lK(4,a,b) + mu*brinkman_inverse_permeability*wl*Nw(b)*Nw(a); } } @@ -1432,8 +1432,8 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Residual contribution Birkman term // Local residue for (int a = 0; a < eNoNw; a++) { - lR(0,a) = lR(0,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u(0)+up(0)); - lR(1,a) = lR(1,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u(1)+up(1)); + lR(0,a) = lR(0,a) + mu*brinkman_inverse_permeability*w*Nw(a)*(u(0)+up(0)); + lR(1,a) = lR(1,a) + mu*brinkman_inverse_permeability*w*Nw(a)*(u(1)+up(1)); } } @@ -1443,7 +1443,7 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, + const Array& bfl, Array& lR, Array3& lK, double brinkman_inverse_permeability, const double urisFactorTotal, const Vector& urisValveVelTermTotal) { #define n_debug_fluid3d_c @@ -1655,7 +1655,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e double kT = 4.0 * pow(ctM/dt,2.0); // If we consider the NSB model, we need to add an extra term inside the computation for the stab parameter - kT = kT + pow(K_inverse_darcy_permeability*mu/rho, 2.0); + kT = kT + pow(brinkman_inverse_permeability*mu/rho, 2.0); // In case of unfitted RIS, compute the delta function at the quad point, // add the additional value to the stabilization param @@ -1682,24 +1682,24 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e rS[1] = mu_x[0]*es[0][1] + mu_x[1]*es[1][1] + mu_x[2]*es[2][1] + mu*d2u2[1]; rS[2] = mu_x[0]*es[0][2] + mu_x[1]*es[1][2] + mu_x[2]*es[2][2] + mu*d2u2[2]; - // up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability*u[0]); - // up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability*u[1]); - // up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability*u[2]); + // up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*brinkman_inverse_permeability*u[0]); + // up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*brinkman_inverse_permeability*u[1]); + // up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*brinkman_inverse_permeability*u[2]); - up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability*u[0] + up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*brinkman_inverse_permeability*u[0] + urisFactorTotal*u[0] - urisValveVelTermTotal[0]); - up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability*u[1] + up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*brinkman_inverse_permeability*u[1] + urisFactorTotal*u[1] - urisValveVelTermTotal[1]); - up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability*u[2] + up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*brinkman_inverse_permeability*u[2] + urisFactorTotal*u[2] - urisValveVelTermTotal[2]); for (int a = 0; a < eNoNw; a++) { double uNx = u[0]*Nwx(0,a) + u[1]*Nwx(1,a) + u[2]*Nwx(2,a); - // T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a); + // T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*brinkman_inverse_permeability*Nw(a); T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - - mu*K_inverse_darcy_permeability*Nw(a) + - mu*brinkman_inverse_permeability*Nw(a) - urisFactorTotal*Nw(a); updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; @@ -1768,7 +1768,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, + const Array& bfl, Array& lR, Array3& lK, double brinkman_inverse_permeability, const double urisFactorTotal, const Vector& urisValveVelTermTotal) { #define n_debug_fluid_3d_m @@ -2001,7 +2001,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e double kT = 4.0 * pow(ctM/dt,2.0); // If we consider the NSB model, we need to add an extra term inside the computation for the stab parameter - kT = kT + pow(K_inverse_darcy_permeability*mu/rho, 2.0); + kT = kT + pow(brinkman_inverse_permeability*mu/rho, 2.0); // In case of unfitted RIS, compute the delta function at the quad point, // add the additional value to the stabilization param @@ -2035,15 +2035,15 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e rS[2] = mu_x[0]*es[0][2] + mu_x[1]*es[1][2] + mu_x[2]*es[2][2] + mu*d2u2[2]; double up[3] = {}; - // up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability * u[0]); - // up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability * u[1]); - // up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability * u[2]); + // up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*brinkman_inverse_permeability * u[0]); + // up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*brinkman_inverse_permeability * u[1]); + // up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*brinkman_inverse_permeability * u[2]); - up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability * u[0] + up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*brinkman_inverse_permeability * u[0] + urisFactorTotal * u[0] - urisValveVelTermTotal[0]); - up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability * u[1] + up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*brinkman_inverse_permeability * u[1] + urisFactorTotal * u[1] - urisValveVelTermTotal[1]); - up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability * u[2] + up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*brinkman_inverse_permeability * u[2] + urisFactorTotal * u[2] - urisValveVelTermTotal[2]); double tauC, tauB, pa; @@ -2121,11 +2121,11 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e uaNx[a] = uNx[a]; } - // T1 = -rho*uNx[a] + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a); + // T1 = -rho*uNx[a] + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*brinkman_inverse_permeability*Nw(a); T1 = -rho*uNx[a] + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - - mu*K_inverse_darcy_permeability*Nw(a) + - mu*brinkman_inverse_permeability*Nw(a) - urisFactorTotal*Nw(a); updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; @@ -2161,8 +2161,8 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // dRm_a1/du_b1 double T2 = (mu + tauC)*rM[0][0] + esNx[0][a]*mu_g*esNx[0][b] - rho*tauM*uaNx[a]*updu[0][0][b]; lK(0,a,b) = lK(0,a,b) + wl*(T2 + T1); - // lK(0,a,b) = lK(0,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); - lK(0,a,b) = lK(0,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) + // lK(0,a,b) = lK(0,a,b) + mu*brinkman_inverse_permeability*wl*Nw(b)*Nw(a); + lK(0,a,b) = lK(0,a,b) + mu*brinkman_inverse_permeability*wl*Nw(b)*Nw(a) + urisFactorTotal*wl*Nw(b)*Nw(a); // dRm_a1/du_b2 @@ -2180,8 +2180,8 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // dRm_a2/du_b2 T2 = (mu + tauC)*rM[1][1] + esNx[1][a]*mu_g*esNx[1][b] - rho*tauM*uaNx[a]*updu[1][1][b]; lK(5,a,b) = lK(5,a,b) + wl*(T2 + T1); - // lK(5,a,b) = lK(5,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); - lK(5,a,b) = lK(5,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) + // lK(5,a,b) = lK(5,a,b) + mu*brinkman_inverse_permeability*wl*Nw(b)*Nw(a); + lK(5,a,b) = lK(5,a,b) + mu*brinkman_inverse_permeability*wl*Nw(b)*Nw(a) + urisFactorTotal*wl*Nw(b)*Nw(a); // dRm_a2/du_b3 @@ -2199,8 +2199,8 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // dRm_a3/du_b3; T2 = (mu + tauC)*rM[2][2] + esNx[2][a]*mu_g*esNx[2][b] - rho*tauM*uaNx[a]*updu[2][2][b]; lK(10,a,b) = lK(10,a,b) + wl*(T2 + T1); - // lK(10,a,b) = lK(10,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); - lK(10,a,b) = lK(10,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) + // lK(10,a,b) = lK(10,a,b) + mu*brinkman_inverse_permeability*wl*Nw(b)*Nw(a); + lK(10,a,b) = lK(10,a,b) + mu*brinkman_inverse_permeability*wl*Nw(b)*Nw(a) + urisFactorTotal*wl*Nw(b)*Nw(a); //dmsg << "lK(10,a,b): " << lK(10,a,b); } @@ -2226,11 +2226,11 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Residual contribution Birkman term // Local residue for (int a = 0; a < eNoNw; a++) { - lR(0,a) = lR(0,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[0]+up[0]) + lR(0,a) = lR(0,a) + mu*brinkman_inverse_permeability*w*Nw(a)*(u[0]+up[0]) + w*Nw(a)*(urisFactorTotal*u[0] - urisValveVelTermTotal[0]); - lR(1,a) = lR(1,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[1]+up[1]) + lR(1,a) = lR(1,a) + mu*brinkman_inverse_permeability*w*Nw(a)*(u[1]+up[1]) + w*Nw(a)*(urisFactorTotal*u[1] - urisValveVelTermTotal[1]); - lR(2,a) = lR(2,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[2]+up[2]) + lR(2,a) = lR(2,a) + mu*brinkman_inverse_permeability*w*Nw(a)*(u[2]+up[2]) + w*Nw(a)*(urisFactorTotal*u[2] - urisValveVelTermTotal[2]); } diff --git a/Code/Source/solver/fluid.h b/Code/Source/solver/fluid.h index af2191206..c1432314d 100644 --- a/Code/Source/solver/fluid.h +++ b/Code/Source/solver/fluid.h @@ -26,23 +26,23 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const SolutionStates& s void fluid_2d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeabilityx); + Array& lR, Array3& lK, double brinkman_inverse_permeability); void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability); + Array& lR, Array3& lK, double brinkman_inverse_permeability); void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability, + Array& lR, Array3& lK, double brinkman_inverse_permeability, const double urisFactorTotal, const Vector& urisValveVelTotal); void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability, + Array& lR, Array3& lK, double brinkman_inverse_permeability, const double urisFactorTotal, const Vector& urisValveVelTotal); void get_viscosity(const ComMod& com_mod, const dmnType& lDmn, double& gamma, double& mu, double& mu_s, double& mu_x); @@ -50,4 +50,3 @@ void get_viscosity(const ComMod& com_mod, const dmnType& lDmn, double& gamma, do }; #endif - diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 30e7fcaff..25d704977 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1577,8 +1577,8 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& rtmp = domain_params->source_term.value(); break; - case PhysicalPropertyType::inverse_darcy_permeability: - rtmp = domain_params->inverse_darcy_permeability.value(); + case PhysicalPropertyType::brinkman_inverse_permeability: + rtmp = domain_params->brinkman_inverse_permeability.value(); break; case PhysicalPropertyType::darcy_permeability: diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index f68a5569b..a7c314b3a 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -191,7 +191,7 @@ SetEquationPropertiesMapType set_equation_props = { propL[0][0] = PhysicalPropertyType::fluid_density; propL[1][0] = PhysicalPropertyType::backflow_stab; - propL[2][0] = PhysicalPropertyType::inverse_darcy_permeability; + propL[2][0] = PhysicalPropertyType::brinkman_inverse_permeability; propL[3][0] = PhysicalPropertyType::f_x; propL[4][0] = PhysicalPropertyType::f_y; diff --git a/tests/cases/fluid/driven_cavity_2d_porous/solver.xml b/tests/cases/fluid/driven_cavity_2d_porous/solver.xml index d2ff7d3d2..4baa98921 100644 --- a/tests/cases/fluid/driven_cavity_2d_porous/solver.xml +++ b/tests/cases/fluid/driven_cavity_2d_porous/solver.xml @@ -56,7 +56,7 @@ 1.111111111111111e-06 - 1e11 + 1e11 @@ -65,7 +65,7 @@ 1e-06 - 0.0 + 0.0 @@ -120,4 +120,4 @@ - \ No newline at end of file + From 1d9d56206ba4e735554827cf5a1b26f27e6aaef1 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 1 Sep 2026 16:10:36 -0700 Subject: [PATCH 13/28] fix: validate Darcy material parameters --- Code/Source/solver/darcy.cpp | 42 ++++++++++++++++++++++++- Code/Source/solver/darcy.h | 4 ++- Code/Source/solver/read_files.cpp | 1 + Code/Source/solver/set_equation_props.h | 4 +++ 4 files changed, 49 insertions(+), 2 deletions(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 30319b5c9..8db236991 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -1,10 +1,15 @@ #include "darcy.h" +#include "Core/Exception.h" #include "all_fun.h" #include "mat_fun.h" #include "nn.h" #include "utils.h" +#include +#include +#include + namespace darcy { /* This code implements the Darcy Equation for 2D and 3D @@ -48,6 +53,41 @@ namespace darcy { * - \f$ \Gamma \f$ : Domain boundary */ +void validate_material_properties(const dmnType& domain) +{ + using consts::PhysicalPropertyType; + + const auto validate = [&domain](const char* name, const double value, + const bool is_valid, + const char* requirement) { + if (is_valid) { + return; + } + + std::ostringstream message; + message << std::setprecision(std::numeric_limits::max_digits10) + << "Darcy domain " << domain.Id << " has invalid " << name + << " value " << value << "; expected " << requirement << "."; + svmp::raise(message.str()); + }; + + const double permeability = + domain.prop.at(PhysicalPropertyType::darcy_permeability); + validate("Darcy_permeability", permeability, permeability > 0.0, + "a value greater than zero"); + + const double viscosity = + domain.prop.at(PhysicalPropertyType::darcy_fluid_viscosity); + validate("Darcy_fluid_viscosity", viscosity, viscosity > 0.0, + "a value greater than zero"); + + const double compressibility = + domain.prop.at(PhysicalPropertyType::darcy_media_compressibility); + validate("Darcy_media_compressibility", compressibility, + compressibility >= 0.0, + "a value greater than or equal to zero"); +} + void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR) { for (int a = 0; a < eNoN; a++) { @@ -265,4 +305,4 @@ void darcy_3d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR); void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& solutions); @@ -22,4 +24,4 @@ namespace darcy { const Array& al, const Array& yl, Array& lR, Array3& lK); } -#endif //DARCY_H \ No newline at end of file +#endif //DARCY_H diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 25d704977..601b3fd91 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -10,6 +10,7 @@ #include "ActiveStress.h" #include "all_fun.h" #include "consts.h" +#include "darcy.h" #include "IonicModel.h" #include "read_msh.h" #include "vtk_xml.h" diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index a7c314b3a..491fa4953 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -292,6 +292,10 @@ SetEquationPropertiesMapType set_equation_props = { read_domain(simulation, eq_params, lEq, propL); + for (const auto& domain : lEq.dmn) { + darcy::validate_material_properties(domain); + } + nDOP = {2,1,1,0}; outPuts = {OutputNameType::out_darcyPressure, OutputNameType::out_darcyFlux}; From 99c80d1f56a92f4fdf3a613cbeafe14c606b1afd Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 1 Sep 2026 18:16:47 -0700 Subject: [PATCH 14/28] docs: align the Darcy equation with assembly --- Code/Source/solver/darcy.cpp | 54 ++++++++++++++---------------------- 1 file changed, 21 insertions(+), 33 deletions(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 8db236991..d66ceee26 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -12,45 +12,33 @@ namespace darcy { /* - This code implements the Darcy Equation for 2D and 3D - problems in perfusion of porous media. + This code implements the Darcy equation for perfusion of porous media. ------------------------------------------------------------- Assumptions: - - Homogeneous Permeability - - Homogeneous Density + - Material coefficients are homogeneous within each solver domain - Isotropic Permeability - Assumptions of Stokes Flow - - Steady-State ------------------------------------------------------------- - * Strong form of the Single-Compartment Darcy equation: - * \f[ u = -K\nabla P \f] - * \f[ \nabla \cdot u = \beta_0(P_{source} - P) - \beta_1(P - P_{sink}) \f] - * Note: See equations 8(a)/(b) in https://doi.org/10.1007/s10439-020-02681-z - * - * Combined Strong form: - * \f[ -\nabla \cdot (K \nabla P) - \beta_0(P_{source} - P) + \beta_1(P - P_{sink}) = 0 \f] - * + * Strong form of the Darcy equation assembled by these kernels: + * \f[ + * \rho \beta \frac{\partial p}{\partial t} + * - \nabla \cdot \left(\frac{\rho K}{\mu}\nabla p\right) + * = \rho s. + * \f] + * * Where: - * - \f$ u \f$ : Darcy flux - * - \f$ K \f$ : Permeability tensor - * - \f$ P \f$ : Pressure - * - \f$ P_{source} \f$ : Source pressure (e.g., arterial pressure) - * - \f$ P_{sink} \f$ : Sink pressure (e.g., venous/extraction pressure) - * - \f$ \beta_0 \f$ : Source coupling term (describes conductance of flow entering myocardium) - * - \f$ \beta_1 \f$ : Sink coupling term (describes conductance of flow exiting myocardium) - * - * ------------------------------------------------------------- - * - * Weak form of the Single-Compartment Darcy equation: - * \f[ -\int_{\Omega} (\nabla q \cdot \nabla P) d\Omega - \lambda \int_{\Omega} q P d\Omega = \int_{\Omega} q F d\Omega - \int_{\Gamma} q (\nabla P \cdot n) d\Gamma \f] - * - * Where: - * - \f$ q \f$ : Test function - * - \f$ \lambda \f$ : \f$ \frac{\beta_0 + \beta_1}{K} \f$ - * - \f$ F \f$ : \f$ -\frac{\beta_0 P_{source} + \beta_1 P_{sink}}{K} \f$ - * - \f$ n \f$ : Normal vector to the boundary - * - \f$ \Omega \f$ : Computational domain - * - \f$ \Gamma \f$ : Domain boundary + * - \f$ p \f$ is pressure. + * - \f$ \rho \f$ is fluid density. + * - \f$ \beta \f$ is the porous-media compressibility. + * - \f$ K \f$ is intrinsic permeability. + * - \f$ \mu \f$ is dynamic viscosity. + * - \f$ s \f$ is the volumetric source specified by `Source_term`. + * + * Density is retained in the mass-conservation form because it may differ + * between solver domains, although it cancels after normalization within a + * single homogeneous domain. Viscosity is retained because \f$K/\mu\f$ is + * mobility; folding \f$\mu\f$ into \f$K\f$ would change the meaning and units + * of the permeability input. */ void validate_material_properties(const dmnType& domain) From 33449862931c2099d8b22275e1253c0d34c9020c Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 1 Sep 2026 18:20:33 -0700 Subject: [PATCH 15/28] docs: state the Darcy single-field formulation --- Code/Source/solver/darcy.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index d66ceee26..562ce732a 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -26,6 +26,12 @@ namespace darcy { * = \rho s. * \f] * + * After pressure is solved, Darcy velocity is + * evaluated as the derived field + * \f[ + * \boldsymbol{q} = -\frac{K}{\mu}\nabla p. + * \f] + * * Where: * - \f$ p \f$ is pressure. * - \f$ \rho \f$ is fluid density. From 5b094b41a141071597bf6f4fd268747c94131862 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 1 Sep 2026 18:22:51 -0700 Subject: [PATCH 16/28] docs: retain permeability in the diffusion operator --- Code/Source/solver/darcy.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 562ce732a..57510aca5 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -40,6 +40,10 @@ namespace darcy { * - \f$ \mu \f$ is dynamic viscosity. * - \f$ s \f$ is the volumetric source specified by `Source_term`. * + * The current implementation accepts a positive scalar \f$K\f$ for each + * solver domain. Spatially heterogeneous or tensor-valued permeability is + * currently not implemented. + * * Density is retained in the mass-conservation form because it may differ * between solver domains, although it cancels after normalization within a * single homogeneous domain. Viscosity is retained because \f$K/\mu\f$ is From 7b4ca1c77745d216fdb93d64c75cdc6c82ef6390 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 1 Sep 2026 18:27:29 -0700 Subject: [PATCH 17/28] docs: define Darcy model quantities and ranges --- Code/Source/solver/darcy.cpp | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 57510aca5..bc7bb97ad 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -32,13 +32,19 @@ namespace darcy { * \boldsymbol{q} = -\frac{K}{\mu}\nabla p. * \f] * - * Where: - * - \f$ p \f$ is pressure. - * - \f$ \rho \f$ is fluid density. - * - \f$ \beta \f$ is the porous-media compressibility. - * - \f$ K \f$ is intrinsic permeability. - * - \f$ \mu \f$ is dynamic viscosity. - * - \f$ s \f$ is the volumetric source specified by `Source_term`. + * Model quantities and ranges: + * - \f$p\f$: pressure. + * - \f$\boldsymbol{q}\f$: derived Darcy velocity. + * - \f$K\f$: scalar permeability (`Darcy_permeability`, default + * \f$10^{-15}\f$, \f$K > 0\f$). + * - \f$\mu\f$: dynamic viscosity (`Darcy_fluid_viscosity`, default 1, + * \f$\mu > 0\f$). + * - \f$\rho\f$: fluid density (`Fluid_density`, default 0.5, + * \f$\rho > 0\f$). + * - \f$\beta\f$: compressibility (`Darcy_media_compressibility`, default 0, + * \f$\beta \ge 0\f$). + * - \f$s\f$: volumetric source (`Source_term`, default 0), constant per + * domain. * * The current implementation accepts a positive scalar \f$K\f$ for each * solver domain. Spatially heterogeneous or tensor-valued permeability is @@ -79,6 +85,10 @@ void validate_material_properties(const dmnType& domain) validate("Darcy_fluid_viscosity", viscosity, viscosity > 0.0, "a value greater than zero"); + const double density = domain.prop.at(PhysicalPropertyType::fluid_density); + validate("Fluid_density", density, density > 0.0, + "a value greater than zero"); + const double compressibility = domain.prop.at(PhysicalPropertyType::darcy_media_compressibility); validate("Darcy_media_compressibility", compressibility, From a785875c48a3bc55b470e8ff12a70e39807d55dd Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 1 Sep 2026 18:33:05 -0700 Subject: [PATCH 18/28] docs: document the Darcy namespace and references --- Code/Source/solver/darcy.cpp | 45 ------------------- Code/Source/solver/darcy.h | 85 +++++++++++++++++++++++++++++++++--- 2 files changed, 80 insertions(+), 50 deletions(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index bc7bb97ad..91e967e3d 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -11,51 +11,6 @@ #include namespace darcy { -/* - This code implements the Darcy equation for perfusion of porous media. - ------------------------------------------------------------- - Assumptions: - - Material coefficients are homogeneous within each solver domain - - Isotropic Permeability - - Assumptions of Stokes Flow - ------------------------------------------------------------- - * Strong form of the Darcy equation assembled by these kernels: - * \f[ - * \rho \beta \frac{\partial p}{\partial t} - * - \nabla \cdot \left(\frac{\rho K}{\mu}\nabla p\right) - * = \rho s. - * \f] - * - * After pressure is solved, Darcy velocity is - * evaluated as the derived field - * \f[ - * \boldsymbol{q} = -\frac{K}{\mu}\nabla p. - * \f] - * - * Model quantities and ranges: - * - \f$p\f$: pressure. - * - \f$\boldsymbol{q}\f$: derived Darcy velocity. - * - \f$K\f$: scalar permeability (`Darcy_permeability`, default - * \f$10^{-15}\f$, \f$K > 0\f$). - * - \f$\mu\f$: dynamic viscosity (`Darcy_fluid_viscosity`, default 1, - * \f$\mu > 0\f$). - * - \f$\rho\f$: fluid density (`Fluid_density`, default 0.5, - * \f$\rho > 0\f$). - * - \f$\beta\f$: compressibility (`Darcy_media_compressibility`, default 0, - * \f$\beta \ge 0\f$). - * - \f$s\f$: volumetric source (`Source_term`, default 0), constant per - * domain. - * - * The current implementation accepts a positive scalar \f$K\f$ for each - * solver domain. Spatially heterogeneous or tensor-valued permeability is - * currently not implemented. - * - * Density is retained in the mass-conservation form because it may differ - * between solver domains, although it cancels after normalization within a - * single homogeneous domain. Viscosity is retained because \f$K/\mu\f$ is - * mobility; folding \f$\mu\f$ into \f$K\f$ would change the meaning and units - * of the permeability input. -*/ void validate_material_properties(const dmnType& domain) { diff --git a/Code/Source/solver/darcy.h b/Code/Source/solver/darcy.h index 6822d70ac..01888991e 100644 --- a/Code/Source/solver/darcy.h +++ b/Code/Source/solver/darcy.h @@ -1,25 +1,100 @@ -// -// This code implements the Darcy Equation for 2D and 3D -// problems in perfusion of porus media. -// - #ifndef DARCY_H #define DARCY_H #include "ComMod.h" #include "SolutionStates.h" +/** + * @brief Pressure-based Darcy flow in porous media. + * + * This namespace implements the Darcy equation for perfusion of porous media + * with intrinsic dimensions 2 and 3. Material coefficients are homogeneous + * within each solver domain, permeability is isotropic, and Stokes-flow + * assumptions apply. The assembled pressure strong form is + * \f[ + * \rho \beta \frac{\partial p}{\partial t} + * - \nabla \cdot \left(\frac{\rho K}{\mu}\nabla p\right) + * = \rho s. + * \f] + * + * This discretizes the pressure-only strong form. + * Velocity is not an independent unknown. + * After pressure is solved, Darcy velocity is + * evaluated as the derived field + * \f[ + * \boldsymbol{q} = -\frac{K}{\mu}\nabla p. + * \f] + * + * The model quantities and their admissible ranges are: + * - \f$p\f$: pressure unknown. + * - \f$\boldsymbol{q}\f$: derived Darcy velocity. + * - \f$K\f$: configured intrinsic scalar permeability. `Darcy_permeability` + * defaults to \f$10^{-15}\f$ and must satisfy \f$K > 0\f$. + * - \f$\mu\f$: configured dynamic viscosity. `Darcy_fluid_viscosity` defaults + * to 1 and must satisfy \f$\mu > 0\f$. + * - \f$\rho\f$: configured reference fluid density. `Fluid_density` defaults + * to 0.5 and must satisfy \f$\rho > 0\f$. + * - \f$\beta\f$: configured storage/compressibility. + * `Darcy_media_compressibility` defaults to 0 and must satisfy + * \f$\beta \ge 0\f$. + * - \f$s\f$: configured volumetric source provided by `Source_term`; it + * defaults to 0 and is constant within each configured domain. + * + * @par Cardiovascular porous-flow context + * The following works describe future multi-compartment and microcirculation + * model extensions than the single-compartment formulation implemented here: + * - C. Michler et al., "A computationally efficient framework for the + * simulation of cardiac perfusion using a multi-compartment Darcy + * porous-media flow model," DOI + * 10.1002/cnm.2520. + * - G. Montino Pelagi et al., "Modeling cardiac microcirculation for the + * simulation of coronary flow and 3D myocardial perfusion," DOI + * 10.1007/s10237-024-01873-z. + */ namespace darcy { + /// Validate the configured Darcy material coefficients for a domain. + /// @param[in] domain Solver domain whose material properties are checked. void validate_material_properties(const dmnType& domain); + /// Assemble a Darcy boundary contribution into the element residual. + /// @param[in] com_mod Common solver state retained for the common assembly interface. + /// @param[in] eNoN Number of element nodes. + /// @param[in] w Weighted boundary quadrature measure. + /// @param[in] N Shape-function values at the quadrature point. + /// @param[in] h Prescribed boundary flux contribution. + /// @param[in,out] lR Element residual. void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR); + /// Assemble Darcy volume contributions for all supported elements in a mesh. + /// @param[in,out] com_mod Common solver state and assembly interface. + /// @param[in] lM Mesh whose Darcy elements are assembled. + /// @param[in] solutions Solution states used for element-local fields. void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& solutions); + /// Assemble the residual and tangent for an intrinsic two-dimensional element. + /// @param[in] com_mod Common solver state. + /// @param[in] eNoN Number of element nodes. + /// @param[in] w Weighted volume quadrature measure. + /// @param[in] N Shape-function values at the quadrature point. + /// @param[in] Nx Mapped spatial shape-function derivatives. + /// @param[in] al Element-local pressure rates. + /// @param[in] yl Element-local pressure state. + /// @param[in,out] lR Element residual. + /// @param[in,out] lK Element tangent matrix. void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, const Array& al, const Array& yl, Array& lR, Array3& lK); + /// Assemble the residual and tangent for an intrinsic three-dimensional element. + /// @param[in] com_mod Common solver state. + /// @param[in] eNoN Number of element nodes. + /// @param[in] w Weighted volume quadrature measure. + /// @param[in] N Shape-function values at the quadrature point. + /// @param[in] Nx Mapped spatial shape-function derivatives. + /// @param[in] al Element-local pressure rates. + /// @param[in] yl Element-local pressure state. + /// @param[in,out] lR Element residual. + /// @param[in,out] lK Element tangent matrix. void darcy_3d(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, const Array& al, const Array& yl, Array& lR, Array3& lK); } From 5b7d58245e27725b1e4d4ca6de9969f446aeea3e Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 1 Sep 2026 22:35:06 -0700 Subject: [PATCH 19/28] fix: reject embedded-line Darcy use --- Code/Source/solver/darcy.cpp | 16 ++++++++++++---- Code/Source/solver/darcy.h | 7 +++++++ Code/Source/solver/post.cpp | 5 +++++ 3 files changed, 24 insertions(+), 4 deletions(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 91e967e3d..29fcd3300 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -51,6 +51,15 @@ void validate_material_properties(const dmnType& domain) "a value greater than or equal to zero"); } +void validate_element_support(const mshType& mesh) +{ + if (mesh.lFib) { + svmp::raise( + "The Darcy equation supports only 2D and 3D meshes; lFib marks an " + "embedded one-dimensional mesh, which is not supported."); + } +} + void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const double h, Array& lR) { for (int a = 0; a < eNoN; a++) { @@ -60,6 +69,8 @@ void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector& res, const S bool FSIeq = false; auto& eq = com_mod.eq[iEq]; + if (outGrp == OutputNameType::outGrp_darcyFlux) { + darcy::validate_element_support(lM); + } + if (eq.phys == EquationType::phys_FSI) { FSIeq = true; } From 8eb5ec5f17c99b80d3e0258fe9e9ecb3517ff803 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 1 Sep 2026 22:43:36 -0700 Subject: [PATCH 20/28] refactor: clarify Darcy flux reconstruction --- Code/Source/solver/darcy.h | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Code/Source/solver/darcy.h b/Code/Source/solver/darcy.h index e413192f0..fab832561 100644 --- a/Code/Source/solver/darcy.h +++ b/Code/Source/solver/darcy.h @@ -40,6 +40,13 @@ * - \f$s\f$: configured volumetric source provided by `Source_term`; it * defaults to 0 and is constant within each configured domain. * + * @par Darcy flux output + * On supported two- and three-dimensional meshes, the pressure gradient is + * reconstructed in the mesh coordinates and the derived Darcy flux is + * \f[ + * \boldsymbol{q} = -\frac{K}{\mu}\nabla p. + * \f] + * * @par Cardiovascular porous-flow context * The following works describe future multi-compartment and microcirculation * model extensions than the single-compartment formulation implemented here: From a87e72981ffa6de9ee9b9c16d0c95f9285e90d7a Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 10 Sep 2026 09:40:29 -0700 Subject: [PATCH 21/28] fix: initialize pressure using equation DOFs --- Code/Source/solver/initialize.cpp | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index c082867c9..89dad332e 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -992,9 +992,15 @@ void zero_init(Simulation* simulation, SolutionStates& solutions) #ifdef debug_zero_init dmsg << "Initialize Yo to provided P solution"; #endif - for (int a = 0; a < com_mod.tnNo; a++) { - for (int i = 0; i < nsd; i++) { - Yo(nsd,a) = com_mod.Pinit(a); + for (const auto& eq : com_mod.eq) { + const bool is_darcy = eq.phys == consts::EquationType::phys_darcy; + // Skip equations without a pressure unknown. + if (!is_darcy && eq.dof != nsd + 1) { + continue; + } + const int pressure_dof = eq.s + (is_darcy ? 0 : nsd); + for (int a = 0; a < com_mod.tnNo; ++a) { + Yo(pressure_dof,a) = com_mod.Pinit(a); } } } From 3d3d45744618de996c43960e510ac2d3418d9c98 Mon Sep 17 00:00:00 2001 From: Michael Date: Fri, 11 Sep 2026 12:16:22 -0700 Subject: [PATCH 22/28] First set of comments: dropping element validation in inappropriate places and adding parameter dimensions to documentation --- Code/Source/solver/darcy.cpp | 1 - Code/Source/solver/darcy.h | 16 ++++++++-------- Code/Source/solver/post.cpp | 4 ---- 3 files changed, 8 insertions(+), 13 deletions(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 29fcd3300..2db4c258b 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -69,7 +69,6 @@ void b_darcy(ComMod& com_mod, const int eNoN, const double w, const Vector 0\f$. - * - \f$\mu\f$: configured dynamic viscosity. `Darcy_fluid_viscosity` defaults + * - \f$\mu\f$: configured dynamic viscosity [M/L/T]. `Darcy_fluid_viscosity` defaults * to 1 and must satisfy \f$\mu > 0\f$. - * - \f$\rho\f$: configured reference fluid density. `Fluid_density` defaults + * - \f$\rho\f$: configured reference fluid density [M/L^3]. `Fluid_density` defaults * to 0.5 and must satisfy \f$\rho > 0\f$. - * - \f$\beta\f$: configured storage/compressibility. - * `Darcy_media_compressibility` defaults to 0 and must satisfy + * - \f$\beta\f$: configured storage/compressibility [L*T^2/M]. + * `Darcy_compressibility` defaults to 0 and must satisfy * \f$\beta \ge 0\f$. - * - \f$s\f$: configured volumetric source provided by `Source_term`; it + * - \f$s\f$: configured volumetric source provided by `Source_term` [1/T]; it * defaults to 0 and is constant within each configured domain. * * @par Darcy flux output diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index c2cb24095..5f45d272e 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -843,10 +843,6 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S bool FSIeq = false; auto& eq = com_mod.eq[iEq]; - if (outGrp == OutputNameType::outGrp_darcyFlux) { - darcy::validate_element_support(lM); - } - if (eq.phys == EquationType::phys_FSI) { FSIeq = true; } From a9bbd48d3d70a24d613fff64c8d7524097e66a13 Mon Sep 17 00:00:00 2001 From: Michael Date: Fri, 11 Sep 2026 12:49:42 -0700 Subject: [PATCH 23/28] Renaming Darcy_media_compressibility to Darcy_compressibility & adding exceptions from svmp namespace --- Code/Source/solver/Parameters.cpp | 4 ++-- Code/Source/solver/Parameters.h | 2 +- Code/Source/solver/consts.h | 2 +- Code/Source/solver/darcy.cpp | 12 ++++++------ Code/Source/solver/read_files.cpp | 4 ++-- Code/Source/solver/set_equation_props.h | 2 +- 6 files changed, 13 insertions(+), 13 deletions(-) diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index ce2fdba61..9ff3eca0f 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2054,8 +2054,8 @@ DomainParameters::DomainParameters() { set_parameter("Poisson_ratio", 0.3, !required, poisson_ratio); set_parameter("Darcy_permeability", 1e-15, !required, darcy_permeability); - set_parameter("Darcy_media_compressibility", 0.0, !required, - darcy_media_compressibility); + 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); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index b89db62e7..f97336688 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1662,7 +1662,7 @@ class DomainParameters : public ParameterLists Parameter time_step_for_integration; Parameter darcy_permeability; - Parameter darcy_media_compressibility; + Parameter darcy_compressibility; Parameter darcy_fluid_viscosity; // Inverse permeability K^{-1} used in the Brinkman drag term diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index d08145f84..507ae8bdb 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -420,7 +420,7 @@ enum class PhysicalPropertyType ctau_C = 14, brinkman_inverse_permeability = 15, darcy_permeability = 16, - darcy_media_compressibility = 17, + darcy_compressibility = 17, darcy_fluid_viscosity = 18 }; diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index 2db4c258b..fc0caccc1 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -45,8 +45,8 @@ void validate_material_properties(const dmnType& domain) "a value greater than zero"); const double compressibility = - domain.prop.at(PhysicalPropertyType::darcy_media_compressibility); - validate("Darcy_media_compressibility", compressibility, + domain.prop.at(PhysicalPropertyType::Darcy_compressibility); + validate("Darcy_compressibility", compressibility, compressibility >= 0.0, "a value greater than or equal to zero"); } @@ -134,7 +134,7 @@ void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& s auto Nx_g = lM.Nx.slice(g); nn::gnn(eNoN, nsd, insd, Nx_g, xl, Nx, Jac, ksix); if (utils::is_zero(Jac)) { - throw std::runtime_error( + svmp::raise( "[construct_darcy] Jacobian for element " + std::to_string(e) + " is < 0."); } } @@ -147,7 +147,7 @@ void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& s } else if (insd == 2) { darcy_2d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); } else { - throw std::runtime_error("[construct_darcy] insd must be 2 or 3."); + throw std::InvalidArgumentException("[construct_darcy] insd must be 2 or 3."); } } @@ -176,7 +176,7 @@ void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vectordarcy_permeability.value(); break; - case PhysicalPropertyType::darcy_media_compressibility: - rtmp = domain_params->darcy_media_compressibility.value(); + case PhysicalPropertyType::darcy_compressibility: + rtmp = domain_params->darcy_compressibility.value(); break; case PhysicalPropertyType::darcy_fluid_viscosity: diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index 491fa4953..18e04ce97 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -287,7 +287,7 @@ SetEquationPropertiesMapType set_equation_props = { propL[0][0] = PhysicalPropertyType::darcy_permeability; propL[1][0] = PhysicalPropertyType::source_term; propL[2][0] = PhysicalPropertyType::fluid_density; - propL[3][0] = PhysicalPropertyType::darcy_media_compressibility; + propL[3][0] = PhysicalPropertyType::darcy_compressibility; propL[4][0] = PhysicalPropertyType::darcy_fluid_viscosity; read_domain(simulation, eq_params, lEq, propL); From a68f7c40b9034a1d7f890ac90166be8a99b69c3c Mon Sep 17 00:00:00 2001 From: Michael Date: Fri, 11 Sep 2026 12:57:53 -0700 Subject: [PATCH 24/28] Fixing instance of darcy_compressibility --- Code/Source/solver/darcy.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index fc0caccc1..f560a8a17 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -45,7 +45,7 @@ void validate_material_properties(const dmnType& domain) "a value greater than zero"); const double compressibility = - domain.prop.at(PhysicalPropertyType::Darcy_compressibility); + domain.prop.at(PhysicalPropertyType::darcy_compressibility); validate("Darcy_compressibility", compressibility, compressibility >= 0.0, "a value greater than or equal to zero"); @@ -176,7 +176,7 @@ void darcy_2d(ComMod& com_mod, const int eNoN, const double w, const Vector Date: Fri, 11 Sep 2026 13:03:19 -0700 Subject: [PATCH 25/28] fixing incomplete removal of runtime error call from std namespace --- Code/Source/solver/darcy.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index f560a8a17..beaa0460f 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -147,7 +147,7 @@ void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& s } else if (insd == 2) { darcy_2d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); } else { - throw std::InvalidArgumentException("[construct_darcy] insd must be 2 or 3."); + svmp::raise("[construct_darcy] insd must be 2 or 3."); } } From c1597932e731e76ba6239cd140a9565ba3212376 Mon Sep 17 00:00:00 2001 From: Michael Date: Fri, 11 Sep 2026 13:18:49 -0700 Subject: [PATCH 26/28] fixing svmp exception calls --- Code/Source/solver/darcy.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index beaa0460f..d0399ead3 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -134,8 +134,7 @@ void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& s auto Nx_g = lM.Nx.slice(g); nn::gnn(eNoN, nsd, insd, Nx_g, xl, Nx, Jac, ksix); if (utils::is_zero(Jac)) { - svmp::raise( - "[construct_darcy] Jacobian for element " + std::to_string(e) + " is < 0."); + svmp::InternalErrorException("[construct_darcy] Jacobian for element " + std::to_string(e) + " is < 0."); } } @@ -147,7 +146,7 @@ void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& s } else if (insd == 2) { darcy_2d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); } else { - svmp::raise("[construct_darcy] insd must be 2 or 3."); + svmp::raise("[construct_darcy] insd must be 2 or 3."); } } From d9c4cc8c6463ed6dae92bbe4572a77f50cb270b5 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Mon, 14 Sep 2026 19:27:40 -0700 Subject: [PATCH 27/28] L2 projection and volume-weighted flux. --- Code/Source/solver/post.cpp | 69 +++++++++++++++++++++++++++++++++---- 1 file changed, 63 insertions(+), 6 deletions(-) diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 5f45d272e..89481bafd 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -3,7 +3,9 @@ #include "post.h" +#include "Core/Exception.h" #include "FE/Common/FEException.h" +#include "FE/Math/DenseLinearAlgebra.h" #include "all_fun.h" #include "darcy.h" #include "fluid.h" @@ -15,7 +17,10 @@ #include "shells.h" #include "utils.h" #include "vtk_xml.h" +#include +#include #include +#include namespace post { @@ -880,6 +885,13 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S Array Nx(nsd,eNoN); Vector N(eNoN); + // Linear-simplex flux is constant; lumped recovery is already exact. + const bool project_flux = outGrp == OutputNameType::outGrp_darcyFlux && + lM.eType != ElementType::TRI3 && lM.eType != ElementType::TET4; + // DenseLinearAlgebra expects row-major matrices and multiple right-hand sides. + std::vector mass(project_flux ? eNoN * eNoN : 0); + std::vector flux_rhs(project_flux ? eNoN * nsd : 0); + int insd = nsd; if (lM.lFib) { insd = 1; @@ -890,6 +902,10 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S if (cDmn == -1) { continue; } + std::fill(mass.begin(), mass.end(), 0.0); + std::fill(flux_rhs.begin(), flux_rhs.end(), 0.0); + double element_volume = 0.0; + if (lM.eType == ElementType::NRB) { // CALL NRBNNX(lM, e) } @@ -1098,12 +1114,48 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S throw std::runtime_error("Error in the post() function."); } - // Mapping Tau into the nodes by assembling it into a local vector - for (int a = 0; a < eNoN; a++) { - int Ac = lM.IEN(a,e); - sA(Ac) = sA(Ac) + w*N(a); - for (int i = 0; i < maxNSD; i++) { - sF(i,Ac) = sF(i,Ac) + w*N(a)*lRes(i); + if (project_flux) { + // M_ab = integral(N_a N_b), B_ai = integral(N_a q_i). + // Consistent projection avoids zero lumped weights on affine nodes. + element_volume += w; + for (int a = 0; a < eNoN; ++a) { + const double weighted_shape = w * N(a); + for (int b = 0; b < eNoN; ++b) { + mass[a * eNoN + b] += weighted_shape * N(b); + } + for (int i = 0; i < nsd; ++i) { + flux_rhs[a * nsd + i] += weighted_shape * lRes(i); + } + } + } else { + // Mapping Tau into the nodes by assembling it into a local vector + for (int a = 0; a < eNoN; a++) { + int Ac = lM.IEN(a,e); + sA(Ac) = sA(Ac) + w*N(a); + for (int i = 0; i < maxNSD; i++) { + sF(i,Ac) = sF(i,Ac) + w*N(a)*lRes(i); + } + } + } + } + + if (project_flux) { + svmp::check( + std::isfinite(element_volume) && element_volume > 0.0, + "Darcy flux projection requires a positive element volume."); + + // Solve (M / volume) Q = B for volume-weighted flux directly. + // Normalizing M keeps the pivot tolerance independent of element size. + for (double& value : mass) { + value /= element_volume; + } + svmp::FE::math::factor_dense_matrix(mass, eNoN, "Darcy flux mass matrix").solve_in_place(flux_rhs, nsd); + + for (int a = 0; a < eNoN; ++a) { + const int Ac = lM.IEN(a,e); + sA(Ac) += element_volume; + for (int i = 0; i < nsd; ++i) { + sF(i,Ac) += flux_rhs[a * nsd + i]; } } } @@ -1114,6 +1166,11 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S for (int a = 0; a < lM.nNo; a++) { int Ac = lM.gN(a); + if (project_flux) { + svmp::check( + std::isfinite(sA(Ac)) && sA(Ac) > 0.0, + "Darcy flux projection requires a positive nodal volume weight."); + } for (int i = 0; i < maxNSD; i++) { res(i,a) = sF(i,Ac) / sA(Ac); } From ae99df477a2015b554271c9076b42a8542d92323 Mon Sep 17 00:00:00 2001 From: Michael Date: Thu, 17 Sep 2026 12:32:22 -0700 Subject: [PATCH 28/28] dropping context tag --- Code/Source/solver/darcy.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp index d0399ead3..f14500b1e 100644 --- a/Code/Source/solver/darcy.cpp +++ b/Code/Source/solver/darcy.cpp @@ -134,7 +134,7 @@ void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& s auto Nx_g = lM.Nx.slice(g); nn::gnn(eNoN, nsd, insd, Nx_g, xl, Nx, Jac, ksix); if (utils::is_zero(Jac)) { - svmp::InternalErrorException("[construct_darcy] Jacobian for element " + std::to_string(e) + " is < 0."); + svmp::InternalErrorException("Jacobian for element " + std::to_string(e) + " is < 0."); } } @@ -146,7 +146,7 @@ void construct_darcy(ComMod& com_mod, const mshType& lM, const SolutionStates& s } else if (insd == 2) { darcy_2d(com_mod, eNoN, w, N, Nx, al, yl, lR, lK); } else { - svmp::raise("[construct_darcy] insd must be 2 or 3."); + svmp::raise("insd must be 2 or 3."); } }