diff --git a/Code/Source/solver/active_stress.cpp b/Code/Source/solver/ActiveStress.cpp similarity index 98% rename from Code/Source/solver/active_stress.cpp rename to Code/Source/solver/ActiveStress.cpp index c3ee67172..7db495dae 100644 --- a/Code/Source/solver/active_stress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress.h" +#include "ActiveStress.h" bool supports_active_stress(const consts::EquationType eq_type) { return eq_type == consts::EquationType::phys_struct || diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/ActiveStress.h similarity index 100% rename from Code/Source/solver/active_stress.h rename to Code/Source/solver/ActiveStress.h diff --git a/Code/Source/solver/active_stress_nash_panfilov.cpp b/Code/Source/solver/ActiveStressNashPanfilov.cpp similarity index 57% rename from Code/Source/solver/active_stress_nash_panfilov.cpp rename to Code/Source/solver/ActiveStressNashPanfilov.cpp index 31095da3e..caa4def7c 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.cpp +++ b/Code/Source/solver/ActiveStressNashPanfilov.cpp @@ -1,9 +1,9 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_nash_panfilov.h" +#include "ActiveStressNashPanfilov.h" -void NashPanfilov::read_model_specific_parameters( +void ActiveStressNashPanfilov::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { ActiveStressODE::read_model_specific_parameters(params); @@ -15,8 +15,8 @@ void NashPanfilov::read_model_specific_parameters( eta_T = params.get_scalar("eta_T"); } -void NashPanfilov::distribute_model_specific_parameters(const CmMod &cm_mod, - const cmType &cm) { +void ActiveStressNashPanfilov::distribute_model_specific_parameters( + const CmMod &cm_mod, const cmType &cm) { ActiveStressODE::distribute_model_specific_parameters(cm_mod, cm); cm.bcast(cm_mod, &epsilon_0); @@ -27,12 +27,14 @@ void NashPanfilov::distribute_model_specific_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &eta_T); } -void NashPanfilov::init_local(Vector &state) const { state[0] = 0.0; } +void ActiveStressNashPanfilov::init_local(Vector &state) const { + state[0] = 0.0; +} -Vector NashPanfilov::getf(const double t, const Vector &state, - const double calcium, - const double fiber_stretch, - const double fiber_stretch_rate) const { +Vector +ActiveStressNashPanfilov::getf(const double t, const Vector &state, + const double calcium, const double fiber_stretch, + const double fiber_stretch_rate) const { Vector f(1); const double epsilon = @@ -44,10 +46,9 @@ Vector NashPanfilov::getf(const double t, const Vector &state, return f; } -double -NashPanfilov::compute_active_tension_local(const Vector &state, - const double fiber_stretch) const { +double ActiveStressNashPanfilov::compute_active_tension_local( + const Vector &state, const double fiber_stretch) const { return state[0]; } -REGISTER_ACTIVE_STRESS_MODEL("NashPanfilov", NashPanfilov); \ No newline at end of file +REGISTER_ACTIVE_STRESS_MODEL("NashPanfilov", ActiveStressNashPanfilov); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_nash_panfilov.h b/Code/Source/solver/ActiveStressNashPanfilov.h similarity index 97% rename from Code/Source/solver/active_stress_nash_panfilov.h rename to Code/Source/solver/ActiveStressNashPanfilov.h index b01d6f120..ac49775d8 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.h +++ b/Code/Source/solver/ActiveStressNashPanfilov.h @@ -4,7 +4,7 @@ #ifndef ACTIVE_STRESS_NASH_PANFILOV_H #define ACTIVE_STRESS_NASH_PANFILOV_H -#include "active_stress_ode.h" +#include "ActiveStressODE.h" /** * @brief Nash-Panfilov active stress model. @@ -36,7 +36,7 @@ * 1. [Nash, Panfilov (2004)](https://doi.org/10.1016/j.pbiomolbio.2004.01.016) * 2. [Goktepe, Kuhl (2009)](https://doi.org/10.1007/s00466-009-0434-z) */ -class NashPanfilov : public ActiveStressODE { +class ActiveStressNashPanfilov : public ActiveStressODE { public: /// Model label. static inline const std::string label = "NashPanfilov"; @@ -59,7 +59,7 @@ class NashPanfilov : public ActiveStressODE { /** * @brief Constructor. */ - NashPanfilov() : ActiveStressODE(1) {} + ActiveStressNashPanfilov() : ActiveStressODE(1) {} /** * @brief Construct an instance of model parameters. diff --git a/Code/Source/solver/active_stress_ode.cpp b/Code/Source/solver/ActiveStressODE.cpp similarity index 97% rename from Code/Source/solver/active_stress_ode.cpp rename to Code/Source/solver/ActiveStressODE.cpp index 99909cf7f..9b9ca9e17 100644 --- a/Code/Source/solver/active_stress_ode.cpp +++ b/Code/Source/solver/ActiveStressODE.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_ode.h" +#include "ActiveStressODE.h" void ActiveStressODE::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { diff --git a/Code/Source/solver/active_stress_ode.h b/Code/Source/solver/ActiveStressODE.h similarity index 99% rename from Code/Source/solver/active_stress_ode.h rename to Code/Source/solver/ActiveStressODE.h index 6ebfa7d09..470a48ce8 100644 --- a/Code/Source/solver/active_stress_ode.h +++ b/Code/Source/solver/ActiveStressODE.h @@ -4,7 +4,7 @@ #ifndef ACTIVE_STRESS_ODE_H #define ACTIVE_STRESS_ODE_H -#include "active_stress.h" +#include "ActiveStress.h" /** * @brief Abstract ODE-based active stress model. diff --git a/Code/Source/solver/active_stress_regazzoni.cpp b/Code/Source/solver/ActiveStressRegazzoni.cpp similarity index 90% rename from Code/Source/solver/active_stress_regazzoni.cpp rename to Code/Source/solver/ActiveStressRegazzoni.cpp index d7006c083..7c8e88dbf 100644 --- a/Code/Source/solver/active_stress_regazzoni.cpp +++ b/Code/Source/solver/ActiveStressRegazzoni.cpp @@ -1,14 +1,14 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_regazzoni.h" +#include "ActiveStressRegazzoni.h" #include "eigen3/Eigen/Dense" #include #include -void RegazzoniActiveStress::read_model_specific_parameters( +void ActiveStressRegazzoni::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { Kbasic = params.get_scalar("Kbasic"); Koff = params.get_scalar("Koff"); @@ -19,13 +19,15 @@ void RegazzoniActiveStress::read_model_specific_parameters( alphaKd = params.get_scalar("alphaKd"); if (alphaKd > 0.0) svmp::raise( - "RegazzoniActiveStress: alphaKd must be <= 0 (positive values reduce calcium " + "ActiveStressRegazzoni: alphaKd must be <= 0 (positive values reduce " + "calcium " "sensitivity with stretch, reversing length-dependent activation, " "and can produce a zero dissociation constant at physiological " "sarcomere lengths)."); SL0 = params.get_scalar("SL0"); ru_substep = params.get_scalar("ru_substep"); - kd_reference_sarcomere_length = params.get_scalar("kd_reference_sarcomere_length"); + kd_reference_sarcomere_length = + params.get_scalar("kd_reference_sarcomere_length"); r0 = params.get_scalar("r0"); alpha = params.get_scalar("alpha"); @@ -41,7 +43,7 @@ void RegazzoniActiveStress::read_model_specific_parameters( params.get_bool("Disable_force_strain_rate_feedback"); } -void RegazzoniActiveStress::distribute_model_specific_parameters( +void ActiveStressRegazzoni::distribute_model_specific_parameters( const CmMod &cm_mod, const cmType &cm) { cm.bcast(cm_mod, &Kbasic); cm.bcast(cm_mod, &Koff); @@ -66,14 +68,14 @@ void RegazzoniActiveStress::distribute_model_specific_parameters( cm.bcast(cm_mod, &disable_force_strain_rate_feedback_); } -void RegazzoniActiveStress::init_local(Vector &state) const { +void ActiveStressRegazzoni::init_local(Vector &state) const { for (unsigned int i = 0; i < n_state_variables; ++i) state[i] = 0.0; state[ru_index(0, 0, 0, 0)] = 1.0; // == state[0] } -void RegazzoniActiveStress::advance_time_step_local( +void ActiveStressRegazzoni::advance_time_step_local( const double t, const double dt, const double calcium, const double fiber_stretch, const double fiber_stretch_rate, Vector &state) const { @@ -114,7 +116,8 @@ void RegazzoniActiveStress::advance_time_step_local( } // Advance the crossbridge moments (entries 16-19) from the updated RU state. - // The velocity v = -dSL/dt / SL0 reduces to -d(lambda)/dt because SL = SL0 * lambda. + // The velocity v = -dSL/dt / SL0 reduces to -d(lambda)/dt because SL = SL0 * + // lambda. const double velocity = disable_force_strain_rate_feedback_ ? 0.0 : -fiber_stretch_rate; XBArray state_XB; @@ -134,7 +137,7 @@ void RegazzoniActiveStress::advance_time_step_local( state[xb_index(i)] = state_XB[i]; } -double RegazzoniActiveStress::compute_active_tension_local( +double ActiveStressRegazzoni::compute_active_tension_local( const Vector &state, const double fiber_stretch) const { const double sarcomere_length = SL0 * fiber_stretch; @@ -145,8 +148,8 @@ double RegazzoniActiveStress::compute_active_tension_local( fraction_single_overlap(sarcomere_length); } -RegazzoniActiveStress::RUArray -RegazzoniActiveStress::ru_transition_rates_tropomyosin() const { +ActiveStressRegazzoni::RUArray +ActiveStressRegazzoni::ru_transition_rates_tropomyosin() const { RUArray rates_T; for (int TL = 0; TL < 2; ++TL) for (int TR = 0; TR < 2; ++TR) { @@ -167,9 +170,9 @@ RegazzoniActiveStress::ru_transition_rates_tropomyosin() const { return rates_T; } -void RegazzoniActiveStress::ru_forward_euler_substep( - double dt, const RUArray &rates_T, - const BinaryPairArray &rates_C, RUArray &state_RU) const { +void ActiveStressRegazzoni::ru_forward_euler_substep( + double dt, const RUArray &rates_T, const BinaryPairArray &rates_C, + RUArray &state_RU) const { // Probability fluxes from central-unit transitions. RUArray flux_TC; // central tropomyosin RUArray flux_CC; // central troponin @@ -179,8 +182,7 @@ void RegazzoniActiveStress::ru_forward_euler_substep( for (int CC = 0; CC < 2; ++CC) { flux_TC[TL][TC][TR][CC] = state_RU[TL][TC][TR][CC] * rates_T[TL][TC][TR][CC]; - flux_CC[TL][TC][TR][CC] = - state_RU[TL][TC][TR][CC] * rates_C[CC][TC]; + flux_CC[TL][TC][TR][CC] = state_RU[TL][TC][TR][CC] * rates_C[CC][TC]; } // Effective transition rates of the boundary neighbours, obtained from the @@ -240,10 +242,8 @@ void RegazzoniActiveStress::ru_forward_euler_substep( flux_CC[TL][TC][TR][CC] + flux_CC[TL][TC][TR][1 - CC]); } -RegazzoniActiveStress::XBArray RegazzoniActiveStress::xb_implicit_update( - double dt, double velocity, - const RUArray &rates_T, - const RUArray &state_RU, +ActiveStressRegazzoni::XBArray ActiveStressRegazzoni::xb_implicit_update( + double dt, double velocity, const RUArray &rates_T, const RUArray &state_RU, const XBArray &state_XB) const { // Permissivity and the permissive/non-permissive probability fluxes from the // updated RU state. @@ -302,7 +302,8 @@ RegazzoniActiveStress::XBArray RegazzoniActiveStress::xb_implicit_update( return result; } -double RegazzoniActiveStress::fraction_single_overlap(double sarcomere_length) const { +double +ActiveStressRegazzoni::fraction_single_overlap(double sarcomere_length) const { const double SL = sarcomere_length; const double half_single_overlap = (LM - LB) * 0.5; @@ -317,4 +318,4 @@ double RegazzoniActiveStress::fraction_single_overlap(double sarcomere_length) c return 0.0; } -REGISTER_ACTIVE_STRESS_MODEL("Regazzoni", RegazzoniActiveStress); +REGISTER_ACTIVE_STRESS_MODEL("Regazzoni", ActiveStressRegazzoni); diff --git a/Code/Source/solver/active_stress_regazzoni.h b/Code/Source/solver/ActiveStressRegazzoni.h similarity index 80% rename from Code/Source/solver/active_stress_regazzoni.h rename to Code/Source/solver/ActiveStressRegazzoni.h index 27f7427ff..bcfe438ee 100644 --- a/Code/Source/solver/active_stress_regazzoni.h +++ b/Code/Source/solver/ActiveStressRegazzoni.h @@ -4,22 +4,23 @@ #ifndef ACTIVE_STRESS_REGAZZONI_H #define ACTIVE_STRESS_REGAZZONI_H -#include "active_stress.h" +#include "ActiveStress.h" #include /** * @brief Mean-field active stress model (implements the RDQ20-MF formulation). * - * This class implements the mean-field RDQ20-MF sarcomere model of cardiomyocyte - * force generation of Regazzoni, Dede', and Quarteroni (2020), described in [1] - * and validated against the authors' reference implementation [2]. The node-local state has 20 variables: 16 - * regulatory-unit (RU) probabilities (entries 0-15) describing the - * tropomyosin/troponin configuration of a triplet of neighbouring units, and 4 - * crossbridge (XB) moments (entries 16-19). The RU probabilities are advanced - * with an explicit forward-Euler substepping scheme — for every macro time step, - * a number of smaller sub-steps are taken to update the RU states — and the XB - * moments with one implicit-Euler step per time step; the active tension is then + * This class implements the mean-field RDQ20-MF sarcomere model of + * cardiomyocyte force generation of Regazzoni, Dede', and Quarteroni (2020), + * described in [1] and validated against the authors' reference implementation + * [2]. The node-local state has 20 variables: 16 regulatory-unit (RU) + * probabilities (entries 0-15) describing the tropomyosin/troponin + * configuration of a triplet of neighbouring units, and 4 crossbridge (XB) + * moments (entries 16-19). The RU probabilities are advanced with an explicit + * forward-Euler substepping scheme — for every macro time step, a number of + * smaller sub-steps are taken to update the RU states — and the XB moments with + * one implicit-Euler step per time step; the active tension is then * reconstructed from the XB first moments. * * The returned scalar active tension is @@ -27,24 +28,27 @@ * \Tact = a_\text{XB} \, (\mu_P^1 + \mu_N^1) \, \phi(SL)\;, * @f] * where @f$\mu_P^1@f$ and @f$\mu_N^1@f$ are the permissive and non-permissive - * first XB moments (state entries 17 and 19), @f$\phi(SL)@f$ is the single-overlap - * fraction of the sarcomere at sarcomere length @f$SL = SL_0 \, \fiberstretch@f$ - * (with @f$\fiberstretch@f$ the fiber stretch), and @f$a_\text{XB}@f$ is the tension - * upscaling factor. Because @f$\mu_P^1 + \mu_N^1@f$ and @f$\phi(SL)@f$ are - * dimensionless, @f$a_\text{XB}@f$ sets the units of the returned active tension. + * first XB moments (state entries 17 and 19), @f$\phi(SL)@f$ is the + * single-overlap fraction of the sarcomere at sarcomere length @f$SL = SL_0 \, + * \fiberstretch@f$ (with @f$\fiberstretch@f$ the fiber stretch), and + * @f$a_\text{XB}@f$ is the tension upscaling factor. Because @f$\mu_P^1 + + * \mu_N^1@f$ and @f$\phi(SL)@f$ are dimensionless, @f$a_\text{XB}@f$ sets the + * units of the returned active tension. * * **References**: - * 1. [Regazzoni, Dede', Quarteroni (2020)](https://doi.org/10.1371/journal.pcbi.1008294) - * 2. [F. Regazzoni, cardiac-activation reference implementation](https://github.com/FrancescoRegazzoni/cardiac-activation) + * 1. [Regazzoni, Dede', Quarteroni + * (2020)](https://doi.org/10.1371/journal.pcbi.1008294) + * 2. [F. Regazzoni, cardiac-activation reference + * implementation](https://github.com/FrancescoRegazzoni/cardiac-activation) * * @note Although this model is governed by a system of ODEs, it inherits from - * @c ActiveStress rather than @c ActiveStressODE because it requires a customized - * time-stepping scheme to handle the stiffness of the model. + * @c ActiveStress rather than @c ActiveStressODE because it requires a + * customized time-stepping scheme to handle the stiffness of the model. * * @todo Force-strain-rate feedback requires a stabilization strategy for robust * use in coupled electromechanics. This will be addressed in a follow-up PR. */ -class RegazzoniActiveStress : public ActiveStress { +class ActiveStressRegazzoni : public ActiveStress { public: /// Model label, used for factory registration and XML selection. static inline const std::string label = "Regazzoni"; @@ -123,7 +127,7 @@ class RegazzoniActiveStress : public ActiveStress { /** * @brief Constructor. */ - RegazzoniActiveStress() : ActiveStress(n_state_variables) {} + ActiveStressRegazzoni() : ActiveStress(n_state_variables) {} /** * @brief Construct an instance of model parameters. @@ -162,8 +166,8 @@ class RegazzoniActiveStress : public ActiveStress { * * Advances the RU probabilities (entries 0-15) with the forward-Euler * substepping scheme and then the XB moments (entries 16-19) with one - * implicit-Euler step, using the calcium, fiber stretch and fiber-stretch rate - * at the node. + * implicit-Euler step, using the calcium, fiber stretch and fiber-stretch + * rate at the node. */ virtual void advance_time_step_local(const double t, const double dt, const double calcium, @@ -184,7 +188,8 @@ class RegazzoniActiveStress : public ActiveStress { const double fiber_stretch) const override; private: - /// Array indexed over the four binary RU configuration variables (TL, TC, TR, CC). + /// Array indexed over the four binary RU configuration variables (TL, TC, TR, + /// CC). using RUArray = std::array, 2>, 2>, 2>; @@ -225,8 +230,7 @@ class RegazzoniActiveStress : public ActiveStress { * @param[in,out] state_RU The 16 RU-state probabilities, * indexed @c state_RU[TL][TC][TR][CC]. */ - void ru_forward_euler_substep(double dt, - const RUArray &rates_T, + void ru_forward_euler_substep(double dt, const RUArray &rates_T, const BinaryPairArray &rates_C, RUArray &state_RU) const; @@ -248,8 +252,7 @@ class RegazzoniActiveStress : public ActiveStress { * @return Updated crossbridge moments, ordered * @f$[\mu_P^0, \mu_P^1, \mu_N^0, \mu_N^1]@f$. */ - XBArray xb_implicit_update(double dt, double velocity, - const RUArray &rates_T, + XBArray xb_implicit_update(double dt, double velocity, const RUArray &rates_T, const RUArray &state_RU, const XBArray &state_XB) const; @@ -269,24 +272,27 @@ class RegazzoniActiveStress : public ActiveStress { /// @name RU model parameters /// @{ - double Kbasic; ///< Basic tropomyosin transition rate [1/time]. - double Koff; ///< Troponin unbinding rate [1/time]. - double Q; ///< Tropomyosin transition-rate asymmetry factor [-]. - double mu; ///< Calcium-binding cooperativity factor [-]. - double gamma; ///< Nearest-neighbour cooperativity factor [-]. - double Kd0; ///< Calcium dissociation constant at reference length [calcium]. - double alphaKd; ///< Length dependence of the dissociation constant [calcium/length]. - double SL0; ///< Reference sarcomere length [length]; maps stretch to length. + double Kbasic; ///< Basic tropomyosin transition rate [1/time]. + double Koff; ///< Troponin unbinding rate [1/time]. + double Q; ///< Tropomyosin transition-rate asymmetry factor [-]. + double mu; ///< Calcium-binding cooperativity factor [-]. + double gamma; ///< Nearest-neighbour cooperativity factor [-]. + double Kd0; ///< Calcium dissociation constant at reference length [calcium]. + double alphaKd; ///< Length dependence of the dissociation constant + ///< [calcium/length]. + double SL0; ///< Reference sarcomere length [length]; maps stretch to length. double ru_substep; ///< RU forward-Euler substep size [time]. /// Reference sarcomere length [length] used in the length-dependent /// dissociation constant (distinct from the parameter SL0). double kd_reference_sarcomere_length; - double r0; ///< Combined attachment-detachment rate at zero velocity [1/time]. + double r0; ///< Combined attachment-detachment rate at zero velocity [1/time]. double alpha; ///< Coefficient of |v| in r(v) = r0 + alpha * |v| [-]. - double mu0_fP; ///< Permissive influx into the zeroth-moment crossbridge state [1/time]. - double mu1_fP; ///< Permissive influx into the first-moment crossbridge state [1/time]. + double mu0_fP; ///< Permissive influx into the zeroth-moment crossbridge state + ///< [1/time]. + double mu1_fP; ///< Permissive influx into the first-moment crossbridge state + ///< [1/time]. double LA; ///< Thin-filament (actin) length [length]. double LM; ///< Thick-filament (myosin) length [length]. @@ -294,9 +300,10 @@ class RegazzoniActiveStress : public ActiveStress { /// Tension upscaling factor [stress]. /// - /// Because the crossbridge moments and the overlap fraction are dimensionless, - /// a_XB sets the stress unit of the returned active tension. It must be - /// expressed in the same stress unit as the mechanical configuration. + /// Because the crossbridge moments and the overlap fraction are + /// dimensionless, a_XB sets the stress unit of the returned active tension. + /// It must be expressed in the same stress unit as the mechanical + /// configuration. double a_XB; /// Controls force–strain-rate feedback in the XB update. diff --git a/Code/Source/solver/active_stress_uniform_steady.cpp b/Code/Source/solver/ActiveStressUniformSteady.cpp similarity index 57% rename from Code/Source/solver/active_stress_uniform_steady.cpp rename to Code/Source/solver/ActiveStressUniformSteady.cpp index 3c72855de..482ec1371 100644 --- a/Code/Source/solver/active_stress_uniform_steady.cpp +++ b/Code/Source/solver/ActiveStressUniformSteady.cpp @@ -1,16 +1,16 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_uniform_steady.h" +#include "ActiveStressUniformSteady.h" -void UniformSteadyActiveStress::read_model_specific_parameters( +void ActiveStressUniformSteady::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { value = params.get_scalar("Value"); } -void UniformSteadyActiveStress::distribute_model_specific_parameters( +void ActiveStressUniformSteady::distribute_model_specific_parameters( const CmMod &cm_mod, const cmType &cm) { cm.bcast(cm_mod, &value); } -REGISTER_ACTIVE_STRESS_MODEL("UniformSteady", UniformSteadyActiveStress); \ No newline at end of file +REGISTER_ACTIVE_STRESS_MODEL("UniformSteady", ActiveStressUniformSteady); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_uniform_steady.h b/Code/Source/solver/ActiveStressUniformSteady.h similarity index 94% rename from Code/Source/solver/active_stress_uniform_steady.h rename to Code/Source/solver/ActiveStressUniformSteady.h index ce5e2fd41..545a17549 100644 --- a/Code/Source/solver/active_stress_uniform_steady.h +++ b/Code/Source/solver/ActiveStressUniformSteady.h @@ -4,7 +4,7 @@ #ifndef ACTIVE_STRESS_UNIFORM_STEADY_H #define ACTIVE_STRESS_UNIFORM_STEADY_H -#include "active_stress.h" +#include "ActiveStress.h" /** * @brief Uniform and steady active stress model. @@ -15,7 +15,7 @@ * @f] * where @f$g@f$ is a user-defined constant value. */ -class UniformSteadyActiveStress : public ActiveStress { +class ActiveStressUniformSteady : public ActiveStress { public: /// Model label. static inline const std::string label = "UniformSteady"; @@ -33,7 +33,7 @@ class UniformSteadyActiveStress : public ActiveStress { /** * @brief Constructor. */ - UniformSteadyActiveStress() : ActiveStress(/* n_states = */ 0) {} + ActiveStressUniformSteady() : ActiveStress(/* n_states = */ 0) {} /** * @brief Construct an instance of model parameters. diff --git a/Code/Source/solver/active_stress_uniform_unsteady.cpp b/Code/Source/solver/ActiveStressUniformUnsteady.cpp similarity index 67% rename from Code/Source/solver/active_stress_uniform_unsteady.cpp rename to Code/Source/solver/ActiveStressUniformUnsteady.cpp index 708424085..99ed50b89 100644 --- a/Code/Source/solver/active_stress_uniform_unsteady.cpp +++ b/Code/Source/solver/ActiveStressUniformUnsteady.cpp @@ -1,33 +1,33 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_uniform_unsteady.h" +#include "ActiveStressUniformUnsteady.h" #include #include -void UniformUnsteadyActiveStress::init(const unsigned int tnNo) { +void ActiveStressUniformUnsteady::init(const unsigned int tnNo) { ActiveStress::init(tnNo); fourier_interpolation = FourierInterpolation::from_time_series_file( temporal_values_file_path, /* n_components = */ 1, ramp); } -void UniformUnsteadyActiveStress::read_model_specific_parameters( +void ActiveStressUniformUnsteady::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { ramp = params.get_bool("Ramp"); temporal_values_file_path = params.get_string("Temporal_values_file_path"); } -void UniformUnsteadyActiveStress::distribute_model_specific_parameters( +void ActiveStressUniformUnsteady::distribute_model_specific_parameters( const CmMod &cm_mod, const cmType &cm) { cm.bcast(cm_mod, &ramp); cm.bcast(cm_mod, temporal_values_file_path); } -double UniformUnsteadyActiveStress::compute_active_tension_local( +double ActiveStressUniformUnsteady::compute_active_tension_local( const Vector &state, const double fiber_stretch) const { return fourier_interpolation.value(time)[0]; } -REGISTER_ACTIVE_STRESS_MODEL("UniformUnsteady", UniformUnsteadyActiveStress); \ No newline at end of file +REGISTER_ACTIVE_STRESS_MODEL("UniformUnsteady", ActiveStressUniformUnsteady); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_uniform_unsteady.h b/Code/Source/solver/ActiveStressUniformUnsteady.h similarity index 95% rename from Code/Source/solver/active_stress_uniform_unsteady.h rename to Code/Source/solver/ActiveStressUniformUnsteady.h index f28b05165..8a00f5258 100644 --- a/Code/Source/solver/active_stress_uniform_unsteady.h +++ b/Code/Source/solver/ActiveStressUniformUnsteady.h @@ -4,9 +4,9 @@ #ifndef ACTIVE_STRESS_UNIFORM_UNSTEADY_H #define ACTIVE_STRESS_UNIFORM_UNSTEADY_H -#include "FourierInterpolation.h" +#include "ActiveStress.h" -#include "active_stress.h" +#include "FourierInterpolation.h" /** * @brief Uniform and time dependent active stress model. @@ -18,7 +18,7 @@ * @f] * where @f$g(t)@f$ is a user-defined function of time. */ -class UniformUnsteadyActiveStress : public ActiveStress { +class ActiveStressUniformUnsteady : public ActiveStress { public: /// Model label. static inline const std::string label = "UniformUnsteady"; @@ -37,7 +37,7 @@ class UniformUnsteadyActiveStress : public ActiveStress { /** * @brief Constructor. */ - UniformUnsteadyActiveStress() : ActiveStress(/* n_states = */ 0) {} + ActiveStressUniformUnsteady() : ActiveStress(/* n_states = */ 0) {} /** * @brief Construct an instance of model parameters. diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index b37e8b985..e7db6a830 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -223,18 +223,18 @@ set(CSRCS Timer.h ArtificialNeuralNetMaterial.h ArtificialNeuralNetMaterial.cpp - ionic_model.cpp - ionic_aliev_panfilov.cpp - ionic_bueno_orovio.cpp - ionic_fitzhugh_nagumo.cpp - ionic_ttp.cpp - - active_stress.cpp - active_stress_uniform_steady.cpp - active_stress_uniform_unsteady.cpp - active_stress_ode.cpp - active_stress_nash_panfilov.cpp - active_stress_regazzoni.cpp + IonicModel.cpp + IonicModelAlievPanfilov.cpp + IonicModelBuenoOrovio.cpp + IonicModelFitzHughNagumo.cpp + IonicModelTTP.cpp + + ActiveStress.cpp + ActiveStressUniformSteady.cpp + ActiveStressUniformUnsteady.cpp + ActiveStressODE.cpp + ActiveStressNashPanfilov.cpp + ActiveStressRegazzoni.cpp SPLIT.c diff --git a/Code/Source/solver/CepMod.h b/Code/Source/solver/CepMod.h index 6d280228b..950b94896 100644 --- a/Code/Source/solver/CepMod.h +++ b/Code/Source/solver/CepMod.h @@ -12,7 +12,7 @@ #define CEP_MOD_H #include "consts.h" -#include "ionic_model.h" +#include "IonicModel.h" #include "Array.h" #include "Vector.h" diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 2e9f60068..0fecefc45 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -22,7 +22,7 @@ #include "SolutionStates.h" #include "Timer.h" #include "Vector.h" -#include "active_stress.h" +#include "ActiveStress.h" #include "DebugMsg.h" @@ -790,6 +790,10 @@ class cplBCType { public: cplBCType(); + + /// @brief Index of the equation that this condition is associated with. + unsigned int equationIndex = 0; + /// @brief Is multi-domain active bool coupled = false; @@ -824,6 +828,19 @@ class cplBCType consts::CplBCType schm = consts::CplBCType::cplBC_NA; //int schm = cplBC_NA; + /// @brief Absolute floor on the flow-rate perturbation used to + /// finite-difference the coupled-BC tangent dP/dQ in + /// \c set_bc::calc_der_cpl_bc. This is a dimensional quantity, so it must + /// be set consistently with the unit system of the simulation. + double finite_difference_absolute_perturbation = 1.0e-7; + + /// @brief Flow-rate perturbation used to finite-difference the coupled-BC + /// tangent dP/dQ in \c set_bc::calc_der_cpl_bc, relative to the RMS coupled + /// flow rate. The perturbation actually applied is + /// \c max(rms(Q)*finite_difference_relative_perturbation, + /// \c finite_difference_absolute_perturbation). + double finite_difference_relative_perturbation = 1.0e-5; + /// @brief Path to the 0D code binary file std::string binPath; diff --git a/Code/Source/solver/CoupledBoundaryCondition.cpp b/Code/Source/solver/CoupledBoundaryCondition.cpp index 2d2ed353d..5785035db 100644 --- a/Code/Source/solver/CoupledBoundaryCondition.cpp +++ b/Code/Source/solver/CoupledBoundaryCondition.cpp @@ -318,12 +318,16 @@ void CoupledBoundaryCondition::compute_flowrates(ComMod& com_mod, const CmMod& c int nsd = com_mod.nsd; const auto& Yo = solutions.old.get_velocity(); const auto& Yn = solutions.current.get_velocity(); - - Qo_ = all_fun::integ(com_mod, cm_mod, *face_, Yo, 0, solutions, - std::optional(nsd - 1), false, flowrate_cfg_o_); - Qn_ = all_fun::integ(com_mod, cm_mod, *face_, Yn, 0, solutions, - std::optional(nsd - 1), false, flowrate_cfg_n_); - + const unsigned int equation_offset = + com_mod.eq[com_mod.cplBC.equationIndex].s; + + Qo_ = all_fun::integ(com_mod, cm_mod, *face_, Yo, equation_offset, + solutions, equation_offset + nsd - 1, false, + flowrate_cfg_o_, equation_offset); + Qn_ = all_fun::integ(com_mod, cm_mod, *face_, Yn, equation_offset, + solutions, equation_offset + nsd - 1, false, + flowrate_cfg_n_, equation_offset); + if (has_cap_) { const auto [Qo_cap, Qn_cap] = calculate_cap_contribution(com_mod, cm_mod, solutions, flowrate_cfg_o_, flowrate_cfg_n_); @@ -347,12 +351,17 @@ void CoupledBoundaryCondition::compute_pressures(ComMod& com_mod, const CmMod& c double area = face_->area; const auto& Yo = solutions.old.get_velocity(); const auto& Yn = solutions.current.get_velocity(); - - Po_ = all_fun::integ(com_mod, cm_mod, *face_, Yo, nsd, solutions, - std::nullopt, false, flowrate_cfg_o_) / area; - Pn_ = all_fun::integ(com_mod, cm_mod, *face_, Yn, nsd, solutions, - std::nullopt, false, flowrate_cfg_n_) / area; - + const unsigned int equation_offset = + com_mod.eq[com_mod.cplBC.equationIndex].s; + + Po_ = all_fun::integ(com_mod, cm_mod, *face_, Yo, equation_offset + nsd, + solutions, std::nullopt, false, flowrate_cfg_o_, + equation_offset) / + area; + Pn_ = all_fun::integ(com_mod, cm_mod, *face_, Yn, equation_offset + nsd, + solutions, std::nullopt, false, flowrate_cfg_n_, + equation_offset) / + area; } double CoupledBoundaryCondition::get_Qo() const @@ -583,7 +592,8 @@ void CoupledBoundaryCondition::initialize_cap(ComMod& com_mod) namespace { /// @brief Gathers cap-node mesh state on the serial rank (columns 0..n_cap-1 match \p cap_gn order). -void gather_global_mesh_state_serial(ComMod& com_mod, const SolutionStates& solutions, bool gather_Y, int nsd, int tnNo, +void gather_global_mesh_state_serial(ComMod& com_mod, const SolutionStates& solutions, bool gather_Y, + int equation_offset, int nsd, int tnNo, const Vector& cap_gn, CapGlobalMeshState& out) { const auto& Do = solutions.old.get_displacement(); @@ -610,14 +620,16 @@ void gather_global_mesh_state_serial(ComMod& com_mod, const SolutionStates& solu continue; } for (int i = 0; i < nsd; i++) { + // x is geometry (rows 0..nsd-1); Do/Dn are solution rows of the + // coupled equation, hence the equation_offset shift. out.x(i, a) = com_mod.x(i, Ac); - out.Do(i, a) = Do(i, Ac); - out.Dn(i, a) = Dn(i, Ac); + out.Do(i, a) = Do(equation_offset + i, Ac); + out.Dn(i, a) = Dn(equation_offset + i, Ac); } if (gather_Y) { for (int i = 0; i < nsd; i++) { - out.Yo(i, a) = Yo(i, Ac); - out.Yn(i, a) = Yn(i, Ac); + out.Yo(i, a) = Yo(equation_offset + i, Ac); + out.Yn(i, a) = Yn(equation_offset + i, Ac); } } break; @@ -628,7 +640,8 @@ void gather_global_mesh_state_serial(ComMod& com_mod, const SolutionStates& solu /// @brief Gathers cap-node mesh state on the MPI root (columns 0..n_cap-1 match \p cap_gn order). void gather_global_mesh_state_parallel(ComMod& com_mod, const CmMod& cm_mod, cmType& cm, const SolutionStates& solutions, - bool gather_Y, int nsd, int tnNo, int root, int nProcs, const Vector& cap_gn, + bool gather_Y, int equation_offset, int nsd, int tnNo, int root, int nProcs, + const Vector& cap_gn, const std::unordered_map& g_to_cap_col, CapGlobalMeshState& out) { @@ -664,21 +677,23 @@ void gather_global_mesh_state_parallel(ComMod& com_mod, const CmMod& cm_mod, cmT continue; } send_buf(idx++) = static_cast(g); + // x is geometry (rows 0..nsd-1); Do/Dn/Yo/Yn are solution rows of the + // coupled equation, hence the equation_offset shift. for (int i = 0; i < nsd; i++) { send_buf(idx++) = com_mod.x(i, Ac); } for (int i = 0; i < nsd; i++) { - send_buf(idx++) = Do(i, Ac); + send_buf(idx++) = Do(equation_offset + i, Ac); } for (int i = 0; i < nsd; i++) { - send_buf(idx++) = Dn(i, Ac); + send_buf(idx++) = Dn(equation_offset + i, Ac); } if (gather_Y) { for (int i = 0; i < nsd; i++) { - send_buf(idx++) = Yo(i, Ac); + send_buf(idx++) = Yo(equation_offset + i, Ac); } for (int i = 0; i < nsd; i++) { - send_buf(idx++) = Yn(i, Ac); + send_buf(idx++) = Yn(equation_offset + i, Ac); } } } @@ -760,13 +775,19 @@ void CoupledBoundaryCondition::gather_global_mesh_state(ComMod& com_mod, const C const int root = cm_mod.master; const int nProcs = cm.np(); + // The coupled equation's velocity and displacement occupy solution rows + // starting at this offset. The gather must read them there rather than + // assuming the coupled equation is the first one (offset 0). This mirrors + // the offset used by compute_flowrates / compute_pressures. + const int equation_offset = com_mod.eq[com_mod.cplBC.equationIndex].s; + if (cm.seq()) { - gather_global_mesh_state_serial(com_mod, solutions, gather_Y, nsd, tnNo, cap_mesh_global_node_ids_, - cap_global_mesh_state_); + gather_global_mesh_state_serial(com_mod, solutions, gather_Y, equation_offset, nsd, tnNo, + cap_mesh_global_node_ids_, cap_global_mesh_state_); return; } - gather_global_mesh_state_parallel(com_mod, cm_mod, cm, solutions, gather_Y, nsd, tnNo, root, nProcs, + gather_global_mesh_state_parallel(com_mod, cm_mod, cm, solutions, gather_Y, equation_offset, nsd, tnNo, root, nProcs, cap_mesh_global_node_ids_, cap_g_to_cap_col_, cap_global_mesh_state_); } diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 75b9b5447..5fd4557c2 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -92,7 +92,7 @@ bool Integrator::step() { iEqOld = cEq; auto& eq = com_mod.eq[cEq]; - if (com_mod.cplBC.coupled && cEq == 0) { + if (cEq == com_mod.cplBC.equationIndex && com_mod.cplBC.coupled) { #ifdef debug_integrator_step dmsg << "Set coupled BCs " << std::endl; #endif diff --git a/Code/Source/solver/ionic_model.cpp b/Code/Source/solver/IonicModel.cpp similarity index 99% rename from Code/Source/solver/ionic_model.cpp rename to Code/Source/solver/IonicModel.cpp index 6bcfdd0bd..1472d1772 100644 --- a/Code/Source/solver/ionic_model.cpp +++ b/Code/Source/solver/IonicModel.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "ionic_model.h" +#include "IonicModel.h" #include "ComMod.h" #include "Parameters.h" diff --git a/Code/Source/solver/ionic_model.h b/Code/Source/solver/IonicModel.h similarity index 100% rename from Code/Source/solver/ionic_model.h rename to Code/Source/solver/IonicModel.h diff --git a/Code/Source/solver/ionic_aliev_panfilov.cpp b/Code/Source/solver/IonicModelAlievPanfilov.cpp similarity index 98% rename from Code/Source/solver/ionic_aliev_panfilov.cpp rename to Code/Source/solver/IonicModelAlievPanfilov.cpp index 2d0476f41..9914ddbce 100644 --- a/Code/Source/solver/ionic_aliev_panfilov.cpp +++ b/Code/Source/solver/IonicModelAlievPanfilov.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "ionic_aliev_panfilov.h" +#include "IonicModelAlievPanfilov.h" void AlievPanfilov::read_parameters(const IonicModelParameters ¶ms) { IonicModel::read_parameters(params); diff --git a/Code/Source/solver/ionic_aliev_panfilov.h b/Code/Source/solver/IonicModelAlievPanfilov.h similarity index 99% rename from Code/Source/solver/ionic_aliev_panfilov.h rename to Code/Source/solver/IonicModelAlievPanfilov.h index 9f0c956ee..5e1edef8a 100644 --- a/Code/Source/solver/ionic_aliev_panfilov.h +++ b/Code/Source/solver/IonicModelAlievPanfilov.h @@ -4,7 +4,7 @@ #ifndef IONIC_ALIEV_PANFILOV_H #define IONIC_ALIEV_PANFILOV_H -#include "ionic_model.h" +#include "IonicModel.h" #include "Vector.h" diff --git a/Code/Source/solver/ionic_bueno_orovio.cpp b/Code/Source/solver/IonicModelBuenoOrovio.cpp similarity index 99% rename from Code/Source/solver/ionic_bueno_orovio.cpp rename to Code/Source/solver/IonicModelBuenoOrovio.cpp index fd8001dcb..88617071e 100644 --- a/Code/Source/solver/ionic_bueno_orovio.cpp +++ b/Code/Source/solver/IonicModelBuenoOrovio.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "ionic_bueno_orovio.h" +#include "IonicModelBuenoOrovio.h" void BuenoOrovio::read_parameters(const IonicModelParameters ¶ms) { IonicModel::read_parameters(params); diff --git a/Code/Source/solver/ionic_bueno_orovio.h b/Code/Source/solver/IonicModelBuenoOrovio.h similarity index 99% rename from Code/Source/solver/ionic_bueno_orovio.h rename to Code/Source/solver/IonicModelBuenoOrovio.h index fc8868f4f..e59134945 100644 --- a/Code/Source/solver/ionic_bueno_orovio.h +++ b/Code/Source/solver/IonicModelBuenoOrovio.h @@ -4,7 +4,7 @@ #ifndef IONIC_BUENO_OROVIO_H #define IONIC_BUENO_OROVIO_H -#include "ionic_model.h" +#include "IonicModel.h" #include "Vector.h" #include "utils.h" diff --git a/Code/Source/solver/ionic_fitzhugh_nagumo.cpp b/Code/Source/solver/IonicModelFitzHughNagumo.cpp similarity index 97% rename from Code/Source/solver/ionic_fitzhugh_nagumo.cpp rename to Code/Source/solver/IonicModelFitzHughNagumo.cpp index 608196377..06c6585e2 100644 --- a/Code/Source/solver/ionic_fitzhugh_nagumo.cpp +++ b/Code/Source/solver/IonicModelFitzHughNagumo.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "ionic_fitzhugh_nagumo.h" +#include "IonicModelFitzHughNagumo.h" void FitzHughNagumo::read_parameters(const IonicModelParameters ¶ms) { IonicModel::read_parameters(params); diff --git a/Code/Source/solver/ionic_fitzhugh_nagumo.h b/Code/Source/solver/IonicModelFitzHughNagumo.h similarity index 99% rename from Code/Source/solver/ionic_fitzhugh_nagumo.h rename to Code/Source/solver/IonicModelFitzHughNagumo.h index 572aa2ad4..2ea6b4c77 100644 --- a/Code/Source/solver/ionic_fitzhugh_nagumo.h +++ b/Code/Source/solver/IonicModelFitzHughNagumo.h @@ -4,7 +4,7 @@ #ifndef IONIC_FITZHUGH_NAGUMO_H #define IONIC_FITZHUGH_NAGUMO_H -#include "ionic_model.h" +#include "IonicModel.h" #include "Vector.h" #include "utils.h" diff --git a/Code/Source/solver/ionic_ttp.cpp b/Code/Source/solver/IonicModelTTP.cpp similarity index 99% rename from Code/Source/solver/ionic_ttp.cpp rename to Code/Source/solver/IonicModelTTP.cpp index c7cf58152..15e6e9979 100644 --- a/Code/Source/solver/ionic_ttp.cpp +++ b/Code/Source/solver/IonicModelTTP.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "ionic_ttp.h" +#include "IonicModelTTP.h" void TTP::read_parameters(const IonicModelParameters ¶ms) { IonicModel::read_parameters(params); @@ -456,5 +456,4 @@ Vector TTP::getf(const unsigned int zone_id, const Vector &X, return dX; } - REGISTER_IONIC_MODEL("TTP", TTP); \ No newline at end of file diff --git a/Code/Source/solver/ionic_ttp.h b/Code/Source/solver/IonicModelTTP.h similarity index 99% rename from Code/Source/solver/ionic_ttp.h rename to Code/Source/solver/IonicModelTTP.h index 1a90408b2..2602c036c 100644 --- a/Code/Source/solver/ionic_ttp.h +++ b/Code/Source/solver/IonicModelTTP.h @@ -4,7 +4,7 @@ #ifndef IONIC_TTP_H #define IONIC_TTP_H -#include "ionic_model.h" +#include "IonicModel.h" #include "Parameters.h" diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 8f9eed3d6..baa8638ee 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1208,6 +1208,14 @@ svZeroDSolverInterfaceParameters::svZeroDSolverInterfaceParameters() { set_parameter("Initial_flows", 0.0, !required, initial_flows); set_parameter("Initial_pressures", 0.0, !required, initial_pressures); + // Finite-difference perturbation used for the coupled-BC tangent dP/dQ, + // diff = max(rms(Q) * Finite_difference_relative_perturbation, + // Finite_difference_absolute_perturbation). + set_parameter("Finite_difference_absolute_perturbation", 1.0e-7, !required, + finite_difference_absolute_perturbation); + set_parameter("Finite_difference_relative_perturbation", 1.0e-5, !required, + finite_difference_relative_perturbation); + set_parameter("Configuration_file", "", required, configuration_file); set_parameter("Shared_library", "", required, shared_library); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index b18bd4661..a1d18d437 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -657,6 +657,17 @@ class CoupleGenBCParameters : public ParameterLists //---------------------------------- // svZeroDSolverInterfaceParameters //---------------------------------- +/// @brief Parameters for coupling to the svZeroDSolver (0D lumped-parameter solver). +/// +/// XML element: \code {.xml} +/// +/// implicit +/// svzerod_3Dcoupling.json +/// /path/to/libsvzero_interface +/// 1.0e-7 +/// 1.0e-5 +/// +/// \endcode // class svZeroDSolverInterfaceParameters : public ParameterLists { @@ -674,6 +685,9 @@ class svZeroDSolverInterfaceParameters : public ParameterLists Parameter initial_flows; Parameter initial_pressures; + Parameter finite_difference_absolute_perturbation; + Parameter finite_difference_relative_perturbation; + Parameter shared_library; bool value_set = false; diff --git a/Code/Source/solver/all_fun.cpp b/Code/Source/solver/all_fun.cpp index 623b77412..bdb15ab07 100644 --- a/Code/Source/solver/all_fun.cpp +++ b/Code/Source/solver/all_fun.cpp @@ -676,18 +676,10 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array& s, - const SolutionStates& solutions, bool pFlag, MechanicalConfigurationType cfg) -{ +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Vector &s, const SolutionStates &solutions, + bool pFlag, MechanicalConfigurationType cfg, + const unsigned int displacement_index) { using namespace consts; #define n_debug_integ_s #ifdef debug_integ_s @@ -811,7 +803,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co if (!isIB) { // Get normal vector in cfg configuration auto Nx = fs.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, insd, fs.eNoN, Nx, n, solutions, cfg); + n = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, cfg, displacement_index); } // Calculating the Jacobian (encodes area of face element) @@ -835,23 +827,23 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co } result = com_mod.cm.reduce(cm_mod, result); - return result; + return result; } -/// @brief This routine integrates vector field s dotted with the face normal n -/// over the face lFa. For example, if s contains the velocity at each node on -/// the face, this function computed the velocity flux through the face. -/// -/// Reproduces 'FUNCTION IntegV(lFa, s)' -/// -/// @param lFa face type, representing a face on the computational mesh -/// @param s an array containing a vector value for each node in the mesh -/// @param pFlag flag for using Taylor-Hood function space for pressure -/// @param cfg denotes which configuration (reference/timestep 0, old/timestep n, or new/timestep n+1). Default reference. -// -double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, - const Array& s, const SolutionStates& solutions, MechanicalConfigurationType cfg) -{ +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Vector &s, const SolutionStates &solutions, + bool pFlag) { + // The displacement index is not used in the reference configuration. + constexpr unsigned int unused_displacement_index = 0; + return integ(com_mod, cm_mod, lFa, s, solutions, pFlag, + consts::MechanicalConfigurationType::reference, + unused_displacement_index); +} + +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const SolutionStates &solutions, + MechanicalConfigurationType cfg, + const unsigned int displacement_index) { using namespace consts; #define n_debug_integ_V @@ -931,7 +923,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, if (!isIB) { // Get normal vector in cfg configuration auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, solutions, cfg); + n = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, cfg, displacement_index); //CALL GNNB(lFa, e, g, nsd-1, lFa.eNoN, lFa.Nx(:,:,g), n) } else { //CALL GNNIB(lFa, e, g, n) @@ -963,29 +955,23 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, result = cm.reduce(cm_mod, result); - return result; + return result; } -/// @brief This routine integrate s(l:u,:) over the surface faId, where s is an -/// array of scalars or an array of nsd-vectors. This routine calls overloaded -/// functions to integrate scalars, if s is scalar (i.e. l=u), or vectors if s -/// is vector (i.e. l& s, const int l, const SolutionStates& solutions, std::optional uo, bool THflag, MechanicalConfigurationType cfg) -{ +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const SolutionStates &solutions) { + // The displacement index is not used in the reference configuration. + constexpr unsigned int unused_displacement_index = 0; + return integ(com_mod, cm_mod, lFa, s, solutions, + consts::MechanicalConfigurationType::reference, + unused_displacement_index); +} + +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const int l, + const SolutionStates &solutions, std::optional uo, + bool THflag, MechanicalConfigurationType cfg, + const unsigned int displacement_index) { using namespace consts; #define n_debug_integ_g @@ -1040,21 +1026,33 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, vec(n,a) = s(i,a); } } - result = integ(com_mod, cm_mod, lFa, vec, solutions, cfg); - // If s scalar, integrate as scalar + result = + integ(com_mod, cm_mod, lFa, vec, solutions, cfg, displacement_index); + // If s scalar, integrate as scalar } else if (l == u) { Vector sclr(nNo); for (int a = 0; a < nNo; a++) { sclr(a) = s(l,a); } - result = integ(com_mod, cm_mod, lFa, sclr, solutions, flag, cfg); + result = integ(com_mod, cm_mod, lFa, sclr, solutions, flag, cfg, + displacement_index); } else { throw std::runtime_error("Unexpected dof in integ"); } - return result; + return result; } +double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const int l, + const SolutionStates &solutions, std::optional uo, + bool THflag) { + // The displacement index is not used in the reference configuration. + constexpr unsigned int unused_displacement_index = 0; + return integ(com_mod, cm_mod, lFa, s, l, solutions, uo, THflag, + consts::MechanicalConfigurationType::reference, + unused_displacement_index); +} bool is_domain(const ComMod& com_mod, const eqType& eq, const int node, const consts::EquationType phys) { diff --git a/Code/Source/solver/all_fun.h b/Code/Source/solver/all_fun.h index 2aa588cfa..3cf5bb8a8 100644 --- a/Code/Source/solver/all_fun.h +++ b/Code/Source/solver/all_fun.h @@ -35,13 +35,162 @@ namespace all_fun { double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array& s, int l, int u, const SolutionStates& solutions, bool pFlag=false); - double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Vector& s, - const SolutionStates& solutions, bool pFlag=false, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); - - double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& s, - const int l, const SolutionStates& solutions, std::optional uo=std::nullopt, bool THflag=false, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); - - double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& s, const SolutionStates& solutions, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); + /** + * @brief Integrate a scalar field over a boundary face. + * + * Reproduces 'FUNCTION IntegS(lFa, s, pflag)'. The scalar field s is + * integrated over the face, i.e. this computes + * \f[ + * \int_{\Gamma} s \, d\Gamma, + * \f] + * where \f$\Gamma\f$ is the boundary face. + * + * For simulations involving structural displacement, this function allows + * computing the integral in any of the following configurations: + * - reference configuration (the mesh is not displaced); + * - current configuration (the mesh is displaced by the current displacement + * field); + * - old configuration (the mesh is displaced by the displacement field from + * previous time step). + * + * @param[in] com_mod The common module. + * @param[in] cm_mod The communication module containing MPI data. + * @param[in] lFa The boundary face over which the integral is computed. + * @param[in] s The scalar value at each node of the mesh. + * @param[in] solutions The solution states that the displacement fields are + * extracted from. + * @param[in] pFlag Whether to use the Taylor-Hood function space for the + * pressure field. + * @param[in] cfg The configuration in which the integral is computed + * (reference, old or current). + * @param[in] displacement_index The index of the displacement field in the + * solution arrays. This should correspond to the start index of the + * equation that solves for the displacement. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Vector &s, const SolutionStates &solutions, + bool pFlag, consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index); + + /** + * @brief Integrate a scalar field over a boundary face. + * + * This is the overload to use in the general case. The other overload adds + * the ability to integrate over a displaced configuration, which is only + * relevant for simulations involving structural displacement. See it for the + * meaning of the arguments. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Vector &s, const SolutionStates &solutions, + bool pFlag); + + /** + * @brief Integrate one or more components of a field over a boundary face. + * + * Reproduces 'FUNCTION IntegG(lFa, s, l, u, THflag)'. Rows l to u of s are + * integrated over the face. When they span the spatial dimensions the field + * is dotted with the outward surface normal (i.e. a flux is computed), + * \f[ + * \int_{\Gamma} \sum_{i=l}^{u} s_i \, n_i \, d\Gamma, + * \f] + * otherwise the single row l is integrated as a scalar, + * \f[ + * \int_{\Gamma} s_l \, d\Gamma, + * \f] + * where \f$\Gamma\f$ is the boundary face and \f$\mathbf{n}\f$ its outward + * unit normal. + * + * For simulations involving structural displacement, this function allows + * computing the integral in any of the following configurations: + * - reference configuration (the mesh is not displaced); + * - current configuration (the mesh is displaced by the current displacement + * field); + * - old configuration (the mesh is displaced by the displacement field from + * previous time step). + * + * @param[in] com_mod The common module. + * @param[in] cm_mod The communication module containing MPI data. + * @param[in] lFa The boundary face over which the integral is computed. + * @param[in] s The field value at each node of the mesh. + * @param[in] l The first row of s to integrate. + * @param[in] solutions The solution states that the displacement fields are + * extracted from. + * @param[in] uo The last row of s to integrate. Defaults to l, i.e. a single + * component. + * @param[in] THflag Whether to use the Taylor-Hood function space for the + * pressure field. + * @param[in] cfg The configuration in which the integral is computed + * (reference, old or current). + * @param[in] displacement_index The index of the displacement field in the + * solution arrays. This should correspond to the start index of the + * equation that solves for the displacement. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const int l, + const SolutionStates &solutions, std::optional uo, + bool THflag, consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index); + + /** + * @brief Integrate one or more components of a field over a boundary face. + * + * This is the overload to use in the general case. The other overload adds + * the ability to integrate over a displaced configuration, which is only + * relevant for simulations involving structural displacement. See it for the + * meaning of the arguments. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const int l, + const SolutionStates &solutions, std::optional uo, + bool THflag); + + /** + * @brief Integrate the flux of a vector field over a boundary face. + * + * Reproduces 'FUNCTION IntegV(lFa, s)'. The vector field s (one component per + * spatial dimension at each node) is dotted with the outward surface normal + * and integrated over the face, i.e. this computes + * \f[ + * \int_{\Gamma} \mathbf{s} \cdot \mathbf{n} \, d\Gamma, + * \f] + * where \f$\Gamma\f$ is the boundary face and \f$\mathbf{n}\f$ its outward + * unit normal. + * + * For simulations involving structural displacement, this function allows + * computing the integral in any of the following configurations: + * - reference configuration (the mesh is not displaced); + * - current configuration (the mesh is displaced by the current displacement + * field); + * - old configuration (the mesh is displaced by the displacement field from + * previous time step). + * + * @param[in] com_mod The common module. + * @param[in] cm_mod The communication module containing MPI data. + * @param[in] lFa The boundary face over which the integral is computed. + * @param[in] s The vector value at each node of the mesh. + * @param[in] solutions The solution states that the displacement fields are + * extracted from. + * @param[in] cfg The configuration in which the integral is computed + * (reference, old or current). + * @param[in] displacement_index The index of the displacement field in the + * solution arrays. This should correspond to the start index of the + * equation that solves for the displacement. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const SolutionStates &solutions, + consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index); + + /** + * @brief Integrate the flux of a vector field over a boundary face. + * + * This is the overload to use in the general case. The other overload adds + * the ability to integrate over a displaced configuration, which is only + * relevant for simulations involving structural displacement. See it for the + * meaning of the arguments. + */ + double integ(const ComMod &com_mod, const CmMod &cm_mod, const faceType &lFa, + const Array &s, const SolutionStates &solutions); bool is_domain(const ComMod& com_mod, const eqType& eq, const int node, const consts::EquationType phys); diff --git a/Code/Source/solver/baf_ini.cpp b/Code/Source/solver/baf_ini.cpp index e83dd8a64..4f3f05b28 100644 --- a/Code/Source/solver/baf_ini.cpp +++ b/Code/Source/solver/baf_ini.cpp @@ -98,7 +98,7 @@ void baf_ini(Simulation* simulation, SolutionStates& solutions) // cplBC faces are initialized here // - int iEq = 0; + int iEq = com_mod.cplBC.equationIndex; com_mod.cplBC.fa.resize(com_mod.cplBC.nFa); com_mod.cplBC.xn.resize(com_mod.cplBC.nX); @@ -163,7 +163,7 @@ void baf_ini(Simulation* simulation, SolutionStates& solutions) } if (com_mod.cplBC.useSvZeroD) { - svZeroD::init_svZeroD(com_mod, cm_mod); + svZeroD::init_svZeroD(com_mod, cm_mod, simulation->get_chnl_mod().appPath); } if (com_mod.cplBC.useSvOneD) { @@ -342,7 +342,9 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l } else if (btest(lBc.bType, iBC_para)) { Vector center(3); for (int i = 0; i < nsd; i++) { - center(i) = all_fun::integ(com_mod, cm_mod, lFa, com_mod.x, i, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference) / lFa.area; + center(i) = all_fun::integ(com_mod, cm_mod, lFa, com_mod.x, i, solutions, + std::nullopt, false) / + lFa.area; } // gNodes is one if a node located on the boundary (beside iFa) @@ -450,8 +452,8 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l // Normalizing the profile for flux // double tmp = 1.0; - if (btest(lBc.bType, enum_int(BoundaryConditionType::bType_flx))) { - tmp = all_fun::integ(com_mod, cm_mod, lFa, s, solutions, false, consts::MechanicalConfigurationType::reference); + if (btest(lBc.bType, enum_int(BoundaryConditionType::bType_flx))) { + tmp = all_fun::integ(com_mod, cm_mod, lFa, s, solutions, false); if (is_zero(tmp)) { tmp = 1.0; throw std::runtime_error("Face '" + lFa.name + "' used for a BC has no non-zero node."); @@ -493,7 +495,7 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa, const Solution // Vector sA(com_mod.tnNo); sA = 1.0; - double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, consts::MechanicalConfigurationType::reference); + double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false); #ifdef debug_face_ini dmsg << "Face '" << lFa.name << "' area: " << area; #endif @@ -525,7 +527,6 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa, const Solution // Compute integral of normal vector over surface element if (!flag) { - Vector nV(nsd); for (int e = 0; e < lFa.nEl; e++) { if (lFa.eType == ElementType::NRB) { @@ -535,7 +536,7 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa, const Solution for (int g = 0; g < lFa.nG; g++) { auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + const Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); @@ -797,9 +798,8 @@ void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceTyp // CALL NRBNNXB(msh(iM),lFa,e) } for (int g = 0; g < lFa.nG; g++) { - Vector n(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, solutions, consts::MechanicalConfigurationType::reference); + const Vector n = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); diff --git a/Code/Source/solver/cmm.cpp b/Code/Source/solver/cmm.cpp index 3d78b5b96..785a3429b 100644 --- a/Code/Source/solver/cmm.cpp +++ b/Code/Source/solver/cmm.cpp @@ -280,9 +280,8 @@ void cmm_b(ComMod& com_mod, const faceType& lFa, const int e, const Array nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, 3, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 00986613e..111148266 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -18,7 +18,7 @@ #include #include -#include "ionic_model.h" +#include "IonicModel.h" extern "C" { @@ -551,8 +551,10 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &cplBC.useGenBC); cm.bcast(cm_mod, &cplBC.useSvZeroD); cm.bcast(cm_mod, &cplBC.useSvOneD); + cm.bcast(cm_mod, &cplBC.finite_difference_absolute_perturbation); + cm.bcast(cm_mod, &cplBC.finite_difference_relative_perturbation); - if (cplBC.useGenBC) { + if (cplBC.useGenBC) { if (cm.slv(cm_mod)) { cplBC.nX = 0; cplBC.xo.resize(cplBC.nX); diff --git a/Code/Source/solver/eq_assem.cpp b/Code/Source/solver/eq_assem.cpp index 9cea03b30..463a178ca 100644 --- a/Code/Source/solver/eq_assem.cpp +++ b/Code/Source/solver/eq_assem.cpp @@ -78,9 +78,8 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& } for (int g = 0; g < lFa.nG; g++) { - Vector nV(nsd); auto Nx = lFa.Nx.rslice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; @@ -251,9 +250,8 @@ void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const } // Get surface normal vector - Vector nV(nsd); auto Nx_g = lFa.Nx.rslice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoNb, Nx_g, nV, solutions, consts::MechanicalConfigurationType::reference); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx_g, solutions); Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g)*Jac; @@ -325,13 +323,13 @@ void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const S // CALL NRBNNXB(msh(iM),lFa,e) } for (int g = 0; g < lFa.nG; g++) { - Vector n(nsd); auto Nx = lFa.Nx.rslice(g); auto cfg = MechanicalConfigurationType::new_timestep; - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, solutions, cfg); - // + const Vector n = nn::gnnb(com_mod, lFa, e, g, Nx, solutions, cfg, + com_mod.eq[com_mod.cEq].s); + for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); for (int i = 0; i < nsd; i++) { diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 7baf4f330..8decef715 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -3,6 +3,8 @@ // The code here replicates the Fortran code in DISTRIBUTE.f. +#include "Core/Exception.h" + #include "initialize.h" #include "distribute.h" @@ -417,11 +419,24 @@ void initialize(Simulation* simulation, Vector& timeP) nFacesLS = nFacesLS + 1; } - for (auto& bc : com_mod.eq[0].bc) { - // Check for coupled faces (Dir, Neu via cplBC) or Coupled BCs - if (bc.cplBCptr != -1 || utils::btest(bc.bType, static_cast(consts::BoundaryConditionType::bType_Coupled))) { - com_mod.cplBC.coupled = true; - break; + // Check for coupled faces (Dir, Neu via cplBC) or Coupled BCs + for (unsigned int i = 0; i < com_mod.eq.size(); ++i) { + for (auto &bc : com_mod.eq[i].bc) { + + if (bc.cplBCptr != -1 || + utils::btest( + bc.bType, + static_cast(consts::BoundaryConditionType::bType_Coupled))) { + svmp::throw_if( + com_mod.cplBC.coupled && com_mod.cplBC.equationIndex != i, + "Coupled boundary conditions can only be assigned in one equation, " + "but they were assigned in equations " + + std::to_string(com_mod.cplBC.equationIndex) + " and " + + std::to_string(i) + "."); + + com_mod.cplBC.equationIndex = i; + com_mod.cplBC.coupled = true; + } } } diff --git a/Code/Source/solver/nn.cpp b/Code/Source/solver/nn.cpp index 78c6e2207..d3e57bb12 100644 --- a/Code/Source/solver/nn.cpp +++ b/Code/Source/solver/nn.cpp @@ -899,23 +899,22 @@ void gnn(const int eNoN, const int nsd, const int insd, Array& Nxi, Arra } } -/// @brief This routine returns a surface normal vector at element "e" and Gauss point -/// 'g' of face 'lFa' that is the normal weighted by Jac, i.e. -/// Jac = norm(n), the Jacobian of the mapping from parent surface element to -/// reference/old/new configuration. -/// -/// cfg denotes which configuration (reference/timestep 0, old/timestep n, or new/timestep n+1). Default reference -/// -/// Reproduce Fortran 'GNNB'. -// -void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, const int nsd, const int insd, - const int eNoNb, const Array& Nx, Vector& n, const SolutionStates& solutions, MechanicalConfigurationType cfg) -{ +Vector gnnb(const ComMod &com_mod, const faceType &lFa, const int e, + const int g, const Array &Nx, + const SolutionStates &solutions, + consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index) { // Local aliases for displacement arrays const auto& Dn = solutions.current.get_displacement(); const auto& Do = solutions.old.get_displacement(); auto& cm = com_mod.cm; + const int nsd = com_mod.nsd; + const int insd = Nx.nrows(); + const int eNoNb = Nx.ncols(); + + Vector n(nsd); + #define n_debug_gnnb #ifdef debug_gnnb DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -998,6 +997,10 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, if (com_mod.mvMsh) { for (int i = 0; i < lX.nrows(); i++) { // Add mesh displacement + // Notice that this assumes that the mesh displacement is stored + // starting at the nsd+1 index of the solution array. This is enforced + // in read_files, by throwing an exception if the equations are not + // ordered correctly. lX(i,a) = lX(i,a) + Do(i+nsd+1,Ac); } } @@ -1009,13 +1012,13 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, case MechanicalConfigurationType::old_timestep: for (int i = 0; i < lX.nrows(); i++) { // Add displacement at timestep n - lX(i,a) = lX(i,a) + Do(i,Ac); + lX(i, a) = lX(i, a) + Do(displacement_index + i, Ac); } break; case MechanicalConfigurationType::new_timestep: for (int i = 0; i < lX.nrows(); i++) { // Add displacement at timestep n+1 - lX(i,a) = lX(i,a) + Dn(i,Ac); + lX(i, a) = lX(i, a) + Dn(displacement_index + i, Ac); } break; default: @@ -1036,10 +1039,15 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, // Compute adjoining mesh element normal // - Array xXi(nsd,nsd-1); + // Note that here we use the shell element's shape function derivatives + // (msh.Nx), so the relevant intrinsic dimension is that of the shell + // element (nsd - 1, a surface), not that of the boundary edge lFa + // (which is what the local 'insd' holds). + const int msh_insd = nsd - 1; + Array xXi(nsd, msh_insd); for (int a = 0; a < eNoN; a++) { - for (int i = 0; i < insd; i++) { + for (int i = 0; i < msh_insd; i++) { for (int j = 0; j < nsd; j++) { xXi(j,i) = xXi(j,i) + lX(j,a)*msh.Nx(i,a,g); } @@ -1083,7 +1091,7 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, n = -n; } - return; + return n; } else { @@ -1099,24 +1107,36 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, } n = utils::cross(xXi); - } - // Changing the sign if neccessary. 'a' locates on the face and 'b' - // in the interior of the element. v points outward along ba - // - a = ptr(0); - int b = ptr(lFa.eNoN); - Vector v(nsd); + // Changing the sign if neccessary. 'a' locates on the face and 'b' + // in the interior of the element. v points outward along ba + // + a = ptr(0); + int b = ptr(lFa.eNoN); + Vector v(nsd); - for (int i = 0; i < nsd; i++) { - v(i) = lX(i,a) - lX(i,b); - } + for (int i = 0; i < nsd; i++) { + v(i) = lX(i, a) - lX(i, b); + } + + if (n * v < 0.0) { + n = -n; + } - if (n * v < 0.0) { - n = -n; + return n; } } +Vector gnnb(const ComMod &com_mod, const faceType &lFa, const int e, + const int g, const Array &Nx, + const SolutionStates &solutions) { + // The displacement index is not used in the reference configuration. + constexpr unsigned int unused_displacement_index = 0; + return gnnb(com_mod, lFa, e, g, Nx, solutions, + consts::MechanicalConfigurationType::reference, + unused_displacement_index); +} + /// @brief Compute shell kinematics: normal vector, covariant & contravariant basis vectors /// /// Replicates 'SUBROUTINE GNNS(eNoN, Nxi, xl, nV, gCov, gCnv)' defined in NN.f. diff --git a/Code/Source/solver/nn.h b/Code/Source/solver/nn.h index 0bd862dd8..c8c8910f3 100644 --- a/Code/Source/solver/nn.h +++ b/Code/Source/solver/nn.h @@ -38,8 +38,71 @@ namespace nn { void gnn(const int eNoN, const int nsd, const int insd, Array& Nxi, Array& x, Array& Nx, double& Jac, Array& ks); - void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g, const int nsd, const int insd, - const int eNoNb, const Array& Nx, Vector& n, const SolutionStates& solutions, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); + /** + * @brief Return the area-weighted surface normal at a given element and Gauss + * point. + * + * Returns the outward normal at element 'e', Gauss point 'g' of face 'lFa', + * weighted by the surface Jacobian: Jac = norm(n) is the Jacobian of the + * mapping from the parent surface element to the reference/old/new + * configuration. + * + * For simulations involving structural displacement, this function allows + * computing the normal vector in any of the following configurations: + * - reference configuration (the mesh is not displaced); + * - current configuration (the mesh is displaced by the current displacement + * field); + * - old configuration (the mesh is displaced by the displacement field from + * previous time step). + * + * @param[in] com_mod The common module. + * @param[in] lFa The boundary face for which the normal vector is computed. + * @param[in] e The face-local index of the element for which the normal + * vector is computed. + * @param[in] g The Gauss point index for which the normal vector is computed. + * @param[in] Nx The shape function derivatives at the Gauss point. Its shape + * (insd x eNoNb) determines the surface's intrinsic dimension and the + * number of nodes per face element. + * @param[in] solutions The solution states that the displacement fields are + * extracted from. + * @param[in] cfg The configuration in which the normal vector is computed + * (reference, old or current). + * @param[in] displacement_index The index of the displacement field in the + * solution arrays. This should correspond to the start index of the + * equation that solves for the displacement. + * @return The area-weighted outward normal vector. + */ + Vector gnnb(const ComMod &com_mod, const faceType &lFa, const int e, + const int g, const Array &Nx, + const SolutionStates &solutions, + consts::MechanicalConfigurationType cfg, + const unsigned int displacement_index); + + /** + * @brief Return the area-weighted surface normal at a given element and Gauss + * point. + * + * Returns the outward normal at element 'e', Gauss point 'g' of face 'lFa', + * weighted by the surface Jacobian: Jac = norm(n) is the Jacobian of the + * mapping from the parent surface element to the mesh. + * + * This is the overload to use in the general case. The other overload adds + * the ability to compute the normal in a displaced configuration, which is + * only relevant for simulations involving structural displacement. + * + * @param[in] com_mod The common module. + * @param[in] lFa The boundary face for which the normal vector is computed. + * @param[in] e The face-local index of the element for which the normal + * vector is computed. + * @param[in] g The Gauss point index for which the normal vector is computed. + * @param[in] Nx The shape function derivatives at the Gauss point. + * @param[in] solutions The solution states that the displacement fields are + * extracted from. + * @return The area-weighted outward normal vector. + */ + Vector gnnb(const ComMod &com_mod, const faceType &lFa, const int e, + const int g, const Array &Nx, + const SolutionStates &solutions); void gnns(const int nsd, const int eNoN, const Array& Nxi, Array& xl, Vector& nV, Array& gCov, Array& gCnv); diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index befd3b9e4..08899995f 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -7,10 +7,10 @@ #include "Core/Exception.h" #include "FE/Common/FEException.h" -#include "active_stress.h" +#include "ActiveStress.h" #include "all_fun.h" #include "consts.h" -#include "ionic_model.h" +#include "IonicModel.h" #include "read_msh.h" #include "vtk_xml.h" @@ -1669,6 +1669,12 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq) cplBC.useSvZeroD = true; cplbc_type_str = eq_params->svzerodsolver_interface_parameters.coupling_type.value(); cplBC.svzerod_solver_interface.set_data(eq_params->svzerodsolver_interface_parameters); + cplBC.finite_difference_absolute_perturbation = + eq_params->svzerodsolver_interface_parameters + .finite_difference_absolute_perturbation.value(); + cplBC.finite_difference_relative_perturbation = + eq_params->svzerodsolver_interface_parameters + .finite_difference_relative_perturbation.value(); } if (eq_params->svonedsolver_interface_parameters.defined()) { diff --git a/Code/Source/solver/ris.cpp b/Code/Source/solver/ris.cpp index cfa5b81d8..896fb8bd6 100644 --- a/Code/Source/solver/ris.cpp +++ b/Code/Source/solver/ris.cpp @@ -59,7 +59,10 @@ void ris_meanq(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) int iM = RIS.lst(i,0,iProj); int iFa = RIS.lst(i,1,iProj); double tmp = msh[iM].fa[iFa].area; - RIS.meanP(iProj,i) = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference)/tmp; + RIS.meanP(iProj, i) = + all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, + std::nullopt, false) / + tmp; } } @@ -77,7 +80,8 @@ void ris_meanq(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) } int iM = RIS.lst(0,0,iProj); int iFa = RIS.lst(0,1,iProj); - RIS.meanFl(iProj) = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference); + RIS.meanFl(iProj) = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, + 0, solutions, m - 1, false); if (cm.mas(cm_mod)) { std::cout << "For RIS projection: " << iProj << std::endl; @@ -480,8 +484,10 @@ void ris0d_status(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solution sA = 1.0; lFa = msh[iM].fa[iFa]; // such update may be not correct - tmp_new = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, consts::MechanicalConfigurationType::reference); - meanP = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference)/tmp_new; + tmp_new = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false); + meanP = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, + m - 1, false) / + tmp_new; // For the velocity m = nsd; @@ -496,7 +502,8 @@ void ris0d_status(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solution } } - meanFl = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference); + meanFl = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, + solutions, m - 1, false); std::cout << "The average pressure is: " << meanP << std::endl; std::cout << "The pressure from 0D is: " << eq[cEq].bc[iBc].g << std::endl; diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 65f744727..30086dbb4 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -13,11 +13,12 @@ #include "lhsa.h" #include "mat_fun.h" #include "nn.h" +#include "svOneD_interface.h" +#include "svZeroD_interface.h" #include "ustruct.h" #include "utils.h" +#include #include -#include "svZeroD_interface.h" -#include "svOneD_interface.h" namespace set_bc { @@ -44,11 +45,9 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& dmsg.banner(); #endif - const int iEq = 0; - // NOTE: For coupling with svZeroDPlus, absTol needs to be > 1e-8 to be compatible with the default convergence tolerance of svZeroDPlus (1e-8) - // If this is not true, the finite difference computation of bc.r below results in zero because the perturbation is below the svZeroDPlus tolerance - const double absTol = 1.0e-7; - const double relTol = 1.0e-5; + const int iEq = com_mod.cplBC.equationIndex; + const double absTol = com_mod.cplBC.finite_difference_absolute_perturbation; + const double relTol = com_mod.cplBC.finite_difference_relative_perturbation; int nsd = com_mod.nsd; auto& eq = com_mod.eq[iEq]; @@ -144,8 +143,16 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& else { throw std::runtime_error("[calc_der_cpl_bc] Invalid physics type for 0D coupling"); } - cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, fa, Yo, 0, solutions, nsd-1, false, cfg_o); - cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, fa, Yn, 0, solutions, nsd-1, false, cfg_n); + + const unsigned int equation_offset = + com_mod.eq[com_mod.cplBC.equationIndex].s; + cplBC.fa[ptr].Qo = all_fun::integ( + com_mod, cm_mod, fa, Yo, equation_offset, solutions, + equation_offset + nsd - 1, false, cfg_o, equation_offset); + cplBC.fa[ptr].Qn = all_fun::integ( + com_mod, cm_mod, fa, Yn, equation_offset, solutions, + equation_offset + nsd - 1, false, cfg_n, equation_offset); + cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; #ifdef debug_calc_der_cpl_bc @@ -158,8 +165,13 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& // Compute avg pressures at 3D Dirichlet boundaries at timesteps n and n+1 else if (utils::btest(bc.bType, iBC_Dir)) { double area = fa.area; - cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, solutions, std::nullopt, false, MechanicalConfigurationType::reference) / area; - cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, fa, Yn, nsd, solutions, std::nullopt, false, MechanicalConfigurationType::reference) / area; + cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, + solutions, std::nullopt, false) / + area; + cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, fa, Yn, nsd, + solutions, std::nullopt, false) / + area; + cplBC.fa[ptr].Qo = 0.0; cplBC.fa[ptr].Qn = 0.0; #ifdef debug_calc_der_cpl_bc @@ -616,7 +628,7 @@ void rcr_init(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& soluti using namespace consts; - const int iEq = 0; + const int iEq = com_mod.cplBC.equationIndex; int nsd = com_mod.nsd; auto& eq = com_mod.eq[iEq]; auto& cplBC = com_mod.cplBC; @@ -635,8 +647,11 @@ void rcr_init(ComMod& com_mod, const CmMod& cm_mod, const SolutionStates& soluti if (cplBC.initRCR) { auto& fa = com_mod.msh[iM].fa[iFa]; double area = fa.area; - double Qo = all_fun::integ(com_mod, cm_mod, fa, Yo, 0, solutions, nsd-1, false, MechanicalConfigurationType::reference); - double Po = all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, solutions, std::nullopt, false, MechanicalConfigurationType::reference) / area; + double Qo = all_fun::integ(com_mod, cm_mod, fa, Yo, 0, solutions, + nsd - 1, false); + double Po = all_fun::integ(com_mod, cm_mod, fa, Yo, nsd, solutions, + std::nullopt, false) / + area; cplBC.xo[ptr] = Po - (Qo * cplBC.fa[ptr].RCR.Rp); } else { cplBC.xo[ptr] = cplBC.fa[ptr].RCR.Xo; @@ -743,15 +758,13 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) const auto& Yo = solutions.old.get_velocity(); const auto& Do = solutions.old.get_displacement(); - static double absTol = 1.E-8, relTol = 1.E-5; - using namespace consts; const int nsd = com_mod.nsd; const int tnNo = com_mod.tnNo; auto& cplBC = com_mod.cplBC; // Yo, Ao, Do now passed as parameters - const int iEq = 0; + const int iEq = com_mod.cplBC.equationIndex; auto& eq = com_mod.eq[iEq]; // Determine current physics @@ -780,7 +793,7 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) faceType& lFa = com_mod.msh[iM].fa[iFa]; Vector sA(com_mod.tnNo); sA = 1.0; - double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, consts::MechanicalConfigurationType::reference); + double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false); baf_ini_ns::bc_ini(com_mod, cm_mod, eq.bc[iBc], lFa, solutions); int ptr = bc.cplBCptr; @@ -825,16 +838,36 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, const SolutionStates& solutions) else { throw std::runtime_error("[set_bc_cpl] Invalid physics type for 0D coupling"); } - - cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, 0, solutions, nsd-1, false, cfg_o); - cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, 0, solutions, nsd-1, false, cfg_n); + + const unsigned int equation_offset = + com_mod.eq[com_mod.cplBC.equationIndex].s; + cplBC.fa[ptr].Qo = all_fun::integ( + com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, equation_offset, + solutions, equation_offset + nsd - 1, false, cfg_o, + equation_offset); + cplBC.fa[ptr].Qn = all_fun::integ( + com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, equation_offset, + solutions, equation_offset + nsd - 1, false, cfg_n, + equation_offset); + cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; } // Compute avg pressures at 3D Dirichlet boundaries at timesteps n and n+1 else if (utils::btest(bc.bType,iBC_Dir)) { - cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, nsd, solutions, std::nullopt, false, MechanicalConfigurationType::reference) / area; - cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, nsd, solutions, std::nullopt, false, MechanicalConfigurationType::reference) / area; + const unsigned int equation_offset = + com_mod.eq[com_mod.cplBC.equationIndex].s; + cplBC.fa[ptr].Po = + all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, + equation_offset + nsd, solutions, std::nullopt, + false) / + area; + cplBC.fa[ptr].Pn = + all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, + equation_offset + nsd, solutions, std::nullopt, + false) / + area; + cplBC.fa[ptr].Qo = 0.0; cplBC.fa[ptr].Qn = 0.0; } @@ -1394,9 +1427,8 @@ void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const // Gauss integration 1 // for (int g = 0; g < lFa.nG; g++) { - Vector nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoNb, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g) * Jac; @@ -1547,36 +1579,37 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const //h(0) = lBc.g; double Q_3D = all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, - eq.s+nsd-1, false, - consts::MechanicalConfigurationType::reference); - - - h(0) = lBc.g; - //h(0) = lBc.g - lBc.r * std::abs(Q_3D); - - // Backflow kinetic energy correction: when backflow is detected - // (Q < 0), subtract the face-averaged dynamic pressure to further - // reduce the applied traction and damp the incoming flow. - if (Q_3D < 0.0) { - 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); - double beta = eq.dmn[cDmn_local].prop.at( - consts::PhysicalProperyType::backflow_stab); - double A = lFa.area; - if (A > 0.0) { - double u_n = Q_3D / A; // face-averaged normal velocity (< 0) - h(0) -= 0.5 * beta * rho * u_n * u_n; - } - } - } else if (utils::btest(lBc.bType,iBC_res)) { - h(0) = lBc.r * all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, eq.s+nsd-1, false, consts::MechanicalConfigurationType::reference); - - } else if (utils::btest(lBc.bType,iBC_std)) { + eq.s + nsd - 1, false); + + h(0) = lBc.g; + // h(0) = lBc.g - lBc.r * std::abs(Q_3D); + + // Backflow kinetic energy correction: when backflow is detected + // (Q < 0), subtract the face-averaged dynamic pressure to further + // reduce the applied traction and damp the incoming flow. + if (Q_3D < 0.0) { + 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); + double beta = eq.dmn[cDmn_local].prop.at( + consts::PhysicalProperyType::backflow_stab); + double A = lFa.area; + if (A > 0.0) { + double u_n = Q_3D / A; // face-averaged normal velocity (< 0) + h(0) -= 0.5 * beta * rho * u_n * u_n; + } + } + + } else if (utils::btest(lBc.bType, iBC_res)) { + h(0) = lBc.r * all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, solutions, + eq.s + nsd - 1, false); + + } else if (utils::btest(lBc.bType, iBC_std)) { h(0) = lBc.g; - } else if (utils::btest(lBc.bType,iBC_ustd)) { + } else if (utils::btest(lBc.bType, iBC_ustd)) { h = lBc.gt.value(com_mod.time); } else { @@ -1695,9 +1728,8 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const RobinBoundaryCondit lKd = 0.0; for (int g = 0; g < lFa.nG; g++) { - Vector nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); double Jac = utils::norm(nV); nV = nV / Jac; double w = lFa.w(g) * Jac; @@ -1977,9 +2009,8 @@ void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, cons lR = 0.0; for (int g = 0; g < lFa.nG; g++) { - Vector nV(nsd); auto Nx = lFa.Nx.slice(g); - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); + const Vector nV = nn::gnnb(com_mod, lFa, e, g, Nx, solutions); double Jac = utils::norm(nV); double w = lFa.w(g)*Jac; N = lFa.N.col(g); diff --git a/Code/Source/solver/svZeroD_interface.cpp b/Code/Source/solver/svZeroD_interface.cpp index 2a3c8da8f..ad791f5c7 100644 --- a/Code/Source/solver/svZeroD_interface.cpp +++ b/Code/Source/solver/svZeroD_interface.cpp @@ -58,6 +58,11 @@ static int numCoupledSrfs; static bool writeSvZeroD = true; static double svZeroDTime = 0.0; +/// Directory the svZeroD output files (svZeroD_data, Q_svZeroD, P_svZeroD) are +/// written to, set in init_svZeroD() to the simulation results directory +/// (chnl_mod.appPath). Includes the trailing separator. +static std::string svZeroD_output_dir = ""; + int num_output_steps; int system_size; int model_id; @@ -148,7 +153,7 @@ void write_svZeroD_solution(const double* lpn_time, std::vector& lpn_sol std::vector variable_names; variable_names = interface->variable_names_; std::ofstream out_file; - out_file.open("svZeroD_data", std::ios::out | std::ios::app); + out_file.open(svZeroD_output_dir + "svZeroD_data", std::ios::out); out_file<(variable_names[i])<<" "; @@ -156,7 +161,7 @@ void write_svZeroD_solution(const double* lpn_time, std::vector& lpn_sol out_file<<'\n'; } else { std::ofstream out_file; - out_file.open("svZeroD_data", std::ios::out | std::ios::app); + out_file.open(svZeroD_output_dir + "svZeroD_data", std::ios::out | std::ios::app); out_file<<*lpn_time<<" "; for (int i = 0; i < system_size; i++) { out_file<& surfID, double Q[], double P[]) { int nParam = 2; - const char* fileNames[2] = {"Q_svZeroD", "P_svZeroD"}; + const std::string fileNames[2] = {svZeroD_output_dir + "Q_svZeroD", + svZeroD_output_dir + "P_svZeroD"}; std::vector> R(nParam, std::vector(*nSrfs)); if (*nSrfs == 0) return; @@ -226,7 +232,7 @@ void print_svZeroD(int* nSrfs, const std::vector& surfID, double Q[], doubl // init_svZeroD //-------------- // -void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) +void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod, const std::string& appPath) { using namespace consts; @@ -241,6 +247,8 @@ void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) auto& cm = com_mod.cm; double dt = com_mod.dt; + svZeroD_output_dir = appPath; + build_svzero_coupled_bc_idxs(com_mod); if (cplBC.nSvZeroD_coupled_bc == 0) { throw std::runtime_error( diff --git a/Code/Source/solver/svZeroD_interface.h b/Code/Source/solver/svZeroD_interface.h index 1de624f85..7187b7542 100644 --- a/Code/Source/solver/svZeroD_interface.h +++ b/Code/Source/solver/svZeroD_interface.h @@ -17,7 +17,21 @@ void get_coupled_QP(ComMod& com_mod, double QCoupled[], double QnCoupled[], doub void print_svZeroD(int* nSrfs, const std::vector& surfID, double Q[], double P[]); -void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod); +/** + * @brief Set up the svZeroD model and its coupled boundary conditions. + * + * Builds the list of svZeroD-coupled boundaries, creates the svZeroD model + * from the interface data in \c com_mod.cplBC, applies the initial flows and + * pressures, and writes the header of the svZeroD state output file. + * + * @param[in,out] com_mod Simulation data; the coupled boundary conditions and + * \c cplBC bookkeeping are initialized here. + * @param[in] cm_mod MPI communicator data. + * @param[in] appPath Directory the simulation results are written to. The + * svZeroD output files (svZeroD_data, Q_svZeroD, P_svZeroD) are written + * there as well. + */ +void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod, const std::string& appPath); void calc_svZeroD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag); diff --git a/Code/Source/solver/txt.cpp b/Code/Source/solver/txt.cpp index 89f03c907..6598ef6b8 100644 --- a/Code/Source/solver/txt.cpp +++ b/Code/Source/solver/txt.cpp @@ -465,17 +465,22 @@ void write_boundary_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eq if (m == 1) { if (div) { tmp = fa.area; - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference) / tmp; + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, + std::nullopt, false) / + tmp; } else { if (pFlag && lTH) { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, std::nullopt, true, consts::MechanicalConfigurationType::reference); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, + std::nullopt, true); } else { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, + std::nullopt, false); } } } else if (m == nsd) { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, m - 1, + false); } else { throw std::runtime_error("WTXT only accepts 1 and nsd"); }