diff --git a/.gitignore b/.gitignore index 6c19192e5..2eb401315 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,27 @@ build +# CMake build artifacts +CMakeCache.txt +CMakeFiles/ +Makefile +cmake_install.cmake +*-build/ +*-prefix/ + +# Compiled binaries and libraries +Code/bin/ +Code/lib/ + +# Generated headers +Code/Source/Include/simvascular_options.h +Code/Source/Include/simvascular_version.h +Code/ThirdParty/*/simvascular_*.h + +# Vim swap files +*.swp +*.swo +*~ + # System files osx **/.DS_Store diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 94a97d6f7..88ace0580 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -181,6 +181,7 @@ set(CSRCS heatf.h heatf.cpp heats.h heats.cpp initialize.h initialize.cpp + Integrator.h Integrator.cpp l_elas.h l_elas.cpp lhsa.h lhsa.cpp ls.h ls.cpp @@ -191,7 +192,6 @@ set(CSRCS nn.h nn.cpp output.h output.cpp load_msh.h load_msh.cpp - pic.h pic.cpp post.h post.cpp read_files.h read_files.cpp read_msh.h read_msh.cpp diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 91ed18e48..95b2cdd6b 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -40,6 +40,40 @@ class LinearAlgebra; +/** + * @brief Represents solution variables at a single time level + * + * Contains the three primary solution arrays used in time integration: + * A (time derivative), D (integrated variable), and Y (variable) + */ +struct Solution { + Array A; ///< Time derivative (acceleration in structural mechanics) + Array D; ///< Integrated variable (displacement in structural mechanics) + Array Y; ///< Variable (velocity in structural mechanics) + + // Semantic getters for improved readability + Array& get_acceleration() { return A; } + const Array& get_acceleration() const { return A; } + + Array& get_velocity() { return Y; } + const Array& get_velocity() const { return Y; } + + Array& get_displacement() { return D; } + const Array& get_displacement() const { return D; } +}; + +/** + * @brief Holds solution state at old and current time levels + * + * Contains solution arrays at two time levels for time integration: + * - old: Previous converged solution at time n + * - current: Current solution being computed at time n+1 + */ +struct SolutionStates { + Solution old; ///< Previous converged solution at time n (Ao, Do, Yo) + Solution current; ///< Current solution being computed at time n+1 (An, Dn, Yn) +}; + /// @brief Fourier coefficients that are used to specify unsteady BCs // class fcType @@ -1749,18 +1783,6 @@ class ComMod { /// @brief RIS mapping array, with global (total) enumeration std::vector grisMapList; - /// @brief Old time derivative of variables (acceleration); known result at current time step - Array Ao; - - /// @brief New time derivative of variables (acceleration); unknown result at next time step - Array An; - - /// @brief Old integrated variables (displacement) - Array Do; - - /// @brief New integrated variables (displacement) - Array Dn; - /// @brief Residual vector Array R; @@ -1770,12 +1792,6 @@ class ComMod { /// @brief Position vector of mesh nodes (in ref config) Array x; - /// @brief Old variables (velocity); known result at current time step - Array Yo; - - /// @brief New variables (velocity); unknown result at next time step - Array Yn; - /// @brief Body force Array Bf; diff --git a/Code/Source/solver/pic.cpp b/Code/Source/solver/Integrator.cpp similarity index 50% rename from Code/Source/solver/pic.cpp rename to Code/Source/solver/Integrator.cpp index 34a274def..a041f5a1e 100644 --- a/Code/Source/solver/pic.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -1,37 +1,641 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. // SPDX-License-Identifier: BSD-3-Clause -// The code here replicates the Fortran code in PIC.f. -// -// See the publications below, section 4.4 for theory and derivation: -// 1. Bazilevs, et al. "Isogeometric fluid-structure interaction: -// theory, algorithms, and computations.", Computational Mechanics, -// 43 (2008): 3-37. doi: 10.1007/s00466-008-0315-x -// 2. Bazilevs, et al. "Variational multiscale residual-based -// turbulence modeling for large eddy simulation of incompressible -// flows.", CMAME (2007) - -#include "pic.h" - -#include "Simulation.h" +#include "Integrator.h" #include "all_fun.h" +#include "bf.h" #include "cep_ion.h" +#include "contact.h" +#include "eq_assem.h" +#include "fs.h" +#include "ls.h" #include "nn.h" +#include "output.h" +#include "ris.h" +#include "set_bc.h" +#include "ustruct.h" #include "utils.h" -#include "mpi.h" +#include +#include +#include + +using namespace consts; + +//------------------------ +// Integrator Constructor +//------------------------ +Integrator::Integrator(Simulation* simulation, SolutionStates&& solutions) + : simulation_(simulation), solutions_(std::move(solutions)), newton_count_(0) +{ + initialize_arrays(); +} + +//------------------------ +// Integrator Destructor +//------------------------ +Integrator::~Integrator() { + // Arrays will be automatically cleaned up +} + +//------------------------ +// initialize_arrays +//------------------------ +void Integrator::initialize_arrays() { + auto& com_mod = simulation_->com_mod; + int tDof = com_mod.tDof; + int tnNo = com_mod.tnNo; + int nFacesLS = com_mod.nFacesLS; + + Ag_.resize(tDof, tnNo); + Yg_.resize(tDof, tnNo); + Dg_.resize(tDof, tnNo); + solutions_.current.get_acceleration().resize(tDof, tnNo); + solutions_.current.get_displacement().resize(tDof, tnNo); + solutions_.current.get_velocity().resize(tDof, tnNo); + res_.resize(nFacesLS); + incL_.resize(nFacesLS); + + // old solution already initialized via move in constructor + // Initialize current solution from old solution + solutions_.current.get_acceleration() = solutions_.old.get_acceleration(); + solutions_.current.get_displacement() = solutions_.old.get_displacement(); + solutions_.current.get_velocity() = solutions_.old.get_velocity(); +} + +//------------------------ +// step +//------------------------ +/// @brief Execute one Newton iteration loop for the current time step +bool Integrator::step() { + using namespace consts; + + auto& com_mod = simulation_->com_mod; + auto& cm_mod = simulation_->cm_mod; + auto& cep_mod = simulation_->get_cep_mod(); + + int& cTS = com_mod.cTS; + int& cEq = com_mod.cEq; + + #define n_debug_integrator_step + #ifdef debug_integrator_step + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + + // Newton iteration loop + newton_count_ = 1; + int reply; + int iEqOld; + + // Looping over Newton iterations + while (true) { + #ifdef debug_integrator_step + dmsg << "---------- Newton Iteration " + std::to_string(newton_count_) << " -----------" << std::endl; + dmsg << "cEq: " << cEq; + dmsg << "com_mod.eq[cEq].sym: " << com_mod.eq[cEq].sym; + #endif + + istr_ = "_" + std::to_string(cTS) + "_" + std::to_string(newton_count_); + iEqOld = cEq; + auto& eq = com_mod.eq[cEq]; + + if (com_mod.cplBC.coupled && cEq == 0) { + #ifdef debug_integrator_step + dmsg << "Set coupled BCs " << std::endl; + #endif + set_bc::set_bc_cpl(com_mod, cm_mod, solutions_); + set_bc::set_bc_dir(com_mod, solutions_); + } + + // Initiator step for Generalized α-Method (quantities at n+am, n+af). + initiator_step(); + + if (com_mod.Rd.size() != 0) { + com_mod.Rd = 0.0; + com_mod.Kd = 0.0; + } + + // Allocate com_mod.R and com_mod.Val arrays + allocate_linear_system(eq); + + // Compute body forces + set_body_forces(); + + // Assemble equations + assemble_equations(); + + // Treatment of boundary conditions on faces + apply_boundary_conditions(); + + // Synchronize R across processes + if (!eq.assmTLS) { + #ifdef debug_integrator_step + dmsg << "Synchronize R across processes ..." << std::endl; + #endif + all_fun::commu(com_mod, com_mod.R); + } + + // Update residual in displacement equation for USTRUCT phys + #ifdef debug_integrator_step + dmsg << "com_mod.sstEq: " << com_mod.sstEq; + #endif + if (com_mod.sstEq) { + ustruct::ustruct_r(com_mod, Yg_); + } + + // Set the residual of the continuity equation to 0 on edge nodes + if (std::set{Equation_stokes, Equation_fluid, Equation_ustruct, Equation_FSI}.count(eq.phys) != 0) { + #ifdef debug_integrator_step + dmsg << "thood_val_rc ..." << std::endl; + #endif + fs::thood_val_rc(com_mod); + } + + // Treat Neumann boundaries that are not deforming + #ifdef debug_integrator_step + dmsg << "set_bc_undef_neu ..." << std::endl; + #endif + set_bc::set_bc_undef_neu(com_mod); + + // Update residual and increment arrays + update_residual_arrays(eq); + + // Solve equation + solve_linear_system(); + + // Solution is obtained, now updating (Corrector) and check for convergence + bool all_converged = corrector_and_check_convergence(); + + // Check if all equations converged + if (all_converged) { + #ifdef debug_integrator_step + dmsg << ">>> All OK" << std::endl; + dmsg << "iEqOld: " << iEqOld + 1; + #endif + return true; + } + + output::output_result(simulation_, com_mod.timeP, 2, iEqOld); + newton_count_ += 1; + } // End of Newton iteration loop + + return false; +} + +//------------------------ +// initiator_step +//------------------------ +void Integrator::initiator_step() { + #define n_debug_integrator_step + #ifdef debug_integrator_step + DebugMsg dmsg(__func__, simulation_->com_mod.cm.idcm()); + dmsg << "Initiator step ..." << std::endl; + #endif + + initiator(Ag_, Yg_, Dg_); + + // Debug output + Ag_.write("Ag_pic" + istr_); + Yg_.write("Yg_pic" + istr_); + Dg_.write("Dg_pic" + istr_); + solutions_.current.get_velocity().write("solutions_.current.Ypic" + istr_); +} + +//------------------------ +// allocate_linear_system +//------------------------ +void Integrator::allocate_linear_system(eqType& eq) { + #define n_debug_integrator_step + #ifdef debug_integrator_step + DebugMsg dmsg(__func__, simulation_->com_mod.cm.idcm()); + dmsg << "Allocating the RHS and LHS" << std::endl; + #endif + + ls_ns::ls_alloc(simulation_->com_mod, eq); + + // Debug output + simulation_->com_mod.Val.write("Val_alloc" + istr_); +} + +//------------------------ +// set_body_forces +//------------------------ +void Integrator::set_body_forces() { + #define n_debug_integrator_step + #ifdef debug_integrator_step + DebugMsg dmsg(__func__, simulation_->com_mod.cm.idcm()); + dmsg << "Set body forces ..." << std::endl; + #endif + + bf::set_bf(simulation_->com_mod, Dg_); + + // Debug output + simulation_->com_mod.Val.write("Val_bf" + istr_); +} + +//------------------------ +// assemble_equations +//------------------------ +void Integrator::assemble_equations() { + auto& com_mod = simulation_->com_mod; + auto& cep_mod = simulation_->get_cep_mod(); + + #define n_debug_integrator_step + #ifdef debug_integrator_step + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg << "Assembling equation: " << com_mod.eq[com_mod.cEq].sym; + #endif + + for (int iM = 0; iM < com_mod.nMsh; iM++) { + eq_assem::global_eq_assem(com_mod, cep_mod, com_mod.msh[iM], Ag_, Yg_, Dg_, solutions_); + } + + // Debug output + com_mod.R.write("R_as" + istr_); + com_mod.Val.write("Val_as" + istr_); +} + +//------------------------ +// apply_boundary_conditions +//------------------------ +void Integrator::apply_boundary_conditions() { + auto& com_mod = simulation_->com_mod; + auto& cm_mod = simulation_->cm_mod; + + #define n_debug_integrator_step + #ifdef debug_integrator_step + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg << "Apply boundary conditions ..." << std::endl; + #endif + + Yg_.write("Yg_vor_neu" + istr_); + Dg_.write("Dg_vor_neu" + istr_); + + // Apply Neumman or Traction boundary conditions + set_bc::set_bc_neu(com_mod, cm_mod, Yg_, Dg_, solutions_); + + // Apply CMM BC conditions + if (!com_mod.cmmInit) { + set_bc::set_bc_cmm(com_mod, cm_mod, Ag_, Dg_, solutions_); + } + + // Apply weakly applied Dirichlet BCs + set_bc::set_bc_dir_w(com_mod, Yg_, Dg_, solutions_); + + if (com_mod.risFlag) { + ris::ris_resbc(com_mod, Yg_, Dg_, solutions_); + } + + if (com_mod.ris0DFlag) { + ris::ris0d_bc(com_mod, cm_mod, Yg_, Dg_, solutions_); + } + + // Apply contact model and add its contribution to residual + if (com_mod.iCntct) { + contact::construct_contact_pnlty(com_mod, cm_mod, Dg_); + } + + // Debug output + com_mod.Val.write("Val_neu" + istr_); + com_mod.R.write("R_neu" + istr_); + Yg_.write("Yg_neu" + istr_); + Dg_.write("Dg_neu" + istr_); +} + +//------------------------ +// solve_linear_system +//------------------------ +void Integrator::solve_linear_system() { + auto& com_mod = simulation_->com_mod; + auto& eq = com_mod.eq[com_mod.cEq]; + + #define n_debug_integrator_step + #ifdef debug_integrator_step + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg << "Solving equation: " << eq.sym; + #endif + + ls_ns::ls_solve(com_mod, eq, incL_, res_); + + // Debug output + com_mod.Val.write("Val_solve" + istr_); + com_mod.R.write("R_solve" + istr_); +} + +//------------------------ +// corrector_and_check_convergence +//------------------------ +bool Integrator::corrector_and_check_convergence() { + auto& com_mod = simulation_->com_mod; + + #define n_debug_integrator_step + #ifdef debug_integrator_step + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg << "Update corrector ..." << std::endl; + #endif + + corrector(); + + // Debug output + solutions_.current.get_velocity().write("solutions_.current.Ycorrector" + istr_); + + // Check if all equations converged + return std::count_if(com_mod.eq.begin(), com_mod.eq.end(), + [](eqType& eq) { return eq.ok; }) == com_mod.eq.size(); +} + +//------------------------ +// update_residual_arrays +//------------------------ +void Integrator::update_residual_arrays(eqType& eq) { + auto& com_mod = simulation_->com_mod; + int nFacesLS = com_mod.nFacesLS; + double dt = com_mod.dt; + + #define n_debug_integrator_step + #ifdef debug_integrator_step + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg << "Update res() and incL ..." << std::endl; + #endif + + incL_ = 0; + if (eq.phys == Equation_mesh) { + incL_(nFacesLS - 1) = 1; + } + + if (com_mod.cmmInit) { + incL_(nFacesLS - 1) = 1; + } + + for (int iBc = 0; iBc < eq.nBc; iBc++) { + int i = eq.bc[iBc].lsPtr; + if (i != -1) { + // Resistance term for coupled Neumann BC tangent contribution + res_(i) = eq.gam * dt * eq.bc[iBc].r; + incL_(i) = 1; + } + } +} + +//------------------------ +// predictor (picp) +//------------------------ +/// @brief Predictor step for next time step +/// +/// Modifies: +/// pS0 +/// Ad +/// Ao +/// Yo +/// Do +/// An +/// Yn +/// Dn +/// +void Integrator::predictor() +{ + using namespace consts; + + auto& com_mod = simulation_->com_mod; + + #define n_debug_picp + #ifdef debug_picp + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "pstEq: " << com_mod.pstEq; + #endif + + // Variables for prestress calculations + auto& pS0 = com_mod.pS0; + auto& pSn = com_mod.pSn; + + // time derivative of displacement + auto& Ad = com_mod.Ad; + + auto& Ao = solutions_.old.get_acceleration(); + auto& An = solutions_.current.get_acceleration(); + auto& Yo = solutions_.old.get_velocity(); + auto& Yn = solutions_.current.get_velocity(); + auto& Do = solutions_.old.get_displacement(); + auto& Dn = solutions_.current.get_displacement(); + + // Prestress initialization + if (com_mod.pstEq) { + pS0 = pS0 + pSn; + Ao = 0.0; + Yo = 0.0; + Do = 0.0; + } + + // IB treatment: Set dirichlet BC and update traces. For explicit + // coupling, compute FSI forcing and freeze it for the time step. + // For implicit coupling, project IB displacement on background + // mesh and predict quantities at next time step + // + // [NOTE] not implemented. + /* + if (ibFlag) { + // Set IB Dirichlet BCs + CALL IB_SETBCDIR(ib.Yb, ib.Ubo) + + // Update IB location and tracers + CALL IB_UPDATE(Do) + + if (ib.cpld == ibCpld_E) { + // FSI forcing for immersed bodies (explicit coupling) + CALL IB_CALCFFSI(Ao, Yo, Do, ib.Auo, ib.Ubo) + + } else { if (ib.cpld == ibCpld_I) { + // Project IB displacement (Ubn) to background mesh + CALL IB_PRJCTU(Do) + + // Predictor step for implicit coupling + CALL IB_PICP() + } + } + */ + + const auto& dt = com_mod.dt; + #ifdef debug_picp + dmsg << "dt: " << dt; + dmsg << "dFlag: " << com_mod.dFlag; + #endif + + for (int iEq = 0; iEq < com_mod.nEq; iEq++) { + auto& eq = com_mod.eq[iEq]; + int s = eq.s; // start row + int e = eq.e; // end row + + #ifdef debug_picp + dmsg << "----- iEq " << iEq << " -----"; + dmsg << "s: " << s; + dmsg << "e: " << e; + dmsg << "eq.gam: " << eq.gam; + dmsg << "coef: " << coef; + #endif + + // [TODO:DaveP] careful here with s amd e. + double coef = (eq.gam - 1.0) / eq.gam; + for (int i = s; i <= e; i++) { + for (int j = 0; j < Ao.ncols(); j++) { + // eqn 87 of Bazilevs 2007 + An(i,j) = Ao(i,j) * coef; + } + } + + // electrophysiology + if (eq.phys == Equation_CEP) { + cep_ion::cep_integ(simulation_, iEq, e, solutions_); + } + + // eqn 86 of Bazilevs 2007 + Yn.set_rows(s,e, Yo.rows(s,e)); + + if (com_mod.dFlag) { + + // struct, lElas, FSI (struct, mesh) + if (!com_mod.sstEq) { + double coef = dt*dt*(0.5*eq.gam - eq.beta) / (eq.gam - 1.0); + Dn.set_rows(s,e, Do.rows(s,e) + Yn.rows(s,e)*dt + An.rows(s,e)*coef); + + // ustruct, FSI + // + } else { + + if (eq.phys == Equation_ustruct || eq.phys == Equation_FSI) { + double coef = (eq.gam - 1.0) / eq.gam; + Ad = Ad*coef; + Dn.set_rows(s,e, Do.rows(s,e)); + + } else if (eq.phys == Equation_mesh) { + double coef = dt*dt*(0.5*eq.gam - eq.beta) / (eq.gam - 1.0); + Dn.set_rows(s,e, Do.rows(s,e) + Yn.rows(s,e)*dt + An.rows(s,e)*coef); + } + } + } else { + Dn.set_rows(s,e, Do.rows(s,e)); + } + } +} + +//------------------------ +// initiator +//------------------------ +/// @brief Initiator for Generalized α-Method +/// +/// Uses Generalized α− Method for time stepping. +/// +/// Modifes Ag from combination of An and Ao defined by coefs from eq.am, eq.af, +/// Ag = (1 - eq.am) * Ao + eq.am * An +/// Yg = (1 - eq.af) * Yo + eq.af * Yn +/// Dg = (1 - eq.af) * Do + eq.af * Dn +/// +/// Modifies: +/// Ag - acceleration +/// Yg - velocity +/// Dg - displacement +/// +void Integrator::initiator(Array& Ag, Array& Yg, Array& Dg) +{ + using namespace consts; + + auto& com_mod = simulation_->com_mod; + + #define n_debug_initiator + #ifdef debug_initiator + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + + const int cEq = com_mod.cEq; + const int tnNo = com_mod.tnNo; + auto& eq = com_mod.eq[cEq]; + auto& dof = com_mod.dof; + eq.itr = eq.itr + 1; + + // [NOTE] Setting gobal variable 'dof'. + dof = eq.dof; + #ifdef debug_initiator + dmsg << "cEq: " << cEq; + dmsg << "eq.itr: " << eq.itr; + dmsg << "dof: " << dof; + dmsg << "tnNo: " << tnNo; + dmsg << "com_mod.pstEq: " << com_mod.pstEq; + #endif + + const auto& Ao = solutions_.old.get_acceleration(); + const auto& An = solutions_.current.get_acceleration(); + const auto& Do = solutions_.old.get_displacement(); + const auto& Dn = solutions_.current.get_displacement(); + const auto& Yo = solutions_.old.get_velocity(); + const auto& Yn = solutions_.current.get_velocity(); + + for (int i = 0; i < com_mod.nEq; i++) { + auto& eq = com_mod.eq[i]; + int s = eq.s; + int e = eq.e; + Vector coef(4); + coef(0) = 1.0 - eq.am; + coef(1) = eq.am; + coef(2) = 1.0 - eq.af; + coef(3) = eq.af; + #ifdef debug_initiator + dmsg << "s: " << s; + dmsg << "e: " << e; + dmsg << "coef: " << coef[0] << " " << coef[1] << " " << coef[2] << " " << coef[3]; + #endif + + if ((eq.phys == Equation_heatF) && (com_mod.usePrecomp)){ + for (int a = 0; a < tnNo; a++) { + for (int j = 0; j < com_mod.nsd; j++) { + //Ag(j, a) = An(j, a); + //Yg(j, a) = Yn(j, a); + //Dg(j, a) = Dn(j, a); + Ag(j, a) = Ao(j, a) * coef(0) + An(j, a) * coef(1); + Yg(j, a) = Yo(j, a) * coef(2) + Yn(j, a) * coef(3); + Dg(j, a) = Do(j, a) * coef(2) + Dn(j, a) * coef(3); + } + } + for (int a = 0; a < tnNo; a++) { + for (int j = s; j <= e; j++) { + Ag(j, a) = Ao(j, a) * coef(0) + An(j, a) * coef(1); + Yg(j, a) = Yo(j, a) * coef(2) + Yn(j, a) * coef(3); + Dg(j, a) = Do(j, a) * coef(2) + Dn(j, a) * coef(3); + } + } + } else { + for (int a = 0; a < tnNo; a++) { + for (int j = s; j <= e; j++) { + // eqn 89 of Bazilevs 2007 + Ag(j, a) = Ao(j, a) * coef(0) + An(j, a) * coef(1); -#include -#include + // eqn 90 of Bazilevs 2007 + Yg(j, a) = Yo(j, a) * coef(2) + Yn(j, a) * coef(3); -namespace pic { + Dg(j, a) = Do(j, a) * coef(2) + Dn(j, a) * coef(3); + } + } + } + } -/// @brief This is the corrector. Decision for next eqn is also made here (modifies cEq global). + // prestress + if (com_mod.pstEq) { + com_mod.pSn = 0.0; + com_mod.pSa = 0.0; + } +} +//------------------------ +// corrector +//------------------------ +/// @brief Corrector with convergence check +/// +/// Decision for next eqn is also made here (modifies cEq global). /// /// Modifies: /// \code {.cpp} /// com_mod.Ad -/// com_mod.An +/// solutions_.current.A (member variable) /// com_mod.Dn /// com_mod.Yn /// cep_mod.Xion @@ -40,19 +644,19 @@ namespace pic { /// com_mod.pSn /// /// com_mod.cEq -/// eq.FSILS.RI.iNorm +/// eq.FSILS.RI.iNorm /// eq.pNorm /// \endcode // -void picc(Simulation* simulation) +void Integrator::corrector() { using namespace consts; - auto& com_mod = simulation->com_mod; - auto& cep_mod = simulation->get_cep_mod(); + auto& com_mod = simulation_->com_mod; + auto& cep_mod = simulation_->get_cep_mod(); - #define n_debug_picc - #ifdef debug_picc + #define n_debug_corrector + #ifdef debug_corrector DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); #endif @@ -67,10 +671,10 @@ void picc(Simulation* simulation) auto& cEq = com_mod.cEq; auto& eq = com_mod.eq[cEq]; - auto& An = com_mod.An; + auto& An = solutions_.current.get_acceleration(); auto& Ad = com_mod.Ad; - auto& Dn = com_mod.Dn; - auto& Yn = com_mod.Yn; + auto& Dn = solutions_.current.get_displacement(); + auto& Yn = solutions_.current.get_velocity(); auto& pS0 = com_mod.pS0; auto& pSa = com_mod.pSa; @@ -86,7 +690,7 @@ void picc(Simulation* simulation) coef[2] = 1.0 / eq.am; coef[3] = eq.af*coef[0]*coef[2]; - #ifdef debug_picc + #ifdef debug_corrector dmsg << "cEq: " << cEq; dmsg << "s: " << s; dmsg << "e: " << e; @@ -130,23 +734,23 @@ void picc(Simulation* simulation) for (int i = 0; i < e-s+1; i++) { // eqn 94 of Bazilevs 2007 // here, -R contains the acceleration update (obtained from Newton solve))? An(i+s,a) = An(i+s,a) - R(i,a); - + // eqn 95 of Bazilevs 2007 Yn(i+s,a) = Yn(i+s,a) - R(i,a)*coef[0]; - + Dn(i+s,a) = Dn(i+s,a) - R(i,a)*coef[1]; } } } if (std::set{Equation_stokes, Equation_fluid, Equation_ustruct, Equation_FSI}.count(eq.phys) != 0) { - pic_eth(simulation); + corrector_taylor_hood(); } if (eq.phys == Equation_FSI) { int s = com_mod.eq[1].s; int e = com_mod.eq[1].e; - #ifdef debug_picc + #ifdef debug_corrector dmsg << "eq.phys == Equation_FSI "; dmsg << "com_mod.eq[1].sym: " << com_mod.eq[1].sym; dmsg << "s: " << s; @@ -154,8 +758,8 @@ void picc(Simulation* simulation) #endif for (int Ac = 0; Ac < tnNo; Ac++) { - if (all_fun::is_domain(com_mod, eq, Ac, Equation_struct) || - all_fun::is_domain(com_mod, eq, Ac, Equation_ustruct) || + if (all_fun::is_domain(com_mod, eq, Ac, Equation_struct) || + all_fun::is_domain(com_mod, eq, Ac, Equation_ustruct) || all_fun::is_domain(com_mod, eq, Ac, Equation_lElas)) { for (int i = 0; i < e-s+1; i++) { An(i+s,Ac) = An(i,Ac); @@ -167,7 +771,7 @@ void picc(Simulation* simulation) } // Update Xion for cardiac electrophysiology - // + // if (eq.phys == Equation_CEP) { int s = eq.s; for (int a = 0; a < tnNo; a++) { @@ -203,9 +807,6 @@ void picc(Simulation* simulation) } } - // IB treatment - //if (ibFlag) CALL IB_PICC() - // Computes norms and check for convergence of Newton iterations double eps = std::numeric_limits::epsilon(); @@ -215,14 +816,14 @@ void picc(Simulation* simulation) if (utils::is_zero(eq.iNorm)) { eq.iNorm = eq.FSILS.RI.iNorm; - #ifdef debug_picc + #ifdef debug_corrector dmsg << "eq.iNorm: " << eq.iNorm; #endif } if (eq.itr == 1) { eq.pNorm = eq.FSILS.RI.iNorm / eq.iNorm; - #ifdef debug_picc + #ifdef debug_corrector dmsg << "eq.itr: " << eq.itr; dmsg << "eq.pNorm: " << eq.pNorm; #endif @@ -234,7 +835,7 @@ void picc(Simulation* simulation) bool l3 = (r1 <= eq.tol*eq.pNorm); bool l4 = (eq.itr >= eq.minItr); - #ifdef debug_picc + #ifdef debug_corrector dmsg << "eq.itr: " << eq.itr; dmsg << "eq.minItr: " << eq.minItr; dmsg << "r1: " << r1; @@ -246,7 +847,7 @@ void picc(Simulation* simulation) if (l1 || ((l2 || l3) && l4)) { eq.ok = true; - #ifdef debug_picc + #ifdef debug_corrector dmsg << "eq.ok: " << eq.ok; dmsg << "com_mod.eq[0].ok: " << com_mod.eq[0].ok; dmsg << "com_mod.eq[1].ok: " << com_mod.eq[1].ok; @@ -255,7 +856,7 @@ void picc(Simulation* simulation) auto& eqs = com_mod.eq; if (std::count_if(eqs.begin(),eqs.end(),[](eqType& eq){return eq.ok;}) == eqs.size()) { - #ifdef debug_picc + #ifdef debug_corrector dmsg << "all ok"; #endif return; @@ -264,7 +865,7 @@ void picc(Simulation* simulation) if (eq.coupled) { cEq = cEq + 1; - #ifdef debug_picc + #ifdef debug_corrector dmsg << "eq " << " coupled "; dmsg << "1st update cEq: " << cEq; #endif @@ -296,29 +897,30 @@ void picc(Simulation* simulation) cEq = cEq + 1; } } - #ifdef debug_picc + #ifdef debug_corrector dmsg << "eq " << " coupled "; dmsg << "2nd update cEq: " << cEq; #endif } -//--------- -// pic_eth -//--------- -// Pressure correction at edge nodes for Taylor-Hood type element. -// Here, we interpolate pressure at the edge nodes by interpolating -// using a reduced basis (such as P1) applied on element vertices -// (i.e., corner nodes). For e.g., for a P2 element, pressure is -// interpolated at the edge nodes using P1 vertices. -// -// Modifies: com_mod.Yn -// -void pic_eth(Simulation* simulation) +//------------------------ +// corrector_taylor_hood +//------------------------ +/// @brief Pressure correction at edge nodes for Taylor-Hood type element +/// +/// Here, we interpolate pressure at the edge nodes by interpolating +/// using a reduced basis (such as P1) applied on element vertices +/// (i.e., corner nodes). For e.g., for a P2 element, pressure is +/// interpolated at the edge nodes using P1 vertices. +/// +/// Modifies: solutions_.current.Y +/// +void Integrator::corrector_taylor_hood() { using namespace consts; - auto& com_mod = simulation->com_mod; - auto& cep_mod = simulation->get_cep_mod(); + auto& com_mod = simulation_->com_mod; + auto& cep_mod = simulation_->get_cep_mod(); const int nsd = com_mod.nsd; const int tnNo = com_mod.tnNo; @@ -327,7 +929,7 @@ void pic_eth(Simulation* simulation) const auto& eq = com_mod.eq[cEq]; auto& cDmn = com_mod.cDmn; - auto& Yn = com_mod.Yn; + auto& Yn = solutions_.current.get_velocity(); // Check for something ... // @@ -358,9 +960,9 @@ void pic_eth(Simulation* simulation) int eNoNq = msh.fs[1].eNoN; Array xl(nsd,eNoN), xql(nsd,eNoNq), Nqx(nsd,eNoNq); - Vector pl(eNoNq), Nq(eNoNq); + Vector pl(eNoNq), Nq(eNoNq); - Vector xp(nsd), xi0(nsd), xi(nsd); + Vector xp(nsd), xi0(nsd), xi(nsd); Array ksix(nsd,nsd); for (int g = 0; g < msh.fs[1].nG; g++) { @@ -373,8 +975,8 @@ void pic_eth(Simulation* simulation) for (int e = 0; e < msh.nEl; e++) { cDmn = all_fun::domain(com_mod, msh, cEq, e); // setting global cDmn - if ((eq.dmn[cDmn].phys != Equation_stokes) && - (eq.dmn[cDmn].phys != Equation_fluid) && + if ((eq.dmn[cDmn].phys != Equation_stokes) && + (eq.dmn[cDmn].phys != Equation_fluid) && (eq.dmn[cDmn].phys != Equation_ustruct)) { continue; } @@ -403,7 +1005,7 @@ void pic_eth(Simulation* simulation) nn::gnn(eNoNq, nsd, nsd, Nx, xql, Nqx, Jac, ksix); if (utils::is_zero(Jac)) { - throw std::runtime_error("[pic_eth] Jacobian for element " + std::to_string(e) + " is < 0."); + throw std::runtime_error("[corrector_taylor_hood] Jacobian for element " + std::to_string(e) + " is < 0."); } } @@ -438,254 +1040,3 @@ void pic_eth(Simulation* simulation) } } } - -//------ -// pici -//------ -// This is the initiator. -// -// Uses Generalized α− Method for time stepping. -// -// Modifes Ag from combination of An and Ao defined by coefs from eq.am, eq.af, -// Ag = (1 - eq.am) * Ao + eq.am * An -// Yg = (1 - eq.af) * Yo + eq.af * Yn -// Dg = (1 - eq.af) * Do + eq.af * Dn -// -// Modifies: -// Ag - acceleration -// Yg - velocity -// Dg - displacement -// -void pici(Simulation* simulation, Array& Ag, Array& Yg, Array& Dg) -{ - using namespace consts; - - auto& com_mod = simulation->com_mod; - - #define n_debug_pici - #ifdef debug_pici - DebugMsg dmsg(__func__, com_mod.cm.idcm()); - dmsg.banner(); - #endif - - const int cEq = com_mod.cEq; - const int tnNo = com_mod.tnNo; - auto& eq = com_mod.eq[cEq]; - auto& dof = com_mod.dof; - eq.itr = eq.itr + 1; - - // [NOTE] Setting gobal variable 'dof'. - dof = eq.dof; - #ifdef debug_pici - dmsg << "cEq: " << cEq; - dmsg << "eq.itr: " << eq.itr; - dmsg << "dof: " << dof; - dmsg << "tnNo: " << tnNo; - dmsg << "com_mod.pstEq: " << com_mod.pstEq; - #endif - - const auto& Ao = com_mod.Ao; - const auto& An = com_mod.An; - const auto& Do = com_mod.Do; - const auto& Dn = com_mod.Dn; - const auto& Yo = com_mod.Yo; - const auto& Yn = com_mod.Yn; - - for (int i = 0; i < com_mod.nEq; i++) { - auto& eq = com_mod.eq[i]; - int s = eq.s; - int e = eq.e; - Vector coef(4); - coef(0) = 1.0 - eq.am; - coef(1) = eq.am; - coef(2) = 1.0 - eq.af; - coef(3) = eq.af; - #ifdef debug_pici - dmsg << "s: " << s; - dmsg << "e: " << e; - dmsg << "coef: " << coef[0] << " " << coef[1] << " " << coef[2] << " " << coef[3]; - #endif - - if ((eq.phys == Equation_heatF) && (com_mod.usePrecomp)){ - for (int a = 0; a < tnNo; a++) { - for (int j = 0; j < com_mod.nsd; j++) { - //Ag(j, a) = An(j, a); - //Yg(j, a) = Yn(j, a); - //Dg(j, a) = Dn(j, a); - Ag(j, a) = Ao(j, a) * coef(0) + An(j, a) * coef(1); - Yg(j, a) = Yo(j, a) * coef(2) + Yn(j, a) * coef(3); - Dg(j, a) = Do(j, a) * coef(2) + Dn(j, a) * coef(3); - } - } - for (int a = 0; a < tnNo; a++) { - for (int j = s; j <= e; j++) { - Ag(j, a) = Ao(j, a) * coef(0) + An(j, a) * coef(1); - Yg(j, a) = Yo(j, a) * coef(2) + Yn(j, a) * coef(3); - Dg(j, a) = Do(j, a) * coef(2) + Dn(j, a) * coef(3); - } - } - } else { - for (int a = 0; a < tnNo; a++) { - for (int j = s; j <= e; j++) { - // eqn 89 of Bazilevs 2007 - Ag(j, a) = Ao(j, a) * coef(0) + An(j, a) * coef(1); - - // eqn 90 of Bazilevs 2007 - Yg(j, a) = Yo(j, a) * coef(2) + Yn(j, a) * coef(3); - - Dg(j, a) = Do(j, a) * coef(2) + Dn(j, a) * coef(3); - } - } - } - } - - // prestress - if (com_mod.pstEq) { - com_mod.pSn = 0.0; - com_mod.pSa = 0.0; - } -} - -//------ -// picp -//------ -// This is the predictor. -// -// Modifies: -// pS0 -// Ad -// Ao -// Yo -// Do -// An -// Yn -// Dn -// -void picp(Simulation* simulation) -{ - using namespace consts; - - auto& com_mod = simulation->com_mod; - - #define n_debug_picp - #ifdef debug_picp - DebugMsg dmsg(__func__, com_mod.cm.idcm()); - dmsg.banner(); - dmsg << "pstEq: " << com_mod.pstEq; - #endif - - // Variables for prestress calculations - auto& pS0 = com_mod.pS0; - auto& pSn = com_mod.pSn; - - // time derivative of displacement - auto& Ad = com_mod.Ad; - - auto& Ao = com_mod.Ao; - auto& An = com_mod.An; - auto& Yo = com_mod.Yo; - auto& Yn = com_mod.Yn; - auto& Do = com_mod.Do; - auto& Dn = com_mod.Dn; - - // Prestress initialization - if (com_mod.pstEq) { - pS0 = pS0 + pSn; - Ao = 0.0; - Yo = 0.0; - Do = 0.0; - } - - // IB treatment: Set dirichlet BC and update traces. For explicit - // coupling, compute FSI forcing and freeze it for the time step. - // For implicit coupling, project IB displacement on background - // mesh and predict quantities at next time step - // - // [NOTE] not implemented. - /* - if (ibFlag) { - // Set IB Dirichlet BCs - CALL IB_SETBCDIR(ib.Yb, ib.Ubo) - - // Update IB location and tracers - CALL IB_UPDATE(Do) - - if (ib.cpld == ibCpld_E) { - // FSI forcing for immersed bodies (explicit coupling) - CALL IB_CALCFFSI(Ao, Yo, Do, ib.Auo, ib.Ubo) - - } else { if (ib.cpld == ibCpld_I) { - // Project IB displacement (Ubn) to background mesh - CALL IB_PRJCTU(Do) - - // Predictor step for implicit coupling - CALL IB_PICP() - } - } - */ - - const auto& dt = com_mod.dt; - #ifdef debug_picp - dmsg << "dt: " << dt; - dmsg << "dFlag: " << com_mod.dFlag; - #endif - - for (int iEq = 0; iEq < com_mod.nEq; iEq++) { - auto& eq = com_mod.eq[iEq]; - int s = eq.s; // start row - int e = eq.e; // end row - - #ifdef debug_picp - dmsg << "----- iEq " << iEq << " -----"; - dmsg << "s: " << s; - dmsg << "e: " << e; - dmsg << "eq.gam: " << eq.gam; - dmsg << "coef: " << coef; - #endif - - // [TODO:DaveP] careful here with s amd e. - double coef = (eq.gam - 1.0) / eq.gam; - for (int i = s; i <= e; i++) { - for (int j = 0; j < Ao.ncols(); j++) { - // eqn 87 of Bazilevs 2007 - An(i,j) = Ao(i,j) * coef; - } - } - - // electrophysiology - if (eq.phys == Equation_CEP) { - cep_ion::cep_integ(simulation, iEq, e, Do); - } - - // eqn 86 of Bazilevs 2007 - Yn.set_rows(s,e, Yo.rows(s,e)); - - if (com_mod.dFlag) { - - // struct, lElas, FSI (struct, mesh) - if (!com_mod.sstEq) { - double coef = dt*dt*(0.5*eq.gam - eq.beta) / (eq.gam - 1.0); - Dn.set_rows(s,e, Do.rows(s,e) + Yn.rows(s,e)*dt + An.rows(s,e)*coef); - - // ustruct, FSI - // - } else { - - if (eq.phys == Equation_ustruct || eq.phys == Equation_FSI) { - double coef = (eq.gam - 1.0) / eq.gam; - Ad = Ad*coef; - Dn.set_rows(s,e, Do.rows(s,e)); - - } else if (eq.phys == Equation_mesh) { - double coef = dt*dt*(0.5*eq.gam - eq.beta) / (eq.gam - 1.0); - Dn.set_rows(s,e, Do.rows(s,e) + Yn.rows(s,e)*dt + An.rows(s,e)*coef); - } - } - } else { - Dn.set_rows(s,e, Do.rows(s,e)); - } - } -} - -}; - diff --git a/Code/Source/solver/Integrator.h b/Code/Source/solver/Integrator.h new file mode 100644 index 000000000..486056e01 --- /dev/null +++ b/Code/Source/solver/Integrator.h @@ -0,0 +1,206 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef INTEGRATOR_H +#define INTEGRATOR_H + +#include "Array.h" +#include "Vector.h" +#include "Simulation.h" + +// Note: Solution and SolutionStates structs are defined in ComMod.h +// and available via Simulation.h include chain + +/** + * @brief Integrator class encapsulates the Newton iteration loop for time integration + * + * This class handles the nonlinear Newton iteration scheme for solving coupled + * multi-physics equations in svMultiPhysics. It manages: + * - Solution variables (Ag, Yg, Dg) at generalized-alpha time levels + * - Newton iteration loop with convergence checking + * - Linear system assembly and solve + * - Boundary condition application + * + * Related to GitHub issue #442: Encapsulate the Newton iteration in main.cpp + */ +class Integrator { + +public: + /** + * @brief Construct a new Integrator object + * + * @param simulation Pointer to the Simulation object containing problem data + * @param solutions Solution states containing old time level arrays (takes ownership via move) + */ + Integrator(Simulation* simulation, SolutionStates&& solutions); + + /** + * @brief Destroy the Integrator object + */ + ~Integrator(); + + /** + * @brief Execute one time step with Newton iteration loop + * + * Performs the complete Newton iteration sequence including initialization, + * assembly, boundary condition application, linear solve, and convergence check. + * + * @return True if all equations converged, false otherwise + */ + bool step(); + + /** + * @brief Perform predictor step for next time step + * + * Performs predictor step using generalized-alpha method to estimate + * solution at n+1 time level based on current solution at n time level. + * This should be called once per time step before the Newton iteration loop. + */ + void predictor(); + + /** + * @brief Get reference to solution variable Ag (time derivative of variables) + * + * @return Reference to Ag array (acceleration in structural mechanics) + */ + Array& get_Ag() { return Ag_; } + + /** + * @brief Get reference to solution variable Yg (variables) + * + * @return Reference to Yg array (velocity in structural mechanics) + */ + Array& get_Yg() { return Yg_; } + + /** + * @brief Get reference to solution variable Dg (integrated variables) + * + * @return Reference to Dg array (displacement in structural mechanics) + */ + Array& get_Dg() { return Dg_; } + + /** + * @brief Get reference to solution states struct + * + * Provides access to all solution arrays at old (n) and current (n+1) time levels. + * Use this to access An, Dn, Yn (current) and Ao, Do, Yo (old) via: + * auto& solutions = integrator.get_solutions(); + * auto& An = solutions.current.get_acceleration(); + * auto& Do = solutions.old.get_displacement(); + * + * @return Reference to SolutionStates struct containing all solution arrays + */ + SolutionStates& get_solutions() { return solutions_; } + +private: + /** @brief Pointer to the simulation object */ + Simulation* simulation_; + + /** @brief Time derivative of variables (acceleration in structural mechanics) */ + Array Ag_; + + /** @brief Variables (velocity in structural mechanics) */ + Array Yg_; + + /** @brief Integrated variables (displacement in structural mechanics) */ + Array Dg_; + + /** @brief Solution states at old and current time levels */ + SolutionStates solutions_; + + /** @brief Residual vector for face-based quantities */ + Vector res_; + + /** @brief Increment flag for faces in linear solver */ + Vector incL_; + + /** @brief Newton iteration counter for current time step */ + int newton_count_; + + /** @brief Debug output suffix string combining time step and iteration number */ + std::string istr_; + + /** + * @brief Initialize solution arrays for Ag, Yg, Dg based on problem size + */ + void initialize_arrays(); + + /** + * @brief Perform initiator step for Generalized-alpha Method + * + * Computes quantities at intermediate time levels (n+alpha_m, n+alpha_f) + */ + void initiator_step(); + + /** + * @brief Allocate right-hand side (RHS) and left-hand side (LHS) arrays + * + * @param eq Reference to the equation being solved + */ + void allocate_linear_system(eqType& eq); + + /** + * @brief Set body forces for the current time step + */ + void set_body_forces(); + + /** + * @brief Assemble global equations for all meshes + */ + void assemble_equations(); + + /** + * @brief Apply all boundary conditions (Neumann, Dirichlet, CMM, contact, etc.) + */ + void apply_boundary_conditions(); + + /** + * @brief Solve the assembled linear system + */ + void solve_linear_system(); + + /** + * @brief Perform corrector step and check convergence of all equations + * + * @return True if all equations converged, false otherwise + */ + bool corrector_and_check_convergence(); + + /** + * @brief Update residual and increment arrays for linear solver + * + * @param eq Reference to the equation being solved + */ + void update_residual_arrays(eqType& eq); + + /** + * @brief Initiator function for generalized-alpha method (initiator) + * + * Computes solution variables at intermediate time levels using + * generalized-alpha parameters (am, af) for time integration. + * Updates Ag, Yg, Dg based on An, Ao, Yn, Yo, Dn, Do. + * + * @param Ag Time derivative array at generalized-alpha level + * @param Yg Solution variable array at generalized-alpha level + * @param Dg Integrated variable array at generalized-alpha level + */ + void initiator(Array& Ag, Array& Yg, Array& Dg); + + /** + * @brief Corrector function with convergence check (corrector) + * + * Updates solution at n+1 time level and checks convergence of Newton + * iterations. Also handles equation switching for coupled problems. + */ + void corrector(); + + /** + * @brief Pressure correction for Taylor-Hood elements (corrector_taylor_hood) + * + * Interpolates pressure at edge nodes using reduced basis applied + * on element vertices for Taylor-Hood type elements. + */ + void corrector_taylor_hood(); +}; + +#endif // INTEGRATOR_H diff --git a/Code/Source/solver/Simulation.cpp b/Code/Source/solver/Simulation.cpp index 145ae182e..7e6d56523 100644 --- a/Code/Source/solver/Simulation.cpp +++ b/Code/Source/solver/Simulation.cpp @@ -2,6 +2,7 @@ // SPDX-License-Identifier: BSD-3-Clause #include "Simulation.h" +#include "Integrator.h" #include "all_fun.h" #include "load_msh.h" @@ -9,6 +10,7 @@ #include "mpi.h" #include +#include Simulation::Simulation() { @@ -88,3 +90,25 @@ void Simulation::set_module_parameters() fTmp = general.simulation_initialization_file_path.value(); roInf = general.spectral_radius_of_infinite_time_step.value(); } + +/// @brief Initialize the Integrator object after simulation setup is complete +/// +/// This should be called at the end of initialize() after solution states have been +/// fully initialized. The Integrator takes ownership of these solution states. +/// +/// @param solutions Solution states containing old acceleration, displacement, and velocity +void Simulation::initialize_integrator(SolutionStates&& solutions) +{ + integrator_ = std::make_unique(this, std::move(solutions)); +} + +/// @brief Get reference to the Integrator object +/// +/// @return Reference to the Integrator +Integrator& Simulation::get_integrator() +{ + if (!integrator_) { + throw std::runtime_error("Integrator not initialized. Call initialize_integrator() first."); + } + return *integrator_; +} diff --git a/Code/Source/solver/Simulation.h b/Code/Source/solver/Simulation.h index e0ad1b7b5..680cdddfc 100644 --- a/Code/Source/solver/Simulation.h +++ b/Code/Source/solver/Simulation.h @@ -10,6 +10,10 @@ #include "LinearAlgebra.h" #include +#include + +// Forward declaration +class Integrator; class Simulation { @@ -22,6 +26,11 @@ class Simulation { CepMod& get_cep_mod() { return cep_mod; }; ChnlMod& get_chnl_mod() { return chnl_mod; }; ComMod& get_com_mod() { return com_mod; }; + Integrator& get_integrator(); + + // Initialize the Integrator object after simulation setup is complete + // Takes ownership of solution states via move semantics + void initialize_integrator(SolutionStates&& solutions); // Read a solver paramerer input XML file. void read_parameters(const std::string& fileName); @@ -61,6 +70,10 @@ class Simulation { std::string history_file_name; LinearAlgebra* linear_algebra = nullptr; + + private: + // Time integrator for Newton iteration loop + std::unique_ptr integrator_; }; #endif diff --git a/Code/Source/solver/Vector.h b/Code/Source/solver/Vector.h index 1580447f7..4843ea9b6 100644 --- a/Code/Source/solver/Vector.h +++ b/Code/Source/solver/Vector.h @@ -5,6 +5,7 @@ #define VECTOR_H #include +#include #include #include #include diff --git a/Code/Source/solver/all_fun.cpp b/Code/Source/solver/all_fun.cpp index 6d205544a..d8c9f3846 100644 --- a/Code/Source/solver/all_fun.cpp +++ b/Code/Source/solver/all_fun.cpp @@ -276,7 +276,7 @@ global(const ComMod& com_mod, const CmMod& cm_mod, const mshType& lM, const Arra /// @param s an array containing a scalar value for each node in the mesh /// Replicates 'FUNCTION vIntegM(dId, s, l, u, pFlag)' defined in ALLFUN.f. // -double integ(const ComMod& com_mod, const CmMod& cm_mod, int iM, const Array& s) +double integ(const ComMod& com_mod, const CmMod& cm_mod, int iM, const Array& s, const SolutionStates& solutions) { using namespace consts; @@ -288,6 +288,9 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, int iM, const Array& s, int l, int u, bool pFlag) +double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array& s, int l, int u, + const SolutionStates& solutions, bool pFlag) { using namespace consts; @@ -439,6 +443,9 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array& s, 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) { using namespace consts; #define n_debug_integ_s @@ -689,9 +697,9 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co dmsg << "lFa.iM: " << lFa.iM+1; dmsg << "lFa.name: " << lFa.name; dmsg << "lFa.eType: " << lFa.eType; - #endif + #endif - bool flag = pFlag; + bool flag = pFlag; int nsd = com_mod.nsd; int insd = nsd - 1; @@ -803,7 +811,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, cfg); + nn::gnnb(com_mod, lFa, e, g, nsd, insd, fs.eNoN, Nx, n, solutions, cfg); } // Calculating the Jacobian (encodes area of face element) @@ -841,8 +849,8 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co /// @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, MechanicalConfigurationType cfg) +double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, + const Array& s, const SolutionStates& solutions, MechanicalConfigurationType cfg) { using namespace consts; @@ -923,7 +931,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, cfg); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, solutions, cfg); //CALL GNNB(lFa, e, g, nsd-1, lFa.eNoN, lFa.Nx(:,:,g), n) } else { //CALL GNNIB(lFa, e, g, n) @@ -975,8 +983,8 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, /// @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 int l, std::optional uo, bool THflag, MechanicalConfigurationType cfg) +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) { using namespace consts; @@ -1025,21 +1033,21 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, double result = 0.0; // If s vector, integrate as vector (dot with surface normal) - if (u-l+1 == nsd) { + if (u-l+1 == nsd) { Array vec(nsd,nNo); for (int a = 0; a < nNo; a++) { for (int i = l, n = 0; i <= u; i++, n++) { - vec(n,a) = s(i,a); + vec(n,a) = s(i,a); } } - result = integ(com_mod, cm_mod, lFa, vec, cfg); + result = integ(com_mod, cm_mod, lFa, vec, solutions, cfg); // 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, flag, cfg); + result = integ(com_mod, cm_mod, lFa, sclr, solutions, flag, cfg); } else { throw std::runtime_error("Unexpected dof in integ"); } diff --git a/Code/Source/solver/all_fun.h b/Code/Source/solver/all_fun.h index 4e2b8572f..219327523 100644 --- a/Code/Source/solver/all_fun.h +++ b/Code/Source/solver/all_fun.h @@ -29,18 +29,18 @@ namespace all_fun { Array global(const ComMod& com_mod, const CmMod& cm_mod, const mshType& lM, const Array& U); - double integ(const ComMod& com_mod, const CmMod& cm_mod, int iM, const Array& s); + double integ(const ComMod& com_mod, const CmMod& cm_mod, int iM, const Array& s, const SolutionStates& solutions); - double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array& s, int l, int u, - bool pFlag=false); + 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, - bool pFlag=false, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); + 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, 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 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, 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); 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 b10b3100e..3a66d4cb4 100644 --- a/Code/Source/solver/baf_ini.cpp +++ b/Code/Source/solver/baf_ini.cpp @@ -39,8 +39,13 @@ namespace baf_ini_ns { /// /// Replicates 'SUBROUTINE BAFINI()' defined in BAFINIT.f // -void baf_ini(Simulation* simulation) +void baf_ini(Simulation* simulation, SolutionStates& solutions) { + // Local aliases for solution arrays + const auto& Ao = solutions.old.get_acceleration(); + auto& Do = solutions.old.get_displacement(); + auto& Yo = solutions.old.get_velocity(); + using namespace consts; using namespace fsi_linear_solver; @@ -67,7 +72,7 @@ void baf_ini(Simulation* simulation) continue; } auto& face = msh.fa[iFa]; - face_ini(simulation, msh, face); + face_ini(simulation, msh, face, solutions); } if (msh.lShl) { shl_ini(com_mod, cm_mod, com_mod.msh[iM]); @@ -82,10 +87,10 @@ void baf_ini(Simulation* simulation) auto& bc = eq.bc[iBc]; int iFa = bc.iFa; int iM = bc.iM; - bc_ini(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa]); + bc_ini(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa], solutions); if (com_mod.msh[iM].lShl) { - shl_bc_ini(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa], com_mod.msh[iM]); + shl_bc_ini(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa], com_mod.msh[iM], solutions); } } } @@ -133,7 +138,12 @@ void baf_ini(Simulation* simulation) } if (!com_mod.stFileFlag) { - set_bc::rcr_init(com_mod, cm_mod); + // Create temporary SolutionStates for set_bc calls + SolutionStates temp_solutions; + temp_solutions.old.A = Ao; + temp_solutions.old.D = Do; + temp_solutions.old.Y = Yo; + set_bc::rcr_init(com_mod, cm_mod, temp_solutions); } if (com_mod.cplBC.useGenBC) { @@ -145,7 +155,15 @@ void baf_ini(Simulation* simulation) } if (com_mod.cplBC.schm != CplBCType::cplBC_E) { - set_bc::calc_der_cpl_bc(com_mod, cm_mod); + // Create temporary SolutionStates for set_bc calls + SolutionStates temp_solutions; + temp_solutions.old.A = Ao; + temp_solutions.old.D = Do; + temp_solutions.old.Y = Yo; + temp_solutions.current.A = Yo; + temp_solutions.current.Y = Yo; + temp_solutions.current.D = Do; + set_bc::calc_der_cpl_bc(com_mod, cm_mod, temp_solutions); } } @@ -163,7 +181,7 @@ void baf_ini(Simulation* simulation) int iFa = bc.iFa; int iM = bc.iM; bc.lsPtr = 0; - fsi_ls_ini(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa], lsPtr); + fsi_ls_ini(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa], lsPtr, solutions); } } @@ -225,15 +243,18 @@ void baf_ini(Simulation* simulation) // bc_ini //--------- // -void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& lFa) +void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& lFa, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; using namespace utils; auto& cm = com_mod.cm; int nsd = com_mod.nsd; int tnNo = com_mod.tnNo; - + #define n_debug_bc_ini #ifdef debug_bc_ini DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -272,7 +293,7 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l Vector disp(cm.np()); // Just a constant value for Flat profile - if (btest(lBc.bType, iBC_flat)) { + if (btest(lBc.bType, iBC_flat)) { for (int a = 0; a < lFa.nNo; a++) { int Ac = lFa.gN(a); s(Ac) = 1.0; @@ -284,10 +305,10 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l // 3- maximize ew(i).e where e is the unit vector from current // point to the center 4- Use the point i as the diam here // - } else if (btest(lBc.bType, iBC_para)) { + } 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) / lFa.area; + center(i) = all_fun::integ(com_mod, cm_mod, lFa, com_mod.x, i, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference) / lFa.area; } // gNodes is one if a node located on the boundary (beside iFa) @@ -395,8 +416,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); + 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 (is_zero(tmp)) { tmp = 1.0; throw std::runtime_error("Face '" + lFa.name + "' used for a BC has no non-zero node."); @@ -413,8 +434,11 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l // face_ini //---------- // -void face_ini(Simulation* simulation, mshType& lM, faceType& lFa) +void face_ini(Simulation* simulation, mshType& lM, faceType& lFa, SolutionStates& solutions) { + // Local alias for old displacement + auto& Do = solutions.old.get_displacement(); + using namespace consts; auto& com_mod = simulation->com_mod; auto& cm = com_mod.cm; @@ -436,7 +460,7 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa) // Vector sA(com_mod.tnNo); sA = 1.0; - double area = all_fun::integ(com_mod, cm_mod, lFa, sA); + double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, consts::MechanicalConfigurationType::reference); #ifdef debug_face_ini dmsg << "Face '" << lFa.name << "' area: " << area; #endif @@ -478,7 +502,7 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa) 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); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); @@ -549,7 +573,7 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa) if (com_mod.mvMsh) { for (int i = 0; i < nsd; i++) { - xl(i,a) = xl(i,a) + com_mod.Do(i+nsd+1,Ac); + xl(i,a) = xl(i,a) + Do(i+nsd+1,Ac); } } } @@ -661,8 +685,11 @@ void face_ini(Simulation* simulation, mshType& lM, faceType& lFa) // // Replicates 'SUBROUTINE FSILSINI'. // -void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceType& lFa, int& lsPtr) +void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceType& lFa, int& lsPtr, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; using namespace utils; using namespace fsi_linear_solver; @@ -728,7 +755,7 @@ void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceTyp 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); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, solutions, consts::MechanicalConfigurationType::reference); for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); @@ -868,8 +895,11 @@ void set_shl_xien(Simulation* simulation, mshType& lM) // // Reproduces 'SUBROUTINE SHLBCINI(lBc, lFa, lM)'. // -void shl_bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& lFa, mshType& lM) +void shl_bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& lFa, mshType& lM, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; using namespace utils; diff --git a/Code/Source/solver/baf_ini.h b/Code/Source/solver/baf_ini.h index 5197934e6..fd086b981 100644 --- a/Code/Source/solver/baf_ini.h +++ b/Code/Source/solver/baf_ini.h @@ -9,17 +9,17 @@ namespace baf_ini_ns { -void baf_ini(Simulation* simulation); +void baf_ini(Simulation* simulation, SolutionStates& solutions); -void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& lFa); +void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& lFa, const SolutionStates& solutions); -void face_ini(Simulation* simulation, mshType& lm, faceType& la); +void face_ini(Simulation* simulation, mshType& lm, faceType& la, SolutionStates& solutions); -void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceType& lFa, int& lsPtr); +void fsi_ls_ini(ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, const faceType& lFa, int& lsPtr, const SolutionStates& solutions); void set_shl_xien(Simulation* simulation, mshType& mesh); -void shl_bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& lFa, mshType& lM); +void shl_bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& lFa, mshType& lM, const SolutionStates& solutions); void shl_ini(const ComMod& com_mod, const CmMod& cm_mod, mshType& lM); diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index 2f805e66a..25dd1049f 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -141,8 +141,11 @@ void cep_init_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector& Dg) +void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionStates& solutions) { + // Local aliases for solution arrays + const auto& Dg = solutions.old.get_displacement(); + auto& Yo = solutions.old.get_velocity(); static bool IPASS = true; using namespace consts; @@ -164,7 +167,6 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, const Arra auto& cem = cep_mod.cem; auto& eq = com_mod.eq[iEq]; - auto& Yo = com_mod.Yo; auto& Xion = cep_mod.Xion; int nXion = cep_mod.nXion; diff --git a/Code/Source/solver/cep_ion.h b/Code/Source/solver/cep_ion.h index d5831da89..cca100f3b 100644 --- a/Code/Source/solver/cep_ion.h +++ b/Code/Source/solver/cep_ion.h @@ -19,7 +19,7 @@ void cep_init(Simulation* simulation); void cep_init_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector& X, Vector& Xg); -void cep_integ(Simulation* simulation, const int iEq, const int iDof, const Array& Dg); +void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionStates& solutions); void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector& X, Vector& Xg, const double t1, double& yl, const double I4f, const double dt); diff --git a/Code/Source/solver/cmm.cpp b/Code/Source/solver/cmm.cpp index 1f60588d5..ff7b7092f 100644 --- a/Code/Source/solver/cmm.cpp +++ b/Code/Source/solver/cmm.cpp @@ -260,9 +260,9 @@ void cmm_3d(ComMod& com_mod, const int eNoN, const double w, const Vector& al, const Array& dl, - const Array& xl, const Array& bfl, const Vector& pS0l, const Vector& vwp, - const Vector& ptr) +void cmm_b(ComMod& com_mod, const faceType& lFa, const int e, const Array& al, const Array& dl, + const Array& xl, const Array& bfl, const Vector& pS0l, const Vector& vwp, + const Vector& ptr, const SolutionStates& solutions) { const int nsd = com_mod.nsd; const int dof = com_mod.dof; @@ -282,7 +282,7 @@ 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); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, 3, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); double Jac = sqrt(utils::norm(nV)); nV = nV / Jac; double w = lFa.w(g)*Jac; diff --git a/Code/Source/solver/cmm.h b/Code/Source/solver/cmm.h index 1f1086a7a..b6bb80d56 100644 --- a/Code/Source/solver/cmm.h +++ b/Code/Source/solver/cmm.h @@ -14,9 +14,9 @@ void cmm_3d(ComMod& com_mod, const int eNoN, const double w, const Vector& al, const Array& yl, const Array& bfl, const Array& Kxi, Array& lR, Array3& lK); -void cmm_b(ComMod& com_mod, const faceType& lFa, const int e, const Array& al, const Array& dl, - const Array& xl, const Array& bfl, const Vector& pS0l, const Vector& vwp, - const Vector& ptr); +void cmm_b(ComMod& com_mod, const faceType& lFa, const int e, const Array& al, const Array& dl, + const Array& xl, const Array& bfl, const Vector& pS0l, const Vector& vwp, + const Vector& ptr, const SolutionStates& solutions); void bcmmi(ComMod& com_mod, const int eNoN, const int idof, const double w, const Vector& N, const Array& Nxi, const Array& xl, const Array& tfl, Array& lR); diff --git a/Code/Source/solver/eq_assem.cpp b/Code/Source/solver/eq_assem.cpp index 588737e03..86eaa9760 100644 --- a/Code/Source/solver/eq_assem.cpp +++ b/Code/Source/solver/eq_assem.cpp @@ -28,8 +28,11 @@ namespace eq_assem { -void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& hg, const Array& Yg) +void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& hg, const Array& Yg, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + #define n_debug_b_assem_neu_bc #ifdef debug_b_assem_neu_bc DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -54,9 +57,9 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& cDmn = all_fun::domain(com_mod, msh, cEq, Ec); auto cPhys = eq.dmn[cDmn].phys; - Vector ptr(eNoN); - Vector N(eNoN), hl(eNoN); - Array yl(tDof,eNoN), lR(dof,eNoN); + Vector ptr(eNoN); + Vector N(eNoN), hl(eNoN); + Array yl(tDof,eNoN), lR(dof,eNoN); Array3 lK(dof*dof,eNoN,eNoN); for (int a = 0; a < eNoN; a++) { @@ -76,7 +79,7 @@ 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); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); double Jac = sqrt(utils::norm(nV)); nV = nV / Jac; double w = lFa.w(g)*Jac; @@ -156,13 +159,16 @@ void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& /// @param lFa /// @param hg Pressure magnitude /// @param Dg -void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const Vector& hg, const Array& Dg) +void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const Vector& hg, const Array& Dg, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; using namespace utils; #define n_debug_b_neu_folw_p - #ifdef debug_b_neu_folw_p + #ifdef debug_b_neu_folw_p DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); dmsg << "lFa.name: " << lFa.name; @@ -181,7 +187,7 @@ void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const auto& eq = com_mod.eq[cEq]; auto& cDmn = com_mod.cDmn; - #ifdef debug_b_neu_folw_p + #ifdef debug_b_neu_folw_p dmsg << "nsd: " << nsd; dmsg << "eNoN: " << nsd; #endif @@ -191,13 +197,13 @@ void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const cDmn = all_fun::domain(com_mod, msh, cEq, Ec); // Changes global auto cPhys = eq.dmn[cDmn].phys; - Vector ptr(eNoN); - Vector hl(eNoN); - Array xl(nsd,eNoN); + Vector ptr(eNoN); + Vector hl(eNoN); + Array xl(nsd,eNoN); Array dl(tDof,eNoN); - Vector N(eNoN); - Array Nxi(nsd,eNoN); - Array Nx(nsd,eNoN); + Vector N(eNoN); + Array Nxi(nsd,eNoN); + Array Nx(nsd,eNoN); Array lR(dof,eNoN); Array3 lK(dof*dof,eNoN,eNoN); Array3 lKd; @@ -245,7 +251,7 @@ 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); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoNb, Nx_g, nV, solutions, consts::MechanicalConfigurationType::reference); Jac = sqrt(utils::norm(nV)); nV = nV / Jac; double w = lFa.w(g)*Jac; @@ -286,8 +292,12 @@ void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const /// is eventually used in ADDBCMUL() in the linear solver to add the contribution /// from the resistance BC to the matrix-vector product of the tangent matrix and /// an arbitrary vector. -void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa) +void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa, SolutionStates& solutions) { + // Local aliases for displacement arrays + const auto& Dn = solutions.current.get_displacement(); + const auto& Do = solutions.old.get_displacement(); + using namespace consts; using namespace utils; using namespace fsi_linear_solver; @@ -306,8 +316,8 @@ void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa) int iM = lFa.iM; int nNo = lFa.nNo; - Array sVl(nsd,nNo); - Array sV(nsd,tnNo); + Array sVl(nsd,nNo); + Array sV(nsd,tnNo); // Updating the value of the surface integral of the normal vector // using the deformed configuration ('n' = new = timestep n+1) @@ -322,7 +332,7 @@ void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa) auto cfg = MechanicalConfigurationType::new_timestep; - nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, cfg); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, lFa.eNoN, Nx, n, solutions, cfg); // for (int a = 0; a < lFa.eNoN; a++) { int Ac = lFa.IEN(a,e); @@ -347,9 +357,12 @@ void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa) /// /// Ag(tDof,tnNo), Yg(tDof,tnNo), Dg(tDof,tnNo) // -void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array& Ag, - const Array& Yg, const Array& Dg) +void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array& Ag, + const Array& Yg, const Array& Dg, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + #define n_debug_global_eq_assem #ifdef debug_global_eq_assem DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -407,7 +420,7 @@ void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const break; case EquationType::phys_mesh: - mesh::construct_mesh(com_mod, cep_mod, lM, Ag, Dg); + mesh::construct_mesh(com_mod, cep_mod, lM, Ag, Dg, Do); break; case EquationType::phys_CEP: diff --git a/Code/Source/solver/eq_assem.h b/Code/Source/solver/eq_assem.h index 653daa097..7684c2d6a 100644 --- a/Code/Source/solver/eq_assem.h +++ b/Code/Source/solver/eq_assem.h @@ -9,13 +9,13 @@ namespace eq_assem { -void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& hg, const Array& Yg); +void b_assem_neu_bc(ComMod& com_mod, const faceType& lFa, const Vector& hg, const Array& Yg, const SolutionStates& solutions); -void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const Vector& hg, const Array& Dg); +void b_neu_folw_p(ComMod& com_mod, const bcType& lBc, const faceType& lFa, const Vector& hg, const Array& Dg, const SolutionStates& solutions); -void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa); +void fsi_ls_upd(ComMod& com_mod, const bcType& lBc, const faceType& lFa, SolutionStates& solutions); -void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array& Ag, const Array& Yg, const Array& Dg); +void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array& Ag, const Array& Yg, const Array& Dg, const SolutionStates& solutions); }; diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 7f1e1ae84..8f4cdaad1 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -48,12 +48,18 @@ void finalize(Simulation* simulation) /// /// Reprodices 'SUBROUTINE INITFROMBIN(fName, timeP)' defined in INITIALIZE.f. // -void init_from_bin(Simulation* simulation, const std::string& fName, std::array& timeP) +void init_from_bin(Simulation* simulation, const std::string& fName, std::array& timeP, + SolutionStates& solutions) { + // Local aliases for solution arrays + auto& Ao = solutions.old.get_acceleration(); + auto& Do = solutions.old.get_displacement(); + auto& Yo = solutions.old.get_velocity(); + auto& com_mod = simulation->com_mod; auto& cm = com_mod.cm; int task_id = cm.idcm(); - #ifdef debug_init_from_bin + #ifdef debug_init_from_bin DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); #endif @@ -75,12 +81,12 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< com_mod.timeP[1] = timeP[1]; com_mod.timeP[2] = timeP[2]; - #ifdef debug_init_from_bin + #ifdef debug_init_from_bin dmsg << "dFlag: " << dFlag; dmsg << "sstEq: " << sstEq; dmsg << "pstEq: " << pstEq; dmsg << "cepEq: " << cepEq; - dmsg << "cm.tF(): " << cm.tF(cm_mod); + dmsg << "cm.tF(): " << cm.tF(cm_mod); #endif // Open file and position at the location in the file @@ -95,11 +101,9 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< std::array tStamp; auto& cplBC = com_mod.cplBC; auto& Ad = com_mod.Ad; - auto& Ao = com_mod.Ao; - auto& Do = com_mod.Do; + // Ao, Do, Yo now passed as parameters auto& pS0 = com_mod.pS0; auto& Xion = cep_mod.Xion; - auto& Yo = com_mod.Yo; output::read_restart_header(com_mod, tStamp, timeP[0], bin_file); bin_file.read((char*)cplBC.xo.data(), cplBC.xo.msize()); @@ -260,8 +264,14 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< /// /// Reproduces the Fortran 'INITFROMVTU' subroutine. // -void init_from_vtu(Simulation* simulation, const std::string& fName, std::array& timeP) +void init_from_vtu(Simulation* simulation, const std::string& fName, std::array& timeP, + SolutionStates& solutions) { + // Local aliases for solution arrays + auto& Ao = solutions.old.get_acceleration(); + auto& Do = solutions.old.get_displacement(); + auto& Yo = solutions.old.get_velocity(); + auto& com_mod = simulation->com_mod; auto& cm_mod = simulation->cm_mod; auto& cm = com_mod.cm; @@ -285,17 +295,17 @@ void init_from_vtu(Simulation* simulation, const std::string& fName, std::array< Array tmpA, tmpY, tmpD; if (cm.mas(cm_mod)) { - tmpA.resize(tDof,gtnNo); - tmpY.resize(tDof,gtnNo); + tmpA.resize(tDof,gtnNo); + tmpY.resize(tDof,gtnNo); tmpD.resize(tDof,gtnNo); vtk_xml::read_vtus(simulation, tmpA, tmpY, tmpD, fName); - } else { + } else { //ALLOCATE(tmpA(0,0), tmpY(0,0), tmpD(0,0)) - } + } - com_mod.Ao = all_fun::local(com_mod, cm_mod, cm, tmpA); - com_mod.Yo = all_fun::local(com_mod, cm_mod, cm, tmpY); - com_mod.Do = all_fun::local(com_mod, cm_mod, cm, tmpD); + Ao = all_fun::local(com_mod, cm_mod, cm, tmpA); + Yo = all_fun::local(com_mod, cm_mod, cm, tmpY); + Do = all_fun::local(com_mod, cm_mod, cm, tmpD); } /// @brief Initialize or finalize svFSI variables/structures. @@ -621,13 +631,20 @@ void initialize(Simulation* simulation, Vector& timeP) com_mod.ltg, com_mod.rowPtr, com_mod.colPtr, nFacesLS); // Variable allocation and initialization - int tnNo = com_mod.tnNo; - com_mod.Ao.resize(tDof,tnNo); - com_mod.An.resize(tDof,tnNo); - com_mod.Yo.resize(tDof,tnNo); - com_mod.Yn.resize(tDof,tnNo); - com_mod.Do.resize(tDof,tnNo); - com_mod.Dn.resize(tDof,tnNo); + int tnNo = com_mod.tnNo; + + // Create SolutionStates that will be moved to Integrator at end + SolutionStates initial_solutions; + initial_solutions.old.A.resize(tDof, tnNo); + initial_solutions.old.Y.resize(tDof, tnNo); + initial_solutions.old.D.resize(tDof, tnNo); + // An, Yn, Dn will be created in Integrator class + + // Local aliases for convenience (these reference initial_solutions members) + auto& Ao = initial_solutions.old.get_acceleration(); + auto& Yo = initial_solutions.old.get_velocity(); + auto& Do = initial_solutions.old.get_displacement(); + com_mod.Bf.resize(nsd,tnNo); // [TODO] DaveP not implemented? @@ -688,14 +705,14 @@ void initialize(Simulation* simulation, Vector& timeP) flag = false; } - if (flag) { + if (flag) { auto& iniFilePath = com_mod.iniFilePath; auto& timeP = com_mod.timeP; - if (iniFilePath.find(".bin") != std::string::npos) { - init_from_bin(simulation, iniFilePath, timeP); + if (iniFilePath.find(".bin") != std::string::npos) { + init_from_bin(simulation, iniFilePath, timeP, initial_solutions); } else { - init_from_vtu(simulation, iniFilePath, timeP); + init_from_vtu(simulation, iniFilePath, timeP, initial_solutions); } } else { @@ -705,13 +722,13 @@ void initialize(Simulation* simulation, Vector& timeP) if (FILE *file = fopen(fTmp.c_str(), "r")) { fclose(file); auto& timeP = com_mod.timeP; - init_from_bin(simulation, fTmp, timeP); + init_from_bin(simulation, fTmp, timeP, initial_solutions); } else { if (cm.mas(cm_mod)) { std::cout << "WARNING: No '" + fTmp + "' file to restart simulation from; Resetting restart flag to false"; } com_mod.stFileFlag = false; - zero_init(simulation); + zero_init(simulation, initial_solutions); } if (rmsh.isReqd) { @@ -722,13 +739,13 @@ void initialize(Simulation* simulation, Vector& timeP) for (int i = 0; i < rmsh.iNorm.size(); i++) { rmsh.iNorm(i) = com_mod.eq[i].iNorm; } - rmsh.A0 = com_mod.Ao; - rmsh.Y0 = com_mod.Yo; - rmsh.D0 = com_mod.Do; + rmsh.A0 = Ao; + rmsh.Y0 = Yo; + rmsh.D0 = Do; } } else { - zero_init(simulation); + zero_init(simulation, initial_solutions); } } @@ -741,25 +758,22 @@ void initialize(Simulation* simulation, Vector& timeP) com_mod.eq[i].iNorm = rmsh.iNorm[i]; } - com_mod.Ao = all_fun::local(com_mod, cm_mod, cm, rmsh.A0); - com_mod.Yo = all_fun::local(com_mod, cm_mod, cm, rmsh.Y0); - com_mod.Do = all_fun::local(com_mod, cm_mod, cm, rmsh.D0); + Ao = all_fun::local(com_mod, cm_mod, cm, rmsh.A0); + Yo = all_fun::local(com_mod, cm_mod, cm, rmsh.Y0); + Do = all_fun::local(com_mod, cm_mod, cm, rmsh.D0); - rmsh.A0.resize(tDof,tnNo); - rmsh.A0 = com_mod.Ao; + rmsh.A0.resize(tDof,tnNo); + rmsh.A0 = Ao; - rmsh.Y0.resize(tDof,tnNo); - rmsh.Y0 = com_mod.Yo; + rmsh.Y0.resize(tDof,tnNo); + rmsh.Y0 = Yo; - rmsh.D0.resize(tDof,tnNo); - rmsh.D0 = com_mod.Do; + rmsh.D0.resize(tDof,tnNo); + rmsh.D0 = Do; } // resetSim // Initialize new variables - // - com_mod.An = com_mod.Ao; - com_mod.Yn = com_mod.Yo; - com_mod.Dn = com_mod.Do; + // An, Yn, Dn moved to Integrator class (initialized there from Ao, Yo, Do) for (int iM = 0; iM < nMsh; iM++) { if (cm.mas(cm_mod)) { @@ -808,7 +822,7 @@ void initialize(Simulation* simulation, Vector& timeP) for (int iDmn = 0; iDmn < eq.nDmn; iDmn++) { int i = eq.dmn[iDmn].Id; - eq.dmn[iDmn].v = all_fun::integ(com_mod, cm_mod, i, s, 0, 0); + eq.dmn[iDmn].v = all_fun::integ(com_mod, cm_mod, i, s, 0, 0, initial_solutions, false); if (!com_mod.shlEq && !com_mod.cmmInit) { //std = " Volume of domain <"//STR(i)//"> is "// 2 STR(eq(iEq)%dmn(iDmn)%v) //IF (ISZERO(eq(iEq)%dmn(iDmn)%v)) wrn = "<< Volume of "// "domain "//iDmn//" of equation "//iEq//" is zero >>" @@ -818,7 +832,7 @@ void initialize(Simulation* simulation, Vector& timeP) // Preparing faces and BCs // - baf_ini_ns::baf_ini(simulation); + baf_ini_ns::baf_ini(simulation, initial_solutions); // As all the arrays are allocated, call BIN to VTK for conversion // @@ -830,15 +844,28 @@ void initialize(Simulation* simulation, Vector& timeP) // Making sure the old solution satisfies BCs // - // Modifes - // com_mod.Ao - Old time derivative of variables (acceleration) - // com_mod.Yo - Old variables (velocity) - // com_mod.Do - Old integrated variables (dissplacement) + // Modifies Ao, Yo, Do (local variables) // - set_bc::set_bc_dir(com_mod, com_mod.Ao, com_mod.Yo, com_mod.Do); - - // Preparing TXT files - txt_ns::txt(simulation, true); + SolutionStates temp_solutions; + temp_solutions.old.A = Ao; + temp_solutions.old.Y = Yo; + temp_solutions.old.D = Do; + temp_solutions.current.A = Ao; + temp_solutions.current.Y = Yo; + temp_solutions.current.D = Do; + set_bc::set_bc_dir(com_mod, temp_solutions); + // Copy back modified values + Ao = temp_solutions.current.A; + Yo = temp_solutions.current.Y; + Do = temp_solutions.current.D; + + // Preparing TXT files (pass solution states for initial output) + SolutionStates init_txt_solutions; + init_txt_solutions.current.A = Ao; + init_txt_solutions.current.Y = Yo; + init_txt_solutions.current.D = Do; + init_txt_solutions.old.D = Do; + txt_ns::txt(simulation, true, init_txt_solutions); // Printing the first line and initializing timeP int co = 1; @@ -847,6 +874,10 @@ void initialize(Simulation* simulation, Vector& timeP) std::fill(com_mod.rmsh.flag.begin(), com_mod.rmsh.flag.end(), false); com_mod.resetSim = false; + + // Create Integrator now that initial_solutions (Ao, Do, Yo) are fully initialized + // The Integrator takes ownership via move semantics + simulation->initialize_integrator(std::move(initial_solutions)); } //----------- @@ -854,8 +885,13 @@ void initialize(Simulation* simulation, Vector& timeP) //----------- // Initialize state variables Yo, Ao and Do. // -void zero_init(Simulation* simulation) +void zero_init(Simulation* simulation, SolutionStates& solutions) { + // Local aliases for solution arrays + auto& Ao = solutions.old.get_acceleration(); + auto& Do = solutions.old.get_displacement(); + auto& Yo = solutions.old.get_velocity(); + auto& com_mod = simulation->com_mod; auto& cm = com_mod.cm; const int nsd = com_mod.nsd; @@ -875,7 +911,7 @@ void zero_init(Simulation* simulation) for (int a = 0; a < com_mod.tnNo; a++) { // In the future this should depend on the equation type. for (int i = 0; i < nsd; i++) { - com_mod.Yo(i,a) = msh.Ys(i,a,0); + Yo(i,a) = msh.Ys(i,a,0); } } } @@ -889,7 +925,7 @@ void zero_init(Simulation* simulation) #endif for (int a = 0; a < com_mod.tnNo; a++) { for (int i = 0; i < nsd; i++) { - com_mod.Yo(i,a) = com_mod.Vinit(i,a); + Yo(i,a) = com_mod.Vinit(i,a); } } } @@ -900,7 +936,7 @@ void zero_init(Simulation* simulation) #endif for (int a = 0; a < com_mod.tnNo; a++) { for (int i = 0; i < nsd; i++) { - com_mod.Yo(nsd,a) = com_mod.Pinit(a); + Yo(nsd,a) = com_mod.Pinit(a); } } } @@ -911,7 +947,7 @@ void zero_init(Simulation* simulation) #endif for (int a = 0; a < com_mod.tnNo; a++) { for (int i = 0; i < nsd; i++) { - com_mod.Do(i,a) = com_mod.Dinit(i,a); + Do(i,a) = com_mod.Dinit(i,a); } } } diff --git a/Code/Source/solver/initialize.h b/Code/Source/solver/initialize.h index ae3fe309a..cd65d67cd 100644 --- a/Code/Source/solver/initialize.h +++ b/Code/Source/solver/initialize.h @@ -8,13 +8,15 @@ void finalize(Simulation* simulation); -void init_from_bin(Simulation* simulation, const std::string& fName, std::array& timeP); +void init_from_bin(Simulation* simulation, const std::string& fName, std::array& timeP, + SolutionStates& solutions); -void init_from_vtu(Simulation* simulation, const std::string& fName, std::array& timeP); +void init_from_vtu(Simulation* simulation, const std::string& fName, std::array& timeP, + SolutionStates& solutions); void initialize(Simulation* simulation, Vector& timeP); -void zero_init(Simulation* simulation); +void zero_init(Simulation* simulation, SolutionStates& solutions); #endif diff --git a/Code/Source/solver/main.cpp b/Code/Source/solver/main.cpp index 9a0ab3760..70245f2ee 100644 --- a/Code/Source/solver/main.cpp +++ b/Code/Source/solver/main.cpp @@ -8,6 +8,7 @@ // svMultiPhysics XML_FILE_NAME // #include "Simulation.h" +#include "Integrator.h" #include "all_fun.h" #include "bf.h" @@ -18,7 +19,6 @@ #include "initialize.h" #include "ls.h" #include "output.h" -#include "pic.h" #include "read_files.h" #include "read_msh.h" #include "remesh.h" @@ -85,7 +85,7 @@ void read_files(Simulation* simulation, const std::string& file_name) /// @brief Iterate the precomputed state-variables in time using linear interpolation to the current time step size // -void iterate_precomputed_time(Simulation* simulation) { +void iterate_precomputed_time(Simulation* simulation, Array& An, Array& Yn, Array& Ao, Array& Yo, Array& Do) { using namespace consts; auto& com_mod = simulation->com_mod; @@ -104,14 +104,6 @@ void iterate_precomputed_time(Simulation* simulation) { auto& Rd = com_mod.Rd; // Residual of the displacement equation auto& Kd = com_mod.Kd; // LHS matrix for displacement equation - auto& Ao = com_mod.Ao; // Old time derivative of variables (acceleration) - auto& Yo = com_mod.Yo; // Old variables (velocity) - auto& Do = com_mod.Do; // Old integrated variables (dissplacement) - - auto& An = com_mod.An; // New time derivative of variables - auto& Yn = com_mod.Yn; // New variables (velocity) - auto& Dn = com_mod.Dn; // New integrated variables - int& cTS = com_mod.cTS; int& nITs = com_mod.nITs; double& dt = com_mod.dt; @@ -211,11 +203,8 @@ void iterate_solution(Simulation* simulation) dmsg << "cmmInit: " << com_mod.cmmInit; #endif - Array Ag(tDof,tnNo); - Array Yg(tDof,tnNo); - Array Dg(tDof,tnNo); - Vector res(nFacesLS); - Vector incL(nFacesLS); + // Get Integrator object (created at end of initialize()) + auto& integrator = simulation->get_integrator(); // current time step int& cTS = com_mod.cTS; @@ -243,13 +232,17 @@ void iterate_solution(Simulation* simulation) auto& Rd = com_mod.Rd; // Residual of the displacement equation auto& Kd = com_mod.Kd; // LHS matrix for displacement equation - auto& Ao = com_mod.Ao; // Old time derivative of variables (acceleration) - auto& Yo = com_mod.Yo; // Old variables (velocity) - auto& Do = com_mod.Do; // Old integrated variables (displacement) + // Get reference to solution states from integrator + auto& solutions = integrator.get_solutions(); + + // Local aliases for convenience + auto& Ao = solutions.old.get_acceleration(); // Old time derivative of variables (acceleration) + auto& Yo = solutions.old.get_velocity(); // Old variables (velocity) + auto& Do = solutions.old.get_displacement(); // Old integrated variables (displacement) - auto& An = com_mod.An; // New time derivative of variables (acceleration) - auto& Yn = com_mod.Yn; // New variables (velocity) - auto& Dn = com_mod.Dn; // New integrated variables (displacement) + auto& An = solutions.current.get_acceleration(); // New time derivative of variables (acceleration) + auto& Yn = solutions.current.get_velocity(); // New variables (velocity) + auto& Dn = solutions.current.get_displacement(); // New integrated variables (displacement) bool l1 = false; bool l2 = false; @@ -309,7 +302,7 @@ void iterate_solution(Simulation* simulation) // Compute mesh properties to check if remeshing is required // if (com_mod.mvMsh && com_mod.rmsh.isReqd) { - read_msh_ns::calc_mesh_props(com_mod, cm_mod, com_mod.nMsh, com_mod.msh); + read_msh_ns::calc_mesh_props(com_mod, cm_mod, com_mod.nMsh, com_mod.msh, Do); if (com_mod.resetSim) { #ifdef debug_iterate_solution dmsg << "#### resetSim is true " << std::endl; @@ -323,7 +316,7 @@ void iterate_solution(Simulation* simulation) #ifdef debug_iterate_solution dmsg << "Predictor step ... " << std::endl; #endif - pic::picp(simulation); + integrator.predictor(); // Apply Dirichlet BCs strongly // @@ -337,276 +330,23 @@ void iterate_solution(Simulation* simulation) dmsg << "Apply Dirichlet BCs strongly ..." << std::endl; #endif - set_bc::set_bc_dir(com_mod, An, Yn, Dn); + set_bc::set_bc_dir(com_mod, solutions); if (com_mod.urisFlag) {uris::uris_calc_sdf(com_mod);} - iterate_precomputed_time(simulation); + iterate_precomputed_time(simulation, An, Yn, Ao, Yo, Do); - // Inner loop for Newton iteration + // Inner loop for Newton iteration // - int inner_count = 1; - int reply; - int iEqOld; - - // Looping over Newton iterations - while (true) { - #ifdef debug_iterate_solution - dmsg << "---------- Inner Loop " + std::to_string(inner_count) << " -----------" << std::endl; - dmsg << "cEq: " << cEq; - dmsg << "com_mod.eq[cEq].sym: " << com_mod.eq[cEq].sym; - //simulation->com_mod.timer.set_time(); - #endif - //std::cout << "-------------------------------------" << std::endl; - //std::cout << "inner_count: " << inner_count << std::endl; - - auto istr = "_" + std::to_string(cTS) + "_" + std::to_string(inner_count); - iEqOld = cEq; - auto& eq = com_mod.eq[cEq]; - - if (com_mod.cplBC.coupled && cEq == 0) { - #ifdef debug_iterate_solution - dmsg << "Set coupled BCs " << std::endl; - #endif - set_bc::set_bc_cpl(com_mod, cm_mod); - - set_bc::set_bc_dir(com_mod, An, Yn, Dn); - } - - // Initiator step for Generalized α− Method (quantities at n+am, n+af). - // - // Modifes - // Ag((tDof, tnNo) - - // Yg((tDof, tnNo) - - // Dg((tDof, tnNo) - - // - #ifdef debug_iterate_solution - dmsg << "Initiator step ..." << std::endl; - #endif - pic::pici(simulation, Ag, Yg, Dg); - Ag.write("Ag_pic"+ istr); - Yg.write("Yg_pic"+ istr); - Dg.write("Dg_pic"+ istr); - Yn.write("Yn_pic"+ istr); - - if (Rd.size() != 0) { - Rd = 0.0; - Kd = 0.0; - } - - // Allocate com_mod.R and com_mod.Val arrays. - // - // com_mod.R(dof,tnNo) - // com_mod.Val(dof*dof, com_mod.lhs.nnz) - // - // If Trilinos is used then allocate - // com_mod.tls.W(dof,tnNo) - // com_mod.tls.R(dof,tnNo) - // - #ifdef debug_iterate_solution - dmsg << "Allocating the RHS and LHS" << std::endl; - #endif - - ls_ns::ls_alloc(com_mod, eq); - com_mod.Val.write("Val_alloc"+ istr); - - // Compute body forces. If phys is shells or CMM (init), apply - // contribution from body forces (pressure) to residual - // - // Modifes: com_mod.Bf, Dg - // - #ifdef debug_iterate_solution - dmsg << "Set body forces ..." << std::endl; - #endif - - bf::set_bf(com_mod, Dg); - com_mod.Val.write("Val_bf"+ istr); - - // Assemble equations. - // - #ifdef debug_iterate_solution - dmsg << "Assembling equation: " << eq.sym; - #endif - - for (int iM = 0; iM < com_mod.nMsh; iM++) { - eq_assem::global_eq_assem(com_mod, cep_mod, com_mod.msh[iM], Ag, Yg, Dg); - } - com_mod.R.write("R_as"+ istr); - com_mod.Val.write("Val_as"+ istr); - - // Treatment of boundary conditions on faces - // Apply Neumman or Traction boundary conditions - // - // Modifies: com_mod.R - // - #ifdef debug_iterate_solution - dmsg << "Apply Neumman or Traction BCs ... " << std::endl; - #endif - - Yg.write("Yg_vor_neu"+ istr); - Dg.write("Dg_vor_neu"+ istr); - - set_bc::set_bc_neu(com_mod, cm_mod, Yg, Dg); - - com_mod.Val.write("Val_neu"+ istr); - com_mod.R.write("R_neu"+ istr); - Yg.write("Yg_neu"+ istr); - Dg.write("Dg_neu"+ istr); - - // Apply CMM BC conditions - // - if (!com_mod.cmmInit) { - #ifdef debug_iterate_solution - dmsg << "Apply CMM BC conditions ... " << std::endl; - #endif - set_bc::set_bc_cmm(com_mod, cm_mod, Ag, Dg); - } - - // Apply weakly applied Dirichlet BCs - // - #ifdef debug_iterate_solution - dmsg << "Apply weakly applied Dirichlet BCs ... " << std::endl; - #endif - - set_bc::set_bc_dir_w(com_mod, Yg, Dg); - - if (com_mod.risFlag) { - ris::ris_resbc(com_mod, Yg, Dg); - } - - if (com_mod.ris0DFlag) { - ris::ris0d_bc(com_mod, cm_mod, Yg, Dg); - } - - // Apply contact model and add its contribution to residual - // - if (com_mod.iCntct) { - contact::construct_contact_pnlty(com_mod, cm_mod, Dg); - -#if 0 - if (cTS <= 2050) { - Array::write_enabled = true; - com_mod.R.write("R_"+ std::to_string(cTS)); - //exit(0); - } -#endif - } - - // Synchronize R across processes. Note: that it is important - // to synchronize residual, R before treating immersed bodies as - // ib.R is already communicated across processes - // - if (!eq.assmTLS) { - #ifdef debug_iterate_solution - dmsg << "Synchronize R across processes ..." << std::endl; - #endif - all_fun::commu(com_mod, com_mod.R); - } - - // Update residual in displacement equation for USTRUCT phys. - // Note that this step is done only first iteration. Residual - // will be 0 for subsequent iterations - // - // Modifies com_mod.Rd. - // - #ifdef debug_iterate_solution - dmsg << "com_mod.sstEq: " << com_mod.sstEq; - #endif - if (com_mod.sstEq) { - ustruct::ustruct_r(com_mod, Yg); - } - - // Set the residual of the continuity equation and its tangent matrix - // due to variation with pressure to 0 on all the edge nodes. - // - if (std::set{Equation_stokes, Equation_fluid, Equation_ustruct, Equation_FSI}.count(eq.phys) != 0) { - #ifdef debug_iterate_solution - dmsg << "thood_val_rc ..." << std::endl; - #endif - fs::thood_val_rc(com_mod); - } - - // Treat Neumann boundaries that are not deforming - // - #ifdef debug_iterate_solution - dmsg << "set_bc_undef_neu ..." << std::endl; - #endif - - set_bc::set_bc_undef_neu(com_mod); - - // IB treatment: for explicit coupling, simply construct residual. - // - /* [NOTE] not implemented. - if (com_mod.ibFlag) { - if (com_mod.ib.cpld == ibCpld_I) { - //CALL IB_IMPLICIT(Ag, Yg, Dg) - } - // CALL IB_CONSTRUCT() - } - */ - - #ifdef debug_iterate_solution - dmsg << "Update res() and incL ..." << std::endl; - dmsg << "nFacesLS: " << nFacesLS; - #endif - incL = 0; - if (eq.phys == Equation_mesh) { - incL(nFacesLS-1) = 1; - } - - if (com_mod.cmmInit) { - incL(nFacesLS-1) = 1; - } - - for (int iBc = 0; iBc < eq.nBc; iBc++) { - int i = eq.bc[iBc].lsPtr; - if (i != -1) { - // Resistance term for coupled Neumann BC tangent contribution - res(i) = eq.gam * dt * eq.bc[iBc].r; - incL(i) = 1; - } - } - - // Solve equation. - // - // Modifies: com_mod.R, com_mod.Val - // - #ifdef debug_iterate_solution - dmsg << "Solving equation: " << eq.sym; - #endif - - ls_ns::ls_solve(com_mod, eq, incL, res); - - com_mod.Val.write("Val_solve"+ istr); - com_mod.R.write("R_solve"+ istr); - - // Solution is obtained, now updating (Corrector) and check for convergence - // - // Modifies: com_mod.An com_mod.Dn com_mod.Yn cep_mod.Xion com_mod.pS0 com_mod.pSa - // com_mod.pSn com_mod.cEq eq.FSILS.RI.iNorm eq.pNorm - // - #ifdef debug_iterate_solution - dmsg << "Update corrector ..." << std::endl; - #endif - - pic::picc(simulation); - com_mod.Yn.write("Yn_picc"+ istr); - - // Writing out the time passed, residual, and etc. - if (std::count_if(com_mod.eq.begin(),com_mod.eq.end(),[](eqType& eq){return eq.ok;}) == com_mod.eq.size()) { - #ifdef debug_iterate_solution - dmsg << ">>> All OK" << std::endl; - dmsg << "iEqOld: " << iEqOld+1; - #endif - break; - } - output::output_result(simulation, com_mod.timeP, 2, iEqOld); + #ifdef debug_iterate_solution + dmsg << "Starting Newton iteration via Integrator ..." << std::endl; + #endif - inner_count += 1; - } // Inner loop + int iEqOld = cEq; + integrator.step(); #ifdef debug_iterate_solution - dmsg << ">>> End of inner loop " << std::endl; + dmsg << ">>> End of Newton iteration" << std::endl; #endif // IB treatment: interpolate flow data on IB mesh from background @@ -624,7 +364,7 @@ void iterate_solution(Simulation* simulation) */ if (com_mod.risFlag) { - ris::ris_meanq(com_mod, cm_mod); + ris::ris_meanq(com_mod, cm_mod, solutions); ris::ris_status(com_mod, cm_mod); if (cm.mas(cm_mod)) { std::cout << "Iteration: " << com_mod.cTS << std::endl; @@ -642,7 +382,7 @@ void iterate_solution(Simulation* simulation) std::cout << "Valve status just changed. Do not update" << std::endl; } } else { - ris::ris_updater(com_mod, cm_mod); + ris::ris_updater(com_mod, cm_mod, solutions); } // goto label_11; } @@ -654,7 +394,7 @@ void iterate_solution(Simulation* simulation) dmsg << "Saving the TXT files containing ECGs ..." << std::endl; #endif - txt_ns::txt(simulation, false); + txt_ns::txt(simulation, false, solutions); // If remeshing is required then save current solution. // @@ -670,9 +410,9 @@ void iterate_solution(Simulation* simulation) com_mod.rmsh.iNorm(i) = com_mod.eq[i].iNorm; } - com_mod.rmsh.A0 = com_mod.Ao; - com_mod.rmsh.Y0 = com_mod.Yo; - com_mod.rmsh.D0 = com_mod.Do; + com_mod.rmsh.A0 = Ao; + com_mod.rmsh.Y0 = Yo; + com_mod.rmsh.D0 = Do; } } @@ -718,7 +458,7 @@ void iterate_solution(Simulation* simulation) // Saving the result to restart bin file if (l1 || l2) { - output::write_restart(simulation, com_mod.timeP); + output::write_restart(simulation, com_mod.timeP, solutions); } // Writing results into the disk with VTU format @@ -760,7 +500,7 @@ void iterate_solution(Simulation* simulation) // [HZ] Part related to RIS0D if (cEq == 0 && com_mod.ris0DFlag) { - ris::ris0d_status(com_mod, cm_mod); + ris::ris0d_status(com_mod, cm_mod, solutions); } // [HZ] Part related to unfitted RIS @@ -769,12 +509,12 @@ void iterate_solution(Simulation* simulation) for (int iUris = 0; iUris < com_mod.nUris; iUris++) { com_mod.uris[iUris].cnt++; if (com_mod.uris[iUris].clsFlg) { - uris::uris_meanp(com_mod, cm_mod, iUris); + uris::uris_meanp(com_mod, cm_mod, iUris, solutions); // if (com_mod.uris[iUris].cnt == 1) { // // GOTO 11 // The GOTO Statement in the Fortran code // } } else { - uris::uris_meanv(com_mod, cm_mod, iUris); + uris::uris_meanv(com_mod, cm_mod, iUris, solutions); } if (cm.mas(cm_mod)) { std::cout << " URIS surface: " << com_mod.uris[iUris].name << ", count: " << com_mod.uris[iUris].cnt << std::endl; @@ -782,7 +522,7 @@ void iterate_solution(Simulation* simulation) } if (com_mod.mvMsh) { - uris::uris_update_disp(com_mod, cm_mod); + uris::uris_update_disp(com_mod, cm_mod, solutions); } if (cm.mas(cm_mod)) { diff --git a/Code/Source/solver/mesh.cpp b/Code/Source/solver/mesh.cpp index 6a4ac44b3..8e754c088 100644 --- a/Code/Source/solver/mesh.cpp +++ b/Code/Source/solver/mesh.cpp @@ -19,9 +19,9 @@ namespace mesh { -void construct_mesh(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array& Ag, const Array& Dg) +void construct_mesh(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array& Ag, const Array& Dg, const Array& Do) { - #define n_debug_construct_mesh + #define n_debug_construct_mesh #ifdef debug_construct_mesh DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); @@ -37,7 +37,6 @@ void construct_mesh(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const A const auto& eq = com_mod.eq[cEq]; auto& cDmn = com_mod.cDmn; const int nsymd = com_mod.nsymd; - auto& Do = com_mod.Do; auto& pS0 = com_mod.pS0; auto& pSn = com_mod.pSn; auto& pSa = com_mod.pSa; diff --git a/Code/Source/solver/mesh.h b/Code/Source/solver/mesh.h index f6bae13f3..7df11f0be 100644 --- a/Code/Source/solver/mesh.h +++ b/Code/Source/solver/mesh.h @@ -10,7 +10,7 @@ namespace mesh { -void construct_mesh(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array& Ag, const Array& Dg); +void construct_mesh(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array& Ag, const Array& Dg, const Array& Do); }; diff --git a/Code/Source/solver/nn.cpp b/Code/Source/solver/nn.cpp index 6a78ca100..f842b345e 100644 --- a/Code/Source/solver/nn.cpp +++ b/Code/Source/solver/nn.cpp @@ -522,9 +522,12 @@ void gnn(const int eNoN, const int nsd, const int insd, Array& Nxi, Arra /// /// 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, MechanicalConfigurationType cfg) +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) { + // Local aliases for displacement arrays + const auto& Dn = solutions.current.get_displacement(); + const auto& Do = solutions.old.get_displacement(); auto& cm = com_mod.cm; #define n_debug_gnnb @@ -608,7 +611,7 @@ 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 - lX(i,a) = lX(i,a) + com_mod.Do(i+nsd+1,Ac); + lX(i,a) = lX(i,a) + Do(i+nsd+1,Ac); } } else { @@ -619,13 +622,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) + com_mod.Do(i,Ac); + lX(i,a) = lX(i,a) + Do(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) + com_mod.Dn(i,Ac); + lX(i,a) = lX(i,a) + Dn(i,Ac); } break; default: diff --git a/Code/Source/solver/nn.h b/Code/Source/solver/nn.h index 758ac491a..123e8ef31 100644 --- a/Code/Source/solver/nn.h +++ b/Code/Source/solver/nn.h @@ -38,7 +38,7 @@ namespace nn { 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, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); + const int eNoNb, const Array& Nx, Vector& n, const SolutionStates& solutions, consts::MechanicalConfigurationType cfg=consts::MechanicalConfigurationType::reference); 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/output.cpp b/Code/Source/solver/output.cpp index cd575ed4a..cba75d51c 100644 --- a/Code/Source/solver/output.cpp +++ b/Code/Source/solver/output.cpp @@ -172,8 +172,13 @@ void read_restart_header(ComMod& com_mod, std::array& tStamp, double& tim /// @brief Reproduces the Fortran 'WRITERESTART' subroutine. // -void write_restart(Simulation* simulation, std::array& timeP) +void write_restart(Simulation* simulation, std::array& timeP, const SolutionStates& solutions) { + // Local aliases for solution arrays + auto& An = solutions.current.get_acceleration(); + auto& Yn = solutions.current.get_velocity(); + auto& Dn = solutions.current.get_displacement(); + auto& com_mod = simulation->com_mod; #define n_debug_write_restart #ifdef debug_write_restart @@ -193,7 +198,7 @@ void write_restart(Simulation* simulation, std::array& timeP) const bool ibFlag = com_mod.ibFlag; const bool dFlag = com_mod.dFlag; - const bool sstEq = com_mod.sstEq; + const bool sstEq = com_mod.sstEq; const bool pstEq = com_mod.pstEq; const bool cepEq = cep_mod.cepEq; const bool risFlag = com_mod.risFlag; @@ -202,10 +207,7 @@ void write_restart(Simulation* simulation, std::array& timeP) auto& cplBC = com_mod.cplBC; auto& Ad = com_mod.Ad; - auto& An = com_mod.An; - auto& Dn = com_mod.Dn; auto& pS0 = com_mod.pS0; - auto& Yn = com_mod.Yn; auto& Xion = cep_mod.Xion; auto& cem = cep_mod.cem; @@ -395,13 +397,14 @@ void write_restart_header(ComMod& com_mod, std::array& timeP, std::ofs /// /// Reproduces: WRITE(fid, REC=myID) stamp, cTS, time,CPUT()-timeP(1), eq.iNorm, cplBC.xn, Yn, An, Dn // -void write_results(ComMod& com_mod, const std::array& timeP, const std::string& fName, const bool sstEq) +void write_results(ComMod& com_mod, const std::array& timeP, const std::string& fName, const bool sstEq, const SolutionStates& solutions) { - int cTS = com_mod.cTS; + // Local aliases for solution arrays + auto& An = solutions.current.get_acceleration(); + auto& Yn = solutions.current.get_velocity(); + auto& Dn = solutions.current.get_displacement(); - auto& An = com_mod.An; - auto& Dn = com_mod.Dn; - auto& Yn = com_mod.Yn; + int cTS = com_mod.cTS; auto& stamp = com_mod.stamp; diff --git a/Code/Source/solver/output.h b/Code/Source/solver/output.h index f0285d948..737678964 100644 --- a/Code/Source/solver/output.h +++ b/Code/Source/solver/output.h @@ -15,11 +15,11 @@ void output_result(Simulation* simulation, std::array& timeP, const i void read_restart_header(ComMod& com_mod, std::array& tStamp, double& timeP, std::ifstream& restart_file); -void write_restart(Simulation* simulation, std::array& timeP); +void write_restart(Simulation* simulation, std::array& timeP, const SolutionStates& solutions); void write_restart_header(ComMod& com_mod, std::array& timeP, std::ofstream& restart_file); -void write_results(ComMod& com_mod, const std::array& timeP, const std::string& fName, const bool sstEq); +void write_results(ComMod& com_mod, const std::array& timeP, const std::string& fName, const bool sstEq, const SolutionStates& solutions); }; diff --git a/Code/Source/solver/pic.h b/Code/Source/solver/pic.h deleted file mode 100644 index af7aaf397..000000000 --- a/Code/Source/solver/pic.h +++ /dev/null @@ -1,22 +0,0 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. -// SPDX-License-Identifier: BSD-3-Clause - -#include "Simulation.h" - -#ifndef PIC_H -#define PIC_H - -namespace pic { - -void picc(Simulation* simulation); - -void pic_eth(Simulation* simulation); - -void pici(Simulation* simulation, Array& Ag, Array& Yg, Array& Dg); - -void picp(Simulation* simulation); - -}; - -#endif - diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 4721b673a..bb0a6c166 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -1164,12 +1164,20 @@ void ppbin2vtk(Simulation* simulation) continue; } + // Create local SolutionStates for loading bin file + const int tDof = com_mod.tDof; + const int tnNo = com_mod.tnNo; + SolutionStates temp_solutions; + temp_solutions.old.A.resize(tDof, tnNo); + temp_solutions.old.Y.resize(tDof, tnNo); + temp_solutions.old.D.resize(tDof, tnNo); + std::array rtmp; - init_from_bin(simulation, fName, rtmp); + init_from_bin(simulation, fName, rtmp, temp_solutions); bool lAve = false; - vtk_xml::write_vtus(simulation, com_mod.Ao, com_mod.Yo, com_mod.Do, lAve); + vtk_xml::write_vtus(simulation, temp_solutions.old.A, temp_solutions.old.Y, temp_solutions.old.D, lAve); } } diff --git a/Code/Source/solver/read_msh.cpp b/Code/Source/solver/read_msh.cpp index 0fdb19000..18d062d85 100644 --- a/Code/Source/solver/read_msh.cpp +++ b/Code/Source/solver/read_msh.cpp @@ -52,9 +52,9 @@ std::map> check_element_conn /// @brief Calculate element Aspect Ratio of a given mesh // -void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag) +void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const Array& Do) { - #define n_debug_calc_elem_ar + #define n_debug_calc_elem_ar #ifdef debug_calc_elem_ar DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); @@ -65,7 +65,7 @@ void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag if (lM.eType != ElementType::TET4 && lM.eType != ElementType::TRI3) { // wrn = "AR is computed for TRI and TET elements only" - return; + return; } const int nsd = com_mod.nsd; @@ -85,7 +85,7 @@ void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag xl.set_col(a, com_mod.x.col(Ac)); if (com_mod.mvMsh) { for (int i = 0; i < nsd; i++) { - dol(i,a) = com_mod.Do(i+nsd+1,Ac); + dol(i,a) = Do(i+nsd+1,Ac); } } } @@ -147,10 +147,10 @@ void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag /// @brief Calculate element Jacobian of a given mesh. // -void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag) +void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const Array& Do) { - #define n_debug_calc_elem_jac - #ifdef debug_calc_elem_jac + #define n_debug_calc_elem_jac + #ifdef debug_calc_elem_jac DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); dmsg << "lM.nEl: " << lM.nEl; @@ -161,7 +161,7 @@ void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rfla #endif using namespace consts; const int nsd = com_mod.nsd; - rflag = false; + rflag = false; // Careful here, lM.nEl can be 0. // @@ -182,8 +182,8 @@ void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rfla xl.set_col(a, com_mod.x.col(Ac)); if (com_mod.mvMsh) { for (int i = 0; i < nsd; i++) { - dol(i,a) = com_mod.Do(i+nsd+1,Ac); - //dol(i,a) = com_mod.Do(i+nsd,Ac); + dol(i,a) = Do(i+nsd+1,Ac); + //dol(i,a) = Do(i+nsd,Ac); } } } @@ -267,7 +267,7 @@ void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rfla /// @brief Calculate element Skewness of a given mesh. // -void calc_elem_skew(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag) +void calc_elem_skew(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const Array& Do) { #define n_debug_calc_elem_skew #ifdef debug_calc_elem_skew @@ -299,7 +299,7 @@ void calc_elem_skew(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rfl xl.set_col(a, com_mod.x.col(Ac)); if (com_mod.mvMsh) { for (int i = 0; i < nsd; i++) { - dol(i,a) = com_mod.Do(i+nsd+1,Ac); + dol(i,a) = Do(i+nsd+1,Ac); } } } @@ -355,7 +355,7 @@ void calc_elem_skew(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rfl #endif } -void calc_mesh_props(ComMod& com_mod, const CmMod& cm_mod, const int nMesh, std::vector& mesh) +void calc_mesh_props(ComMod& com_mod, const CmMod& cm_mod, const int nMesh, std::vector& mesh, const Array& Do) { #define n_debug_calc_mesh_props #ifdef debug_calc_mesh_props @@ -366,11 +366,11 @@ void calc_mesh_props(ComMod& com_mod, const CmMod& cm_mod, const int nMesh, std: using namespace consts; auto& rmsh = com_mod.rmsh; #ifdef debug_calc_mesh_props - dmsg << "nMesh: " << nMesh; - dmsg << "resetSim: " << com_mod.resetSim; + dmsg << "nMesh: " << nMesh; + dmsg << "resetSim: " << com_mod.resetSim; dmsg << "com_mod.cTS: " << com_mod.cTS; - dmsg << "rmsh.fTS: " << rmsh.fTS; - dmsg << "rmsh.freq: " << rmsh.freq; + dmsg << "rmsh.fTS: " << rmsh.fTS; + dmsg << "rmsh.freq: " << rmsh.freq; #endif for (int iM = 0; iM < nMesh; iM++) { @@ -378,9 +378,9 @@ void calc_mesh_props(ComMod& com_mod, const CmMod& cm_mod, const int nMesh, std: dmsg << "----- mesh " + mesh[iM].name << " -----"; #endif bool flag = false; - calc_elem_jac(com_mod, cm_mod, mesh[iM], flag); - calc_elem_skew(com_mod, cm_mod, mesh[iM], flag); - calc_elem_ar(com_mod, cm_mod, mesh[iM], flag); + calc_elem_jac(com_mod, cm_mod, mesh[iM], flag, Do); + calc_elem_skew(com_mod, cm_mod, mesh[iM], flag, Do); + calc_elem_ar(com_mod, cm_mod, mesh[iM], flag, Do); rmsh.flag[iM] = flag; #ifdef debug_calc_mesh_props dmsg << "mesh[iM].flag: " << rmsh.flag[iM]; diff --git a/Code/Source/solver/read_msh.h b/Code/Source/solver/read_msh.h index 600f174fd..bff955768 100644 --- a/Code/Source/solver/read_msh.h +++ b/Code/Source/solver/read_msh.h @@ -21,11 +21,11 @@ namespace read_msh_ns { Vector gN; }; - void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag); - void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag); - void calc_elem_skew(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag); + void calc_elem_ar(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const Array& Do); + void calc_elem_jac(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const Array& Do); + void calc_elem_skew(ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag, const Array& Do); - void calc_mesh_props(ComMod& com_mod, const CmMod& cm_mod, const int nMesh, std::vector& mesh); + void calc_mesh_props(ComMod& com_mod, const CmMod& cm_mod, const int nMesh, std::vector& mesh, const Array& Do); void calc_nbc(mshType& mesh, faceType& face); diff --git a/Code/Source/solver/remesh.cpp b/Code/Source/solver/remesh.cpp index c8db1ed1c..d557ec520 100644 --- a/Code/Source/solver/remesh.cpp +++ b/Code/Source/solver/remesh.cpp @@ -1888,14 +1888,8 @@ void remesh_restart(Simulation* simulation) com_mod.idMap.clear(); com_mod.cmmBdry.clear(); com_mod.iblank.clear(); - com_mod.Ao.clear(); - com_mod.An.clear(); - com_mod.Do.clear(); - com_mod.Dn.clear(); com_mod.R.clear(); com_mod.Val.clear(); - com_mod.Yo.clear(); - com_mod.Yn.clear(); com_mod.Bf.clear(); cplBC.nFa = 0; diff --git a/Code/Source/solver/ris.cpp b/Code/Source/solver/ris.cpp index d55a768d9..d53ff688f 100644 --- a/Code/Source/solver/ris.cpp +++ b/Code/Source/solver/ris.cpp @@ -12,9 +12,15 @@ namespace ris { -/// @brief This subroutine computes the mean pressure and flux on the ris surface -void ris_meanq(ComMod& com_mod, CmMod& cm_mod) +/// @brief This subroutine computes the mean pressure and flux on the ris surface +void ris_meanq(ComMod& com_mod, CmMod& cm_mod, SolutionStates& solutions) { + // Local aliases for solution arrays + auto& An = solutions.current.get_acceleration(); + auto& Yn = solutions.current.get_velocity(); + auto& Dn = solutions.current.get_displacement(); + const auto& Do = solutions.old.get_displacement(); + #define n_debug_ris_meanq #ifdef debug_ris_meanq DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -31,11 +37,7 @@ void ris_meanq(ComMod& com_mod, CmMod& cm_mod) const int nsd = com_mod.nsd; const int cEq = com_mod.cEq; - - auto& An = com_mod.An; auto& Ad = com_mod.Ad; - auto& Dn = com_mod.Dn; - auto& Yn = com_mod.Yn; Array tmpV(maxNSD, com_mod.tnNo); @@ -58,17 +60,17 @@ void ris_meanq(ComMod& com_mod, CmMod& cm_mod) 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)/tmp; + 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; } } // For the velocity - m = nsd; + m = nsd; s = eq[iEq].s; e = s + m - 1; for (int iProj = 0; iProj < nPrj; iProj++) { - // tmpV[0:m,:] = Yn[s:e,:]; + // tmpV[0:m,:] = Yn[s:e,:]; for (int i = 0; i < m; i++) { for (int j = 0; j < Yn.ncols(); j++) { tmpV(i,j) = Yn(s+i,j); @@ -76,11 +78,11 @@ void ris_meanq(ComMod& com_mod, CmMod& cm_mod) } 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, m-1); + RIS.meanFl(iProj) = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference); if (cm.mas(cm_mod)) { std::cout << "For RIS projection: " << iProj << std::endl; - std::cout << " The average pressure is: " << RIS.meanP(iProj,0) << ", " + std::cout << " The average pressure is: " << RIS.meanP(iProj,0) << ", " << RIS.meanP(iProj,1) << std::endl; std::cout << " The average flow is: " << RIS.meanFl(iProj) << std::endl; } @@ -88,8 +90,11 @@ void ris_meanq(ComMod& com_mod, CmMod& cm_mod) } /// @brief Weak treatment of RIS resistance boundary conditions -void ris_resbc(ComMod& com_mod, const Array& Yg, const Array& Dg) +void ris_resbc(ComMod& com_mod, const Array& Yg, const Array& Dg, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; #define n_debug_ris_resbc #ifdef debug_ris_resbc @@ -136,7 +141,7 @@ void ris_resbc(ComMod& com_mod, const Array& Yg, const Array& Dg if (cPhys == EquationType::phys_fluid) { // Build the correct BC - set_bc::set_bc_dir_wl(com_mod, lBc, msh[iM], msh[iM].fa[iFa], Yg, Dg); + set_bc::set_bc_dir_wl(com_mod, lBc, msh[iM], msh[iM].fa[iFa], Yg, Dg, solutions); } lBc.gx.clear(); } @@ -145,17 +150,25 @@ void ris_resbc(ComMod& com_mod, const Array& Yg, const Array& Dg } -void setbc_ris(ComMod& com_mod, const bcType& lBc, const mshType& lM, const faceType& lFa, - const Array& Yg, const Array& Dg) +void setbc_ris(ComMod& com_mod, const bcType& lBc, const mshType& lM, const faceType& lFa, + const Array& Yg, const Array& Dg, const SolutionStates& solutions) { // [HZ] looks not needed in the current implementation } -/// @brief This subroutine updates the resistance and activation flag for the -/// closed and open configurations of the RIS surfaces -void ris_updater(ComMod& com_mod, CmMod& cm_mod) +/// @brief This subroutine updates the resistance and activation flag for the +/// closed and open configurations of the RIS surfaces +void ris_updater(ComMod& com_mod, CmMod& cm_mod, SolutionStates& solutions) { + // Local aliases for solution arrays + auto& An = solutions.current.get_acceleration(); + auto& Yn = solutions.current.get_velocity(); + auto& Dn = solutions.current.get_displacement(); + auto& Ao = solutions.old.get_acceleration(); + auto& Yo = solutions.old.get_velocity(); + auto& Do = solutions.old.get_displacement(); + #define n_debug_ris_updater #ifdef debug_ris_updater DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -180,13 +193,13 @@ void ris_updater(ComMod& com_mod, CmMod& cm_mod) std::cout << "RIS Proj " << iProj << ": Going from close to open." << std::endl; } RIS.nbrIter(iProj) = 0; - // I needed to update the state variables when the valve + // I needed to update the state variables when the valve // goes from close to open to prevent the valve goes back // to close at the next iteration. This is needed only for // close to open and cannot be used for open to close. - com_mod.Ao = com_mod.An; - com_mod.Yo = com_mod.Yn; - if (com_mod.dFlag) {com_mod.Do = com_mod.Dn;} + Ao = An; + Yo = Yn; + if (com_mod.dFlag) {Do = Dn;} com_mod.cplBC.xo = com_mod.cplBC.xn; } } else { @@ -347,14 +360,18 @@ void clean_r_ris(ComMod& com_mod) // [HZ] looks not needed in the current implementation } -void setbcdir_ris(ComMod& com_mod, Array& lA, Array& lY, Array& lD) +void setbcdir_ris(ComMod& com_mod, SolutionStates& solutions) { // [HZ] looks not needed in the current implementation } /// RIS0D code -void ris0d_bc(ComMod& com_mod, CmMod& cm_mod, const Array& Yg, const Array& Dg) +void ris0d_bc(ComMod& com_mod, CmMod& cm_mod, const Array& Yg, const Array& Dg, SolutionStates& solutions) { + // Local aliases for solution arrays + auto& Yn = solutions.current.get_velocity(); + const auto& Do = solutions.old.get_displacement(); + using namespace consts; #define n_debug_ris0d_bc @@ -390,23 +407,29 @@ void ris0d_bc(ComMod& com_mod, CmMod& cm_mod, const Array& Yg, const Arr lBc.eDrn.resize(nsd); lBc.eDrn = 0; - // Apply bc Dir + // Apply bc Dir lBc.gx.resize(msh[iM].fa[iFa].nNo); lBc.gx = 1.0; - set_bc::set_bc_dir_wl(com_mod, lBc, msh[iM], msh[iM].fa[iFa], Yg, Dg); + set_bc::set_bc_dir_wl(com_mod, lBc, msh[iM], msh[iM].fa[iFa], Yg, Dg, solutions); lBc.gx.clear(); lBc.eDrn.clear(); } else { - // Apply Neu bc - set_bc::set_bc_neu_l(com_mod, cm_mod, eq[cEq].bc[iBc], msh[iM].fa[iFa], Yg, Dg); + // Apply Neu bc + set_bc::set_bc_neu_l(com_mod, cm_mod, eq[cEq].bc[iBc], msh[iM].fa[iFa], Yg, Dg, solutions); } } } -void ris0d_status(ComMod& com_mod, CmMod& cm_mod)//, const Array& Yg, const Array& Dg) +void ris0d_status(ComMod& com_mod, CmMod& cm_mod, SolutionStates& solutions) { + // Local aliases for solution arrays + auto& An = solutions.current.get_acceleration(); + auto& Yn = solutions.current.get_velocity(); + auto& Dn = solutions.current.get_displacement(); + const auto& Do = solutions.old.get_displacement(); + using namespace consts; #define n_debug_status @@ -422,11 +445,7 @@ void ris0d_status(ComMod& com_mod, CmMod& cm_mod)//, const Array& Yg, co const int nsd = com_mod.nsd; const int cEq = com_mod.cEq; - - auto& An = com_mod.An; auto& Ad = com_mod.Ad; - auto& Dn = com_mod.Dn; - auto& Yn = com_mod.Yn; bcType lBc; faceType lFa; @@ -460,8 +479,8 @@ void ris0d_status(ComMod& com_mod, CmMod& cm_mod)//, const Array& Yg, co 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); - meanP = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, m-1)/tmp_new; + 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; // For the velocity m = nsd; @@ -476,7 +495,7 @@ void ris0d_status(ComMod& com_mod, CmMod& cm_mod)//, const Array& Yg, co } } - meanFl = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, m-1); + meanFl = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference); 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/ris.h b/Code/Source/solver/ris.h index 5d55d3db1..26bb88ac1 100644 --- a/Code/Source/solver/ris.h +++ b/Code/Source/solver/ris.h @@ -8,12 +8,12 @@ namespace ris { -void ris_meanq(ComMod& com_mod, CmMod& cm_mod); -void ris_resbc(ComMod& com_mod, const Array& Yg, const Array& Dg); -void setbc_ris(ComMod& com_mod, const bcType& lBc, const mshType& lM, const faceType& lFa, - const Array& Yg, const Array& Dg); +void ris_meanq(ComMod& com_mod, CmMod& cm_mod, SolutionStates& solutions); +void ris_resbc(ComMod& com_mod, const Array& Yg, const Array& Dg, const SolutionStates& solutions); +void setbc_ris(ComMod& com_mod, const bcType& lBc, const mshType& lM, const faceType& lFa, + const Array& Yg, const Array& Dg, const SolutionStates& solutions); -void ris_updater(ComMod& com_mod, CmMod& cm_mod); +void ris_updater(ComMod& com_mod, CmMod& cm_mod, SolutionStates& solutions); void ris_status(ComMod& com_mod, CmMod& cm_mod); void doassem_ris(ComMod& com_mod, const int d, const Vector& eqN, @@ -23,11 +23,11 @@ void doassem_velris(ComMod& com_mod, const int d, const Array& eqN, const Array3& lK, const Array& lR); void clean_r_ris(ComMod& com_mod); -void setbcdir_ris(ComMod& com_mod, Array& lA, Array& lY, Array& lD); +void setbcdir_ris(ComMod& com_mod, SolutionStates& solutions); // TODO: RIS 0D code -void ris0d_bc(ComMod& com_mod, CmMod& cm_mod, const Array& Yg, const Array& Dg); -void ris0d_status(ComMod& com_mod, CmMod& cm_mod); //, const Array& Yg, const Array& Dg); +void ris0d_bc(ComMod& com_mod, CmMod& cm_mod, const Array& Yg, const Array& Dg, SolutionStates& solutions); +void ris0d_status(ComMod& com_mod, CmMod& cm_mod, SolutionStates& solutions); }; diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 254229874..720a75e4e 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -25,10 +25,18 @@ namespace set_bc { /// as well as the resistance matrix M ~ dP/dQ from 0D using finite difference. /// Updates the pressure or flowrates stored in cplBC.fa[i].y and the resistance /// matrix M ~ dP/dQ stored in eq.bc[iBc].r. -/// @param com_mod -/// @param cm_mod -void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) +/// @param com_mod +/// @param cm_mod +void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, SolutionStates& solutions) { + // Local aliases for solution arrays + auto& An = solutions.current.get_acceleration(); + auto& Yn = solutions.current.get_velocity(); + auto& Dn = solutions.current.get_displacement(); + const auto& Ao = solutions.old.get_acceleration(); + const auto& Yo = solutions.old.get_velocity(); + const auto& Do = solutions.old.get_displacement(); + using namespace consts; #define n_debug_calc_der_cpl_bc @@ -110,22 +118,22 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) 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, com_mod.Yo, 0, nsd-1, false, cfg_o); - cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yn, 0, nsd-1, false, cfg_n); + 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); cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; - #ifdef debug_calc_der_cpl_bc + #ifdef debug_calc_der_cpl_bc dmsg << "iBC_Neu "; dmsg << "cplBC.fa[ptr].Qo: " << cplBC.fa[ptr].Qo; dmsg << "cplBC.fa[ptr].Qn: " << cplBC.fa[ptr].Qn; #endif } - // Compute avg pressures at 3D Dirichlet boundaries at timesteps n and n+1 + // 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, com_mod.Yo, nsd) / area; - cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yn, nsd) / 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].Qo = 0.0; cplBC.fa[ptr].Qn = 0.0; #ifdef debug_calc_der_cpl_bc @@ -525,8 +533,13 @@ void RCR_Integ_X(ComMod& com_mod, const CmMod& cm_mod, int istat) /// /// Replaces 'SUBROUTINE RCRINIT()' // -void rcr_init(ComMod& com_mod, const CmMod& cm_mod) +void rcr_init(ComMod& com_mod, const CmMod& cm_mod, SolutionStates& solutions) { + // Local aliases for old solution arrays + const auto& Ao = solutions.old.get_acceleration(); + const auto& Do = solutions.old.get_displacement(); + const auto& Yo = solutions.old.get_velocity(); + using namespace consts; const int iEq = 0; @@ -541,15 +554,15 @@ void rcr_init(ComMod& com_mod, const CmMod& cm_mod) int ptr = bc.cplBCptr; if (!utils::btest(bc.bType, iBC_RCR)) { - continue; + continue; } if (ptr != -1) { 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, com_mod.Yo, 0, nsd-1); - double Po = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, nsd) / 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; cplBC.xo[ptr] = Po - (Qo * cplBC.fa[ptr].RCR.Rp); } else { cplBC.xo[ptr] = cplBC.fa[ptr].RCR.Xo; @@ -560,7 +573,7 @@ void rcr_init(ComMod& com_mod, const CmMod& cm_mod) /// @brief Below defines the SET_BC methods for the Coupled Momentum Method (CMM) // -void set_bc_cmm(ComMod& com_mod, const CmMod& cm_mod, const Array& Ag, const Array& Dg ) +void set_bc_cmm(ComMod& com_mod, const CmMod& cm_mod, const Array& Ag, const Array& Dg, SolutionStates& solutions) { using namespace consts; @@ -571,7 +584,7 @@ void set_bc_cmm(ComMod& com_mod, const CmMod& cm_mod, const Array& Ag, c auto& bc = eq.bc[iBc]; if (!utils::btest(bc.bType,iBC_CMM)) { - continue; + continue; } int iFa = bc.iFa; @@ -581,11 +594,11 @@ void set_bc_cmm(ComMod& com_mod, const CmMod& cm_mod, const Array& Ag, c throw std::runtime_error("[set_bc_cmm] CMM equation is formulated for tetrahedral elements (volume) and triangular (surface) elements"); } - set_bc_cmm_l(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Ag, Dg); + set_bc_cmm_l(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Ag, Dg, solutions); } } -void set_bc_cmm_l(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& Ag, const Array& Dg ) +void set_bc_cmm_l(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& Ag, const Array& Dg, SolutionStates& solutions) { using namespace consts; @@ -637,7 +650,7 @@ void set_bc_cmm_l(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, con vwp = vwp / 3.0; // Add CMM BCs contributions to the LHS/RHS - cmm::cmm_b(com_mod, lFa, e, al, dl, xl, bfl, pSl, vwp, ptr); + cmm::cmm_b(com_mod, lFa, e, al, dl, xl, bfl, pSl, vwp, ptr, solutions); } } @@ -645,8 +658,16 @@ void set_bc_cmm_l(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, con /// @brief Coupled BC quantities are computed here. /// Reproduces the Fortran 'SETBCCPL()' subrotutine. // -void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) +void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, SolutionStates& solutions) { + // Local aliases for solution arrays + auto& An = solutions.current.get_acceleration(); + auto& Yn = solutions.current.get_velocity(); + auto& Dn = solutions.current.get_displacement(); + const auto& Ao = solutions.old.get_acceleration(); + 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; @@ -654,8 +675,7 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) const int nsd = com_mod.nsd; const int tnNo = com_mod.tnNo; auto& cplBC = com_mod.cplBC; - auto& Yo = com_mod.Yo; - auto& Yn = com_mod.Yn; + // Yo, Ao, Do now passed as parameters const int iEq = 0; auto& eq = com_mod.eq[iEq]; @@ -666,10 +686,10 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) auto cfg_o = MechanicalConfigurationType::reference; auto cfg_n = MechanicalConfigurationType::reference; - // If coupling scheme is implicit, calculate updated pressure and flowrate + // If coupling scheme is implicit, calculate updated pressure and flowrate // from 0D, as well as resistance from 0D using finite difference. - if (cplBC.schm == CplBCType::cplBC_I) { - calc_der_cpl_bc(com_mod, cm_mod); + if (cplBC.schm == CplBCType::cplBC_I) { + calc_der_cpl_bc(com_mod, cm_mod, solutions); // If coupling scheme is semi-implicit or explicit, only calculated updated // pressure and flowrate from 0D @@ -685,8 +705,8 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) 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); - baf_ini_ns::bc_ini(com_mod, cm_mod, eq.bc[iBc], lFa); + double area = all_fun::integ(com_mod, cm_mod, lFa, sA, solutions, false, consts::MechanicalConfigurationType::reference); + baf_ini_ns::bc_ini(com_mod, cm_mod, eq.bc[iBc], lFa, solutions); int ptr = bc.cplBCptr; @@ -720,15 +740,15 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) 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, nsd-1, false, cfg_o); - cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, 0, nsd-1, false, cfg_n); + 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); 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) / area; - cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, nsd) / area; + 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; cplBC.fa[ptr].Qo = 0.0; cplBC.fa[ptr].Qn = 0.0; } @@ -771,8 +791,16 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) /// /// Reproduces 'SUBROUTINE SETBCDIR(lA, lY, lD)' // -void set_bc_dir(ComMod& com_mod, Array& lA, Array& lY, Array& lD) +void set_bc_dir(ComMod& com_mod, SolutionStates& solutions) { + // Local aliases for solution arrays + auto& lA = solutions.current.get_acceleration(); + auto& lY = solutions.current.get_velocity(); + auto& lD = solutions.current.get_displacement(); + const auto& Yo = solutions.old.get_velocity(); + const auto& Ao = solutions.old.get_acceleration(); + const auto& Do = solutions.old.get_displacement(); + using namespace consts; #define n_set_bc_dir @@ -932,8 +960,8 @@ void set_bc_dir(ComMod& com_mod, Array& lA, Array& lY, Array& lA, Array& lY, Array& lA, Array& lY, Array& lA, Array& lY, Array& Yg, const Array& Dg) +void set_bc_dir_w(ComMod& com_mod, const Array& Yg, const Array& Dg, SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; const int cEq = com_mod.cEq; @@ -1058,14 +1089,17 @@ void set_bc_dir_w(ComMod& com_mod, const Array& Yg, const Array& if (!bc.weakDir) { continue; } - set_bc_dir_wl(com_mod, bc, com_mod.msh[iM], com_mod.msh[iM].fa[iFa], Yg, Dg); + set_bc_dir_wl(com_mod, bc, com_mod.msh[iM], com_mod.msh[iM].fa[iFa], Yg, Dg, solutions); } } /// @brief Reproduces Fortran 'SETBCDIRWL'. // -void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const faceType& lFa, const Array& Yg, const Array& Dg) +void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const faceType& lFa, const Array& Yg, const Array& Dg, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; #define n_debug_set_bc_dir_wl @@ -1252,7 +1286,7 @@ void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const 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); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoNb, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); double Jac = sqrt(utils::norm(nV)); nV = nV / Jac; double w = lFa.w(g) * Jac; @@ -1293,8 +1327,12 @@ void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const /// @brief Set outlet BCs. // -void set_bc_neu(ComMod& com_mod, const CmMod& cm_mod, const Array& Yg, const Array& Dg) +void set_bc_neu(ComMod& com_mod, const CmMod& cm_mod, const Array& Yg, const Array& Dg, SolutionStates& solutions) { + // Local aliases for solution arrays + auto& Yn = solutions.current.get_velocity(); + const auto& Do = solutions.old.get_displacement(); + using namespace consts; #define n_debug_set_bc_neu @@ -1325,18 +1363,22 @@ void set_bc_neu(ComMod& com_mod, const CmMod& cm_mod, const Array& Yg, c dmsg << "iFa: " << iFa+1; dmsg << "name: " << com_mod.msh[iM].fa[iFa].name; #endif - set_bc_neu_l(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa], Yg, Dg); + set_bc_neu_l(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa], Yg, Dg, solutions); - } else if (utils::btest(bc.bType,iBC_trac)) { - set_bc_trac_l(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa]); + } else if (utils::btest(bc.bType,iBC_trac)) { + set_bc_trac_l(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa], solutions); } } } /// @brief Set Neumann BC // -void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa, const Array& Yg, const Array& Dg) +void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa, const Array& Yg, const Array& Dg, SolutionStates& solutions) { + // Local aliases for solution arrays + auto& Yn = solutions.current.get_velocity(); + const auto& Do = solutions.old.get_displacement(); + using namespace consts; #define n_debug_set_bc_neu_l @@ -1349,7 +1391,6 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const auto& eq = com_mod.eq[cEq]; int tnNo = com_mod.tnNo; int nsd = com_mod.nsd; - auto& Yn = com_mod.Yn; int nNo = lFa.nNo; Vector h(1), rtmp(1); @@ -1385,7 +1426,7 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const } } else if (utils::btest(lBc.bType,iBC_res)) { - h(0) = lBc.r * all_fun::integ(com_mod, cm_mod, lFa, Yn, eq.s, eq.s+nsd-1); + 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)) { h(0) = lBc.g; @@ -1426,10 +1467,10 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const // Add Neumann BCs contribution to the residual (and tangent if flwP) // if (lBc.flwP) { - eq_assem::b_neu_folw_p(com_mod, lBc, lFa, hg, Dg); + eq_assem::b_neu_folw_p(com_mod, lBc, lFa, hg, Dg, solutions); } else { - eq_assem::b_assem_neu_bc(com_mod, lFa, hg, Yg); + eq_assem::b_assem_neu_bc(com_mod, lFa, hg, Yg, solutions); } @@ -1439,20 +1480,23 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const // a follower pressure load (struct/ustruct) or a moving mesh (FSI) if (utils::btest(lBc.bType, iBC_res)) { if (lBc.flwP || com_mod.mvMsh) { - eq_assem::fsi_ls_upd(com_mod, lBc, lFa); + eq_assem::fsi_ls_upd(com_mod, lBc, lFa, solutions); } } // Now treat Robin BC (stiffness and damping) here if (lBc.robin_bc.is_initialized()) { - set_bc_rbnl(com_mod, lFa, lBc.robin_bc, Yg, Dg); + set_bc_rbnl(com_mod, lFa, lBc.robin_bc, Yg, Dg, solutions); } } /// @brief Set Robin BC contribution to residual and tangent // void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const RobinBoundaryCondition& robin_bc, - const Array& Yg, const Array& Dg) + const Array& Yg, const Array& Dg, SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; #define n_debug_set_bc_rbnl_l @@ -1506,10 +1550,10 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const RobinBoundaryCondit 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); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); double Jac = sqrt(utils::norm(nV)); nV = nV / Jac; - double w = lFa.w(g) * Jac; + double w = lFa.w(g) * Jac; N = lFa.N.col(g); Vector u(nsd), ud(nsd); @@ -1700,12 +1744,15 @@ void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const RobinBoundaryCondit /// @brief Set Traction BC // -void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa) +void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa, SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + using namespace consts; - #define n_debug_set_bc_trac_l - #ifdef debug_set_bc_trac_l + #define n_debug_set_bc_trac_l + #ifdef debug_set_bc_trac_l DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); #endif @@ -1786,7 +1833,7 @@ void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, cons 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); + nn::gnnb(com_mod, lFa, e, g, nsd, nsd-1, eNoN, Nx, nV, solutions, consts::MechanicalConfigurationType::reference); double Jac = sqrt(utils::norm(nV)); double w = lFa.w(g)*Jac; N = lFa.N.col(g); diff --git a/Code/Source/solver/set_bc.h b/Code/Source/solver/set_bc.h index 5f3be8657..9cfda19b5 100644 --- a/Code/Source/solver/set_bc.h +++ b/Code/Source/solver/set_bc.h @@ -12,33 +12,33 @@ namespace set_bc { -void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod); +void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod, SolutionStates& solutions); void cplBC_Integ_X(ComMod& com_mod, const CmMod& cm_mod, const bool RCRflag); void genBC_Integ_X(ComMod& com_mod, const CmMod& cm_mod, const std::string& genFlag); -void rcr_init(ComMod& com_mod, const CmMod& cm_mod); +void rcr_init(ComMod& com_mod, const CmMod& cm_mod, SolutionStates& solutions); void RCR_Integ_X(ComMod& com_mod, const CmMod& cm_mod, int istat); -void set_bc_cmm(ComMod& com_mod, const CmMod& cm_mod, const Array& Ag, const Array& Dg); -void set_bc_cmm_l(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& Ag, const Array& Dg ); +void set_bc_cmm(ComMod& com_mod, const CmMod& cm_mod, const Array& Ag, const Array& Dg, SolutionStates& solutions); +void set_bc_cmm_l(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& Ag, const Array& Dg, SolutionStates& solutions); -void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod); +void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod, SolutionStates& solutions); -void set_bc_dir(ComMod& com_mod, Array& lA, Array& lY, Array& lD); +void set_bc_dir(ComMod& com_mod, SolutionStates& solutions); void set_bc_dir_l(ComMod& com_mod, const bcType& lBc, const faceType& lFa, Array& lA, Array& lY, int lDof); -void set_bc_dir_w(ComMod& com_mod, const Array& Yg, const Array& Dg); -void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const faceType& lFa, const Array& Yg, const Array& Dg); +void set_bc_dir_w(ComMod& com_mod, const Array& Yg, const Array& Dg, SolutionStates& solutions); +void set_bc_dir_wl(ComMod& com_mod, const bcType& lBc, const mshType& lM, const faceType& lFa, const Array& Yg, const Array& Dg, const SolutionStates& solutions); -void set_bc_neu(ComMod& com_mod, const CmMod& cm_mod, const Array& Yg, const Array& Dg); -void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa, const Array& Yg, const Array& Dg); +void set_bc_neu(ComMod& com_mod, const CmMod& cm_mod, const Array& Yg, const Array& Dg, SolutionStates& solutions); +void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa, const Array& Yg, const Array& Dg, SolutionStates& solutions); void set_bc_rbnl(ComMod& com_mod, const faceType& lFa, const RobinBoundaryCondition& robin_bc, - const Array& Yg, const Array& Dg); + const Array& Yg, const Array& Dg, SolutionStates& solutions); -void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa); +void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa, SolutionStates& solutions); void set_bc_undef_neu(ComMod& com_mod); diff --git a/Code/Source/solver/txt.cpp b/Code/Source/solver/txt.cpp index 4b52cca07..40a6e9157 100644 --- a/Code/Source/solver/txt.cpp +++ b/Code/Source/solver/txt.cpp @@ -140,8 +140,14 @@ void create_volume_integral_file(const ComMod& com_mod, CmMod& cm_mod, const eqT // init_write - if true then this is the start of the simulation and // so create a new file to initialize output. // -void txt(Simulation* simulation, const bool init_write) +void txt(Simulation* simulation, const bool init_write, const SolutionStates& solutions) { + // Local aliases for solution arrays + auto& An = solutions.current.get_acceleration(); + auto& Yn = solutions.current.get_velocity(); + auto& Dn = solutions.current.get_displacement(); + const auto& Do = solutions.old.get_displacement(); + using namespace consts; using namespace utils; @@ -269,7 +275,7 @@ void txt(Simulation* simulation, const bool init_write) case OutputNameType::outGrp_A: for (int j = 0; j < tnNo; j++) { for (int i = 0; i < l; i++) { - tmpV(i,j) = com_mod.An(i+s,j); + tmpV(i,j) = An(i+s,j); } } break; @@ -277,7 +283,7 @@ void txt(Simulation* simulation, const bool init_write) case OutputNameType::outGrp_Y: for (int j = 0; j < tnNo; j++) { for (int i = 0; i < l; i++) { - tmpV(i,j) = com_mod.Yn(i+s,j); + tmpV(i,j) = Yn(i+s,j); } } @@ -289,7 +295,7 @@ void txt(Simulation* simulation, const bool init_write) case OutputNameType::outGrp_D: for (int j = 0; j < tnNo; j++) { for (int i = 0; i < l; i++) { - tmpV(i,j) = com_mod.Dn(i+s,j); + tmpV(i,j) = Dn(i+s,j); } } break; @@ -297,7 +303,7 @@ void txt(Simulation* simulation, const bool init_write) case OutputNameType::outGrp_WSS: case OutputNameType::outGrp_vort: case OutputNameType::outGrp_trac: - post::all_post(simulation, tmpV, com_mod.Yn, com_mod.Dn, oGrp, iEq); + post::all_post(simulation, tmpV, Yn, Dn, oGrp, iEq); for (int a = 0; a < tnNo; a++) { auto vec = tmpV.col(a, {0,nsd-1}); tmpV(0,a) = sqrt(norm(vec)); @@ -310,13 +316,13 @@ void txt(Simulation* simulation, const bool init_write) case OutputNameType::outGrp_divV: case OutputNameType::outGrp_J: case OutputNameType::outGrp_mises: - post::all_post(simulation, tmpV, com_mod.Yn, com_mod.Dn, oGrp, iEq); + post::all_post(simulation, tmpV, Yn, Dn, oGrp, iEq); break; case OutputNameType::outGrp_absV: for (int a = 0; a < tnNo; a++) { for (int i = 0; i < l; i++) { - tmpV(i,a) = com_mod.Yn(i,a) - com_mod.Yn(i+nsd+1,a); + tmpV(i,a) = Yn(i,a) - Yn(i+nsd+1,a); } } break; @@ -367,10 +373,10 @@ void txt(Simulation* simulation, const bool init_write) } else { if (output_options.boundary_integral) { - write_boundary_integral_data(com_mod, cm_mod, eq, l, boundary_file_name, tmpV, div, pflag); + write_boundary_integral_data(com_mod, cm_mod, eq, l, boundary_file_name, tmpV, div, pflag, solutions); } if (output_options.volume_integral) { - write_volume_integral_data(com_mod, cm_mod, eq, l, volume_file_name, tmpV, div, pflag); + write_volume_integral_data(com_mod, cm_mod, eq, l, volume_file_name, tmpV, div, pflag, solutions); } } } @@ -407,9 +413,12 @@ void txt(Simulation* simulation, const bool init_write) /// /// NOTE: Be carefu of a potential indexing problem here because 'm' is a length and not an index. // -void write_boundary_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, - const std::string file_name, const Array& tmpV, const bool div, const bool pFlag) +void write_boundary_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, + const std::string file_name, const Array& tmpV, const bool div, const bool pFlag, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + #define n_debug_write_boundary_integral_data #ifdef debug_write_boundary_integral_data DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -451,17 +460,17 @@ 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) / tmp; + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference) / tmp; } else { if (pFlag && lTH) { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, true); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, std::nullopt, true, consts::MechanicalConfigurationType::reference); } else { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, std::nullopt, false, consts::MechanicalConfigurationType::reference); } } } else if (m == nsd) { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, m-1); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, solutions, m-1, false, consts::MechanicalConfigurationType::reference); } else { throw std::runtime_error("WTXT only accepts 1 and nsd"); } @@ -482,9 +491,12 @@ void write_boundary_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eq /// /// NOTE: Be carefu of a potential indexing problem here because 'm' is a length and not an index. // -void write_volume_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, - const std::string file_name, const Array& tmpV, const bool div, const bool pFlag) +void write_volume_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, + const std::string file_name, const Array& tmpV, const bool div, const bool pFlag, const SolutionStates& solutions) { + // Local alias for old displacement + const auto& Do = solutions.old.get_displacement(); + #define n_debug_write_volume_integral_data #ifdef debug_write_volume_integral_data DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -517,9 +529,9 @@ void write_volume_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqTy if (div) { tmp = dmn.v; - tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 0, m-1) / tmp; + tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 0, m-1, solutions, pFlag) / tmp; } else { - tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 0, m-1, pFlag); + tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 0, m-1, solutions, pFlag); } if (com_mod.cm.mas(cm_mod)) { diff --git a/Code/Source/solver/txt.h b/Code/Source/solver/txt.h index c957c0c66..b70a223b9 100644 --- a/Code/Source/solver/txt.h +++ b/Code/Source/solver/txt.h @@ -16,13 +16,13 @@ void create_boundary_integral_file(const ComMod& com_mod, CmMod& cm_mod, const e void create_volume_integral_file(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const std::string& file_name); -void txt(Simulation* simulation, const bool flag); +void txt(Simulation* simulation, const bool flag, const SolutionStates& solutions); -void write_boundary_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, - const std::string file_name, const Array& tmpV, const bool div, const bool pFlag); +void write_boundary_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, + const std::string file_name, const Array& tmpV, const bool div, const bool pFlag, const SolutionStates& solutions); -void write_volume_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, - const std::string file_name, const Array& tmpV, const bool div, const bool pFlag); +void write_volume_integral_data(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, + const std::string file_name, const Array& tmpV, const bool div, const bool pFlag, const SolutionStates& solutions); }; diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp index cadbb34f6..4b073348a 100644 --- a/Code/Source/solver/uris.cpp +++ b/Code/Source/solver/uris.cpp @@ -17,9 +17,12 @@ namespace uris { -/// @brief This subroutine computes the mean pressure and flux on the -/// immersed surface -void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { +/// @brief This subroutine computes the mean pressure and flux on the +/// immersed surface +void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionStates& solutions) { + // Local aliases for solution arrays + auto& Yn = solutions.current.get_velocity(); + auto& Do = solutions.old.get_displacement(); #define n_debug_uris_meanp #ifdef debug_uris_meanp DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -41,7 +44,6 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { // auto& An = com_mod.An; // auto& Ad = com_mod.Ad; // auto& Dn = com_mod.Dn; - auto& Yn = com_mod.Yn; // Let's conpute the mean pressure in the two regions of the fluid mesh // For the moment let's define a flag IdSubDmn(size the number of elements) @@ -68,7 +70,7 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } for (int iM = 0; iM < com_mod.nMsh; iM++) { - volU += all_fun::integ(com_mod, cm_mod, iM, sUPS); + volU += all_fun::integ(com_mod, cm_mod, iM, sUPS, solutions); } @@ -84,7 +86,7 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } for (int iM = 0; iM < com_mod.nMsh; iM++) { - volD += all_fun::integ(com_mod, cm_mod, iM, sDST); + volD += all_fun::integ(com_mod, cm_mod, iM, sDST, solutions); } // Print volume messages. @@ -106,7 +108,7 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { tmpV(0,j) = Yn(s,j)*sUPS(j); } for (int iM = 0; iM < com_mod.nMsh; iM++) { - meanPU += all_fun::integ(com_mod, cm_mod, iM, tmpV); + meanPU += all_fun::integ(com_mod, cm_mod, iM, tmpV, solutions); } meanPU = meanPU / volU; @@ -114,7 +116,7 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { tmpV(0,j) = Yn(s,j)*sDST(j); } for (int iM = 0; iM < com_mod.nMsh; iM++) { - meanPD += all_fun::integ(com_mod, cm_mod,iM, tmpV); + meanPD += all_fun::integ(com_mod, cm_mod,iM, tmpV, solutions); } meanPD = meanPD / volD; @@ -149,9 +151,12 @@ void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } -/// @brief This subroutine computes the mean velocity in the fluid elements -/// near the immersed surface -void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { +/// @brief This subroutine computes the mean velocity in the fluid elements +/// near the immersed surface +void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionStates& solutions) { + // Local aliases for solution arrays + auto& Yn = solutions.current.get_velocity(); + auto& Do = solutions.old.get_displacement(); #define n_debug_uris_meanv #ifdef debug_uris_meanv DebugMsg dmsg(__func__, com_mod.cm.idcm()); @@ -173,7 +178,6 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { // auto& An = com_mod.An; // auto& Ad = com_mod.Ad; // auto& Dn = com_mod.Dn; - auto& Yn = com_mod.Yn; // Let's compute the neighboring region below the valve normal. When // the valve is open, this region should roughly be valve oriface. @@ -193,7 +197,7 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } for (int iM = 0; iM < com_mod.nMsh; iM++) { - volI += all_fun::integ(com_mod, cm_mod, iM, sImm); + volI += all_fun::integ(com_mod, cm_mod, iM, sImm, solutions); } if (cm.mas(cm_mod)) { std::cout << "volume inside " << volI << " for: " << uris_obj.name << std::endl; @@ -219,7 +223,7 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { double meanV = 0.0; for (int iM = 0; iM < com_mod.nMsh; iM++) { - meanV += all_fun::integ(com_mod, cm_mod, iM, tmpVNrm)/volI; + meanV += all_fun::integ(com_mod, cm_mod, iM, tmpVNrm, solutions)/volI; } if (cm.mas(cm_mod)) { @@ -245,10 +249,12 @@ void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { } -/// @brief This subroutine computes the displacement of the immersed +/// @brief This subroutine computes the displacement of the immersed /// surface with fem projection -void uris_update_disp(ComMod& com_mod, CmMod& cm_mod) { - #define n_debug_uris_update_disp +void uris_update_disp(ComMod& com_mod, CmMod& cm_mod, SolutionStates& solutions) { + // Local alias for solution array + const auto& Do = solutions.old.get_displacement(); + #define n_debug_uris_update_disp #ifdef debug_uris_update_disp DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); @@ -313,10 +319,10 @@ void uris_update_disp(ComMod& com_mod, CmMod& cm_mod) { d = 0.0; for (int a = 0; a < mesh.eNoN; a++) { int Ac = mesh.IEN(a, iEln); - //We have to use Do because Dn contains the result coming from the solid - d(0) += N(a)*com_mod.Do(nsd+1, Ac); - d(1) += N(a)*com_mod.Do(nsd+2, Ac); - d(2) += N(a)*com_mod.Do(nsd+3, Ac); + //We have to use Do because Dn contains the result coming from the solid + d(0) += N(a)*Do(nsd+1, Ac); + d(1) += N(a)*Do(nsd+2, Ac); + d(2) += N(a)*Do(nsd+3, Ac); } // update uris disp localYd.set_col(nd, d); diff --git a/Code/Source/solver/uris.h b/Code/Source/solver/uris.h index cb83ba299..90678a854 100644 --- a/Code/Source/solver/uris.h +++ b/Code/Source/solver/uris.h @@ -9,11 +9,11 @@ namespace uris { -void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris); // done +void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionStates& solutions); // done -void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris); // done +void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris, const SolutionStates& solutions); // done -void uris_update_disp(ComMod& com_mod, CmMod& cm_mod); +void uris_update_disp(ComMod& com_mod, CmMod& cm_mod, SolutionStates& solutions); void uris_find_tetra(ComMod& com_mod, CmMod& cm_mod, const int iUris);