diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index f4f239fed..3f57b4daf 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -212,6 +212,7 @@ set(CSRCS cmm.h cmm.cpp consts.h consts.cpp contact.h contact.cpp + darcy.h darcy.cpp distribute.h distribute.cpp eq_assem.h eq_assem.cpp fluid.h fluid.cpp diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index d1b5f8685..fb7084d1c 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 f0922a258..9ff3eca0f 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2052,6 +2052,11 @@ DomainParameters::DomainParameters() { set_parameter("Penalty_parameter", 0.0, !required, penalty_parameter); set_parameter("Poisson_ratio", 0.3, !required, poisson_ratio); + + set_parameter("Darcy_permeability", 1e-15, !required, darcy_permeability); + set_parameter("Darcy_compressibility", 0.0, !required, + darcy_compressibility); + set_parameter("Darcy_fluid_viscosity", 1.0, !required, darcy_fluid_viscosity); set_parameter("Relative_tolerance", 1e-4, !required, relative_tolerance); set_parameter("Shell_thickness", 0.0, !required, shell_thickness); @@ -2060,8 +2065,8 @@ DomainParameters::DomainParameters() { set_parameter("Time_step_for_integration", 0.0, !required, time_step_for_integration); - set_parameter("Inverse_darcy_permeability", 0.0, !required, - inverse_darcy_permeability); + set_parameter("Brinkman_inverse_permeability", 0.0, !required, + brinkman_inverse_permeability); // Ionic model parameters. IonicModelFactory::visit( diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index aecffc9db..f97336688 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1661,8 +1661,13 @@ class DomainParameters : public ParameterLists Parameter source_term; Parameter time_step_for_integration; - // Inverse of Darcy permeability. Default value of 0.0 for Navier-Stokes and non-zero for Navier-Stokes-Brinkman - Parameter inverse_darcy_permeability; + Parameter darcy_permeability; + Parameter darcy_compressibility; + Parameter darcy_fluid_viscosity; + + // Inverse permeability K^{-1} used in the Brinkman drag term + // mu K^{-1} u. A value of zero disables Brinkman drag. + Parameter brinkman_inverse_permeability; }; /// @brief The RemesherParameters class stores parameters for the @@ -1783,9 +1788,6 @@ class EquationParameters : public ParameterLists // and only then is the mesh equation solved. Parameter 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/cmm.cpp b/Code/Source/solver/cmm.cpp index 785a3429b..a6e8f827c 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; @@ -429,10 +429,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 ; @@ -444,7 +444,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; @@ -497,7 +497,7 @@ void cmm_stiffness(ComMod& com_mod, const Array& Nxi, const Array& Nxi, const Array 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..507ae8bdb 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -287,9 +287,11 @@ enum class EquationType phys_CMM = 209, phys_CEP = 210, phys_ustruct = 211, // Nonlinear elastodynamics using mixed VMS-stabilized formulation - phys_stokes = 212 + phys_stokes = 212, + phys_darcy = 213 }; +constexpr auto Equation_darcy = EquationType::phys_darcy; constexpr auto Equation_CMM = EquationType::phys_CMM; constexpr auto Equation_CEP = EquationType::phys_CEP; constexpr auto Equation_fluid = EquationType::phys_fluid; @@ -347,6 +349,7 @@ enum class OutputNameType { outGrp_activeTensionFibers = 529, outGrp_activeTensionSheets = 530, outGrp_activeTensionNormal = 531, + outGrp_darcyFlux = 532, out_velocity = 599, out_pressure = 598, @@ -380,7 +383,9 @@ enum class OutputNameType { out_fibStretchRate = 570, out_activeTensionFibers = 569, out_activeTensionSheets = 568, - out_activeTensionNormal = 567 + out_activeTensionNormal = 567, + out_darcyPressure = 566, + out_darcyFlux = 565 }; /// @brief Simulation output file types. @@ -396,7 +401,7 @@ extern const std::map output_type_name_to_type; /// @brief Possible physical properties. Current maxNPror is 20. // -enum class PhysicalProperyType +enum class PhysicalPropertyType { NA = 0, fluid_density = 1, @@ -413,7 +418,10 @@ enum class PhysicalProperyType shell_thickness = 12, ctau_M = 13, // stabilization coeffs. for USTRUCT (momentum, continuity) ctau_C = 14, - inverse_darcy_permeability = 15 + brinkman_inverse_permeability = 15, + darcy_permeability = 16, + darcy_compressibility = 17, + darcy_fluid_viscosity = 18 }; enum class PreconditionerType diff --git a/Code/Source/solver/darcy.cpp b/Code/Source/solver/darcy.cpp new file mode 100644 index 000000000..f14500b1e --- /dev/null +++ b/Code/Source/solver/darcy.cpp @@ -0,0 +1,277 @@ +#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 { + +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 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_compressibility); + validate("Darcy_compressibility", compressibility, + compressibility >= 0.0, + "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++) { + 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; + const int insd = nsd; + #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; + } + + // 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)) { + svmp::InternalErrorException("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 { + svmp::raise("insd must be 2 or 3."); + } + } + + eq.linear_algebra->assemble(com_mod, eNoN, ptr, lK, lR); + } +} + +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(PhysicalPropertyType::darcy_permeability); + double source = dmn.prop.at(PhysicalPropertyType::source_term); + double beta_0 = dmn.prop.at(PhysicalPropertyType::darcy_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; + + #ifdef debug_darcy_2d + dmsg << "k: " << k; + dmsg << "source: " << source; + dmsg << "T1: " << T1; + dmsg << "i: " << i; + dmsg << "wl: " << wl; + #endif + + double p_dot = 0.0; + Vector Px(nsd); + + for (int a = 0; a < eNoN; a++) { + p_dot = p_dot + N(a)*al(i,a); + Px(0) = Px(0) + Nx(0,a)*yl(i,a); + 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 * 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) + + 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(PhysicalPropertyType::darcy_permeability); + double source = dmn.prop.at(PhysicalPropertyType::source_term); + double beta_0 = dmn.prop.at(PhysicalPropertyType::darcy_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; + + #ifdef debug_darcy_3d + dmsg << "k: " << k; + dmsg << "source: " << source; + dmsg << "T1: " << T1; + dmsg << "i: " << i; + dmsg << "wl: " << wl; + #endif + + double p_dot = 0.0; + Vector Px(nsd); + + for (int a = 0; a < eNoN; a++) { + p_dot = p_dot + N(a) * al(i,a); + Px(0) = Px(0) + Nx(0,a) * yl(i,a); + 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 * 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) + + Nx(2,a)*Nx(2,b))); + } + } +} + +} diff --git a/Code/Source/solver/darcy.h b/Code/Source/solver/darcy.h new file mode 100644 index 000000000..bde5482cd --- /dev/null +++ b/Code/Source/solver/darcy.h @@ -0,0 +1,116 @@ +#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 [M/L/T^2]. + * - \f$\boldsymbol{q}\f$: derived Darcy velocity [L/T]. + * - \f$K\f$: configured intrinsic scalar permeability [L^2]. `Darcy_permeability` + * defaults to \f$10^{-15}\f$ and must satisfy \f$K > 0\f$. + * - \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 [M/L^3]. `Fluid_density` defaults + * to 0.5 and must satisfy \f$\rho > 0\f$. + * - \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` [1/T]; 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: + * - 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); + + /// Reject mesh types that Darcy assembly and flux output do not implement. + /// @param[in] mesh Mesh whose element type is checked. + /// @note `mshType::lFib` denotes a one-dimensional mesh embedded in the + /// ambient geometry, not a myocardial material fiber direction. Darcy is + /// currently limited to intrinsic dimensions 2 and 3. + void validate_element_support(const mshType& mesh); + + /// 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); +} + +#endif //DARCY_H diff --git a/Code/Source/solver/eq_assem.cpp b/Code/Source/solver/eq_assem.cpp index 463a178ca..1358e23d0 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" @@ -110,6 +111,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; @@ -406,6 +411,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/fluid.cpp b/Code/Source/solver/fluid.cpp index 0d58a2725..c0ef6eb65 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -65,7 +65,7 @@ void b_fluid(ComMod& com_mod, const int eNoN, const double w, const Vector& 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; @@ -794,12 +794,12 @@ void fluid_2d_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]; 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; @@ -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; @@ -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; @@ -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 @@ -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; @@ -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 @@ -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; @@ -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/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& 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) { @@ -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); } } } diff --git a/Code/Source/solver/l_elas.cpp b/Code/Source/solver/l_elas.cpp index a29c658fb..70ad66fbb 100644 --- a/Code/Source/solver/l_elas.cpp +++ b/Code/Source/solver/l_elas.cpp @@ -168,12 +168,12 @@ void l_elas_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/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..89481bafd 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -3,8 +3,11 @@ #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" #include "fs.h" #include "initialize.h" @@ -14,7 +17,10 @@ #include "shells.h" #include "utils.h" #include "vtk_xml.h" +#include +#include #include +#include namespace post { @@ -853,7 +859,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); @@ -879,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; @@ -889,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) } @@ -960,7 +977,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); @@ -980,7 +997,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) { @@ -1011,6 +1028,28 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const S } } + // 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]; + const double viscosity = + eq.dmn[cDmn].prop[PhysicalPropertyType::darcy_fluid_viscosity]; + const double mobility = permeability / viscosity; + const int equation_index = eq.s; + + Vector grad_p(nsd); + + 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); + } + // Strain tensor invariants calculation // } else if (outGrp == OutputNameType::outGrp_stInv) { @@ -1075,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]; } } } @@ -1091,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); } @@ -1260,8 +1340,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 08899995f..6c8c02070 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" @@ -206,8 +207,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); } @@ -1513,59 +1514,59 @@ 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: - if (lEq.phys == EquationType::phys_CMM) { + case PhysicalPropertyType::fluid_density: + if (lEq.phys == EquationType::phys_CMM || lEq.phys == EquationType::phys_darcy) { rtmp = domain_params->fluid_density.value(); } else { rtmp = domain_params->density.value(); } 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,12 +1574,24 @@ 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: - rtmp = domain_params->inverse_darcy_permeability.value(); + case PhysicalPropertyType::brinkman_inverse_permeability: + rtmp = domain_params->brinkman_inverse_permeability.value(); + break; + + case PhysicalPropertyType::darcy_permeability: + rtmp = domain_params->darcy_permeability.value(); + break; + + case PhysicalPropertyType::darcy_compressibility: + rtmp = domain_params->darcy_compressibility.value(); + break; + + case PhysicalPropertyType::darcy_fluid_viscosity: + rtmp = domain_params->darcy_fluid_viscosity.value(); break; } @@ -1728,7 +1741,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; @@ -2240,8 +2253,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 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) 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..18e04ce97 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::brinkman_inverse_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); @@ -274,6 +274,36 @@ 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] = PhysicalPropertyType::darcy_permeability; + propL[1][0] = PhysicalPropertyType::source_term; + propL[2][0] = PhysicalPropertyType::fluid_density; + propL[3][0] = PhysicalPropertyType::darcy_compressibility; + propL[4][0] = PhysicalPropertyType::darcy_fluid_viscosity; + + 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}; + + // Set solver parameters. + read_ls(simulation, eq_params, SolverType::lSolver_CG, lEq); +} }, + + //---------------------------// // phys_FSI // //---------------------------// @@ -291,48 +321,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. @@ -386,13 +416,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); @@ -425,20 +455,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}; @@ -463,14 +493,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); @@ -503,11 +533,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); @@ -539,14 +569,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); @@ -594,15 +624,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 35cfbd8c4..b5f6a8b41 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_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/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 6598ef6b8..4bcd14d89 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_darcyFlux: 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 34b8b22e8..be9addb88 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1165,7 +1165,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_darcyFlux: 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/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 + 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); }