From 318e3a6dd1b553556d25e97e67aff195cd6ef6c9 Mon Sep 17 00:00:00 2001 From: Zachary Sexton <47196674+zasexton@users.noreply.github.com> Date: Tue, 2 Apr 2024 11:48:02 -0700 Subject: [PATCH] Precomputed velocities #202 (#208) * including functionality to read in solution fields and store them within the mesh * adding function for reading in precomputed solutions and parameters to mesh loading (#202) * adding distribution functionality for Array3 and testing AD integration * adding linear temporal interpolation for precomputed state-variables and adding a velocity output for scalar advection (#202) * fixing initialize check for seperate fluid equation in AD simulation (#202) * removing print statements for checking * removing print statements and adding error checking for precomputed solution dimension consistency (#202) * adding test case for precomputed velocity solutions to the heatf equation, removing commented lines for testing in heatf.cpp, fixing spelling mistakes in ComMod.h (#202) * adding Doxygen compatible comments for new functions in vtk_xml.cpp and vtk_xml_parser.cpp, encapsulated precomputed linear interpolation time advancement in a function to remove from main.cpp (#202) * fixing minor formating in main.cpp (#202) --- Code/Source/svFSI/ComMod.h | 15 ++- Code/Source/svFSI/Parameters.cpp | 5 +- Code/Source/svFSI/Parameters.h | 6 +- Code/Source/svFSI/Simulation.cpp | 8 ++ Code/Source/svFSI/all_fun.cpp | 61 ++++++++++++ Code/Source/svFSI/all_fun.h | 5 + Code/Source/svFSI/distribute.cpp | 32 ++++++- Code/Source/svFSI/initialize.cpp | 40 +++++++- Code/Source/svFSI/load_msh.cpp | 3 + Code/Source/svFSI/main.cpp | 94 ++++++++++++++++++- Code/Source/svFSI/pic.cpp | 32 +++++-- Code/Source/svFSI/read_files.cpp | 7 +- Code/Source/svFSI/set_equation_props.h | 6 +- Code/Source/svFSI/vtk_xml.cpp | 75 ++++++++++++++- Code/Source/svFSI/vtk_xml.h | 2 + Code/Source/svFSI/vtk_xml_parser.cpp | 75 +++++++++++++++ Code/Source/svFSI/vtk_xml_parser.h | 1 + .../cases/fluid/precomputed_dye_AD/README.md | 40 ++++++++ .../mesh/mesh-complete.exterior.vtp | 3 + .../mesh/mesh-complete.mesh.vtu | 3 + .../mesh/mesh-surfaces/lumen_inlet.vtp | 3 + .../mesh/mesh-surfaces/lumen_outlet.vtp | 3 + .../mesh/mesh-surfaces/lumen_wall.vtp | 3 + .../mesh/walls_combined.vtp | 3 + .../precomputed_velocity.vtu | 3 + .../fluid/precomputed_dye_AD/result_001.vtu | 3 + .../cases/fluid/precomputed_dye_AD/svFSI.xml | 87 +++++++++++++++++ tests/conftest.py | 1 + tests/test_fluid.py | 3 + 29 files changed, 599 insertions(+), 23 deletions(-) create mode 100644 tests/cases/fluid/precomputed_dye_AD/README.md create mode 100644 tests/cases/fluid/precomputed_dye_AD/mesh/mesh-complete.exterior.vtp create mode 100644 tests/cases/fluid/precomputed_dye_AD/mesh/mesh-complete.mesh.vtu create mode 100644 tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_inlet.vtp create mode 100644 tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_outlet.vtp create mode 100644 tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_wall.vtp create mode 100644 tests/cases/fluid/precomputed_dye_AD/mesh/walls_combined.vtp create mode 100644 tests/cases/fluid/precomputed_dye_AD/precomputed_velocity.vtu create mode 100644 tests/cases/fluid/precomputed_dye_AD/result_001.vtu create mode 100644 tests/cases/fluid/precomputed_dye_AD/svFSI.xml diff --git a/Code/Source/svFSI/ComMod.h b/Code/Source/svFSI/ComMod.h index 2123299a..dcd0536d 100644 --- a/Code/Source/svFSI/ComMod.h +++ b/Code/Source/svFSI/ComMod.h @@ -937,6 +937,10 @@ class mshType /// davep double Nxx(:,:,:) Array3 Nxx; + /// @brief Solution field (displacement, velocity, pressure, etc.) for a known, potentially + /// time-varying, quantity of interest across a mesh + Array3 Ys; + /// @brief Mesh Name std::string name; @@ -1372,7 +1376,8 @@ class ComMod { /// @brief Postprocess step - convert bin to vtk bool bin2VTK = false; - + /// @brief Whether to use precomputed state-variable solutions + bool usePrecomp = false; //----- int members -----// /// @brief Current domain @@ -1448,6 +1453,9 @@ class ComMod { /// @brief Time step size double dt = 0.0; + /// @brief Time step size of the precomputed state-variables + double precompDt = 0.0; + /// @brief Time double time = 0.0; @@ -1466,6 +1474,11 @@ class ComMod { /// @brief Stop_trigger file name std::string stopTrigName; + /// @brief Precomputed state-variable file name + std::string precompFileName; + + /// @brief Precomputed state-variable field name + std::string precompFieldName; // ALLOCATABLE DATA /// @brief Column pointer (for sparse LHS matrix structure) diff --git a/Code/Source/svFSI/Parameters.cpp b/Code/Source/svFSI/Parameters.cpp index cac21c9b..13ed2884 100644 --- a/Code/Source/svFSI/Parameters.cpp +++ b/Code/Source/svFSI/Parameters.cpp @@ -1803,9 +1803,12 @@ GeneralSimulationParameters::GeneralSimulationParameters() set_parameter("Starting time step", 0, !required, starting_time_step); set_parameter("Time_step_size", 0.0, required, time_step_size); - + set_parameter("Precomputed_time_step_size", 0.0, !required, precomputed_time_step_size); set_parameter("Verbose", false, !required, verbose); set_parameter("Warning", false, !required, warning); + set_parameter("Use_precomputed_solution", false, !required, use_precomputed_solution); + set_parameter("Precomputed_solution_file_path", "", !required, precomputed_solution_file_path); + set_parameter("Precomputed_solution_field_name", "", !required, precomputed_solution_field_name); } void GeneralSimulationParameters::print_parameters() diff --git a/Code/Source/svFSI/Parameters.h b/Code/Source/svFSI/Parameters.h index efb7bcc7..697f7435 100644 --- a/Code/Source/svFSI/Parameters.h +++ b/Code/Source/svFSI/Parameters.h @@ -1207,9 +1207,11 @@ class GeneralSimulationParameters : public ParameterLists Parameter start_averaging_from_zero; Parameter verbose; Parameter warning; + Parameter use_precomputed_solution; Parameter spectral_radius_of_infinite_time_step; Parameter time_step_size; + Parameter precomputed_time_step_size; Parameter increment_in_saving_restart_files; Parameter increment_in_saving_vtk_files; @@ -1223,7 +1225,9 @@ class GeneralSimulationParameters : public ParameterLists Parameter restart_file_name; Parameter searched_file_name_to_trigger_stop; Parameter save_results_in_folder; - Parameter simulation_initialization_file_path; + Parameter simulation_initialization_file_path; + Parameter precomputed_solution_file_path; + Parameter precomputed_solution_field_name; }; /// @brief The FaceParameters class is used to store parameters for the diff --git a/Code/Source/svFSI/Simulation.cpp b/Code/Source/svFSI/Simulation.cpp index 1e375d61..39d8c1e1 100644 --- a/Code/Source/svFSI/Simulation.cpp +++ b/Code/Source/svFSI/Simulation.cpp @@ -100,6 +100,14 @@ void Simulation::set_module_parameters() com_mod.stFileIncr = general.increment_in_saving_restart_files.value(); com_mod.rmsh.isReqd = general.simulation_requires_remeshing.value(); + com_mod.usePrecomp = general.use_precomputed_solution.value(); + com_mod.precompFileName = general.precomputed_solution_file_path.value(); + com_mod.precompFieldName = general.precomputed_solution_field_name.value(); + com_mod.precompDt = general.precomputed_time_step_size.value(); + if ((com_mod.precompDt == 0.0) && (com_mod.usePrecomp)) { + std::cout << "Precomputed time step size is zero. Setting to simulation time step size." << std::endl; + com_mod.precompDt = com_mod.dt; + } // Set simulation parameters. nTs = general.number_of_time_steps.value(); fTmp = general.simulation_initialization_file_path.value(); diff --git a/Code/Source/svFSI/all_fun.cpp b/Code/Source/svFSI/all_fun.cpp index 1285894c..89362b3f 100644 --- a/Code/Source/svFSI/all_fun.cpp +++ b/Code/Source/svFSI/all_fun.cpp @@ -1095,6 +1095,67 @@ local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Array +local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Array3& U){ + if (com_mod.ltg.size() == 0) { + throw std::runtime_error("ltg is not set yet"); + } + + Array3 local_array; + int m; + int n; + int r; + + if (cm.mas(cm_mod)) { + m = U.nrows(); // nsd + r = U.ncols(); // tnNo + n = U.nslices(); // time + if (U.ncols() != com_mod.gtnNo) { + throw std::runtime_error("local_rv is only specified for vector with size gtnNo"); + } + } + + if (cm.seq()) { + local_array.resize(m, com_mod.gtnNo, U.nslices()); + local_array = U; + return local_array; + } + + cm.bcast(cm_mod, &m); + cm.bcast(cm_mod, &n); + cm.bcast(cm_mod, &r); + + local_array.resize(m, com_mod.tnNo, n); + Vector tmpU(m * com_mod.gtnNo * n); + + if (cm.mas(cm_mod)) { + for (int a = 0; a < com_mod.gtnNo; a++) { + int s = m * a; + for (int i = 0; i < n; i++) { + int e = i * m * (com_mod.gtnNo); + for (int j = 0; j < m; j++) { + tmpU(j+s+e) = U(j, a, i); + } + } + } + } + + cm.bcast(cm_mod, tmpU); + + for (int a = 0; a < com_mod.tnNo; a++) { + int Ac = com_mod.ltg[a]; + int s = m * Ac; + for (int i = 0; i < n; i++) { + int e = i * m * (com_mod.gtnNo); + for (int j = 0; j < m; j++) { + local_array(j, a, i) = tmpU(j+s+e); + } + } + } + + return local_array; +} + Vector mkc(const ComMod& com_mod, Vector& U) { diff --git a/Code/Source/svFSI/all_fun.h b/Code/Source/svFSI/all_fun.h index 89377920..681646d8 100644 --- a/Code/Source/svFSI/all_fun.h +++ b/Code/Source/svFSI/all_fun.h @@ -31,7 +31,9 @@ #ifndef ALL_FUN_H #define ALL_FUN_H +#include "Array3.h" #include "Array.h" +#include "Vector.h" #include "ComMod.h" #include "consts.h" @@ -70,8 +72,11 @@ namespace all_fun { double jacobian(ComMod& com_mod, const int nDim, const int eNoN, const Array& x, const Array&Nxi); Vector local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Vector& u); + Array local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Array& u); + Array3 local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Array3& u); + Vector mkc(const ComMod& com_mod, Vector& U); Array mkc(const ComMod& com_mod, Array& U); diff --git a/Code/Source/svFSI/distribute.cpp b/Code/Source/svFSI/distribute.cpp index 5420d99d..c1f9bdfd 100644 --- a/Code/Source/svFSI/distribute.cpp +++ b/Code/Source/svFSI/distribute.cpp @@ -309,6 +309,7 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &com_mod.startTS); cm.bcast(cm_mod, &com_mod.nEq); cm.bcast(cm_mod, &com_mod.dt); + cm.bcast(cm_mod, &com_mod.precompDt); cm.bcast(cm_mod, &com_mod.zeroAve); cm.bcast(cm_mod, &com_mod.cmmInit); @@ -320,6 +321,7 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &simulation->cep_mod.cepEq); + cm.bcast(cm_mod, &com_mod.usePrecomp); if (com_mod.rmsh.isReqd) { auto& rmsh = com_mod.rmsh; cm.bcast_enum(cm_mod, &rmsh.method); @@ -1423,6 +1425,13 @@ void part_face(Simulation* simulation, mshType& lM, faceType& lFa, faceType& gFa /// @brief Reproduces the Fortran 'PARTMSH' subroutine. +/// Parameters for the part_msh function: +/// @param[in] simulation A pointer to the simulation object. +/// @param[in] iM The mesh index. +/// @param[in] lM The local mesh data. +/// @param[in] gmtl The global to local map. +/// @param[in] nP The number of processors. +/// @param[in] wgt The weights. // void part_msh(Simulation* simulation, int iM, mshType& lM, Vector& gmtl, int nP, Vector& wgt) { @@ -2017,6 +2026,27 @@ void part_msh(Simulation* simulation, int iM, mshType& lM, Vector& gmtl, in // Now scattering the sorted lM%INN to all processors MPI_SCATTERV(tempIEN, sCount, disp, mpint, lM%INN, nEl*insd, mpint, master, cm%com(), ierr) */ - } + } + // If necessary, distribute precomputed state-variable data. + // + flag = (lM.Ys.size() != 0); + cm.bcast(cm_mod, &flag); + if (flag){ + #ifdef dbg_part_msh + dmsg << "Distributing precomputed state-variable data " << " ..."; + #endif + Array3 tmpYs; + int nsYs = lM.Ys.nslices(); + if (cm.mas(cm_mod)) { + tmpYs.resize(lM.Ys.nrows(), lM.Ys.ncols(), nsYs); + tmpYs = lM.Ys; + lM.Ys.clear(); + } else { + tmpYs.clear(); + } + lM.Ys.resize(com_mod.nsd, com_mod.tnNo, nsYs); + lM.Ys = all_fun::local(com_mod, cm_mod, cm, tmpYs); + tmpYs.clear(); + } } diff --git a/Code/Source/svFSI/initialize.cpp b/Code/Source/svFSI/initialize.cpp index 33003c58..a7313ddb 100644 --- a/Code/Source/svFSI/initialize.cpp +++ b/Code/Source/svFSI/initialize.cpp @@ -434,16 +434,39 @@ void initialize(Simulation* simulation, Vector& timeP) eq.dof = nsd; } + // This code checks to see if the heatF equation is accompanied + // by a fluid equation. If not, it adds the appropriate number of + // degrees of freedom based on the precomputed state-variable (velocity) + // data. + if (eq.phys == Equation_heatF) { + bool fflag = false; + for (int jEq = 0; jEq < com_mod.nEq; jEq++) { + if (std::set < EquationType > + {Equation_fluid, Equation_FSI, Equation_CMM, Equation_stokes}.count(com_mod.eq[jEq].phys)) { + fflag = true; + } + } + if (com_mod.usePrecomp) { + if (!fflag) { + tDof = tDof + nsd; + } + } else { + if (!fflag) { + throw std::runtime_error( + "HeatF equation must be accompanied by a fluid equation or precomputed velocity data."); + } + } + } eq.pNorm = std::numeric_limits::max(); eq.af = 1.0 / (1.0 + eq.roInf); eq.beta = 0.25 * pow((1.0 + eq.am - eq.af), 2.0); eq.gam = 0.5 + eq.am - eq.af; // These are indexes into arrays so need to be zero-based. + eq.s = tDof; eq.e = tDof + eq.dof - 1; tDof = eq.e + 1; - if (eq.useTLS) { flag = true; } @@ -805,6 +828,21 @@ void zero_init(Simulation* simulation) dmsg.banner(); #endif + // Initialize precomputed state variables + // + + if (com_mod.usePrecomp) { + for (int l = 0; l < com_mod.nMsh; l++) { + auto& msh = com_mod.msh[l]; + 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); + } + } + } + } + // Load any explicitly provided solution variables // if (com_mod.Vinit.size() != 0) { diff --git a/Code/Source/svFSI/load_msh.cpp b/Code/Source/svFSI/load_msh.cpp index 6442dc4f..f70b6b0f 100644 --- a/Code/Source/svFSI/load_msh.cpp +++ b/Code/Source/svFSI/load_msh.cpp @@ -196,6 +196,9 @@ void read_sv(Simulation* simulation, mshType& mesh, const MeshParameters* mesh_p // Note: This may change element node ordering. // auto &com_mod = simulation->get_com_mod(); + if (com_mod.usePrecomp) { + vtk_xml::read_precomputed_solution_vtu(com_mod.precompFileName, com_mod.precompFieldName, mesh); + } if (com_mod.ichckIEN) { read_msh_ns::check_ien(simulation, mesh); } diff --git a/Code/Source/svFSI/main.cpp b/Code/Source/svFSI/main.cpp index 3a1fca6a..0ef712bb 100644 --- a/Code/Source/svFSI/main.cpp +++ b/Code/Source/svFSI/main.cpp @@ -57,6 +57,8 @@ #include #include #include +#include +#include /// @brief Read in a solver XML file and all mesh and BC data. // @@ -83,9 +85,95 @@ 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) { + using namespace consts; + + auto& com_mod = simulation->com_mod; + auto& cm_mod = simulation->cm_mod; + auto& cm = com_mod.cm; + auto& cep_mod = simulation->get_cep_mod(); + + int nTS = com_mod.nTS; + int stopTS = nTS; + int tDof = com_mod.tDof; + int tnNo = com_mod.tnNo; + int nFacesLS = com_mod.nFacesLS; + int nsd = com_mod.nsd; + + auto& Ad = com_mod.Ad; // Time derivative of displacement + 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; + + if (com_mod.usePrecomp) { +#ifdef debug_iterate_solution + dmsg << "Use precomputed values ..." << std::endl; +#endif + // This loop is used to interpolate between known time values of the precomputed + // state-variable solution + for (int l = 0; l < com_mod.nMsh; l++) { + auto lM = com_mod.msh[l]; + if (lM.Ys.nslices() > 1) { + // If there is only one temporal slice, then the solution is assumed constant + // in time and no interpolation is performed + // If there are multiple temporal slices, then the solution is linearly interpolated + // between the known time values and the current time. + double precompDt = com_mod.precompDt; + double preTT = precompDt * (lM.Ys.nslices() - 1); + double cT = cTS * dt; + double rT = std::fmod(cT, preTT); + int n1, n2; + double alpha; + if (precompDt == dt) { + alpha = 0.0; + if (cTS < lM.Ys.nslices()) { + n1 = cTS - 1; + } else { + n1 = cTS % lM.Ys.nslices() - 1; + } + } else { + n1 = static_cast(rT / precompDt) - 1; + alpha = std::fmod(rT, precompDt); + } + n2 = n1 + 1; + for (int i = 0; i < tnNo; i++) { + for (int j = 0; j < nsd; j++) { + if (alpha == 0.0) { + Yn(j, i) = lM.Ys(j, i, n2); + } else { + Yn(j, i) = (1.0 - alpha) * lM.Ys(j, i, n1) + alpha * lM.Ys(j, i, n2); + } + } + } + } else { + for (int i = 0; i < tnNo; i++) { + for (int j = 0; j < nsd; j++) { + Yn(j, i) = lM.Ys(j, i, 0); + } + } + } + } + } +} + /// @brief Iterate the simulation in time. /// -/// Reproduces the outer and inner loops in Fortan MAIN.f. +/// Reproduces the outer and inner loops in Fortan MAIN.f. // void iterate_solution(Simulation* simulation) { @@ -243,6 +331,8 @@ void iterate_solution(Simulation* simulation) set_bc::set_bc_dir(com_mod, An, Yn, Dn); + iterate_precomputed_time(simulation); + // Inner loop for Newton iteration // int inner_count = 1; @@ -282,7 +372,6 @@ void iterate_solution(Simulation* simulation) #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); @@ -707,6 +796,7 @@ int main(int argc, char *argv[]) #endif read_files(simulation, file_name); + // Distribute data to processors. #ifdef debug_main dmsg << "Distribute data to processors " << " ... "; diff --git a/Code/Source/svFSI/pic.cpp b/Code/Source/svFSI/pic.cpp index de6f6356..4bbcebd7 100644 --- a/Code/Source/svFSI/pic.cpp +++ b/Code/Source/svFSI/pic.cpp @@ -526,12 +526,32 @@ void pici(Simulation* simulation, Array& Ag, Array& Yg, Arrayparameters.equation_parameters[0]; + auto& eq1_params = simulation->parameters.equation_parameters[0]; + auto& general_params = simulation->parameters.general_simulation_parameters; auto eq1_type = eq1_params->type.value(); - if ((eq1_type != "fluid") && (eq1_type != "FSI")) { + if ((eq1_type != "fluid") && (eq1_type != "FSI") && (!general_params.use_precomputed_solution.value())) { throw std::runtime_error("heatF equation has to be specified after fluid/FSI equation"); - } + } } if (eq.phys == EquationType::phys_mesh) { diff --git a/Code/Source/svFSI/set_equation_props.h b/Code/Source/svFSI/set_equation_props.h index 56d8d3b6..c9bc5279 100644 --- a/Code/Source/svFSI/set_equation_props.h +++ b/Code/Source/svFSI/set_equation_props.h @@ -264,8 +264,10 @@ SetEquationPropertiesMapType set_equation_props = { read_domain(simulation, eq_params, lEq, propL); - nDOP = {2,1,1,0}; - outPuts = {OutputType::out_temperature, OutputType::out_heatFlux}; + nDOP = {3,1,1,0}; + outPuts = {OutputType::out_temperature, + OutputType::out_heatFlux, + OutputType::out_velocity}; // Set solver parameters. read_ls(simulation, eq_params, SolverType::lSolver_GMRES, lEq); diff --git a/Code/Source/svFSI/vtk_xml.cpp b/Code/Source/svFSI/vtk_xml.cpp index a9ac6999..2d5804e0 100644 --- a/Code/Source/svFSI/vtk_xml.cpp +++ b/Code/Source/svFSI/vtk_xml.cpp @@ -593,6 +593,53 @@ void read_vtu(const std::string& file_name, mshType& mesh) #endif } + +//---------- +// read_precomputed_solution_vtu +//---------- + +/// @brief Read a mesh file with state-variable solution fields from a .vtu or .vtp file. +/// +/// Mesh variables set +/// mesh.Ys = precomputed state-variable solutions (e.g. velocity, pressure, etc.) + +void read_precomputed_solution_vtu(const std::string& file_name, const std::string& field_name, mshType& mesh) +{ + using namespace vtk_xml_parser; + + if (FILE *file = fopen(file_name.c_str(), "r")) { + fclose(file); + } else { + throw std::runtime_error("The VTU mesh file '" + file_name + "' can't be read."); + } + + // Read data from a VTK file. + // + #define n_read_vtu_use_VtkData + #ifdef read_vtu_use_VtkData + auto vtk_data = VtkData::create_reader(file_name); + int num_elems = vtk_data->num_elems(); + int np_elem = vtk_data->np_elem(); + + // Set mesh data. + mesh.nEl = num_elems; + mesh.eNoN = np_elem; + mesh.IEN = vtk_data->get_connectivity(); + mesh.x = vtk_data->get_points(); + + delete vtk_data; + + #else + + auto file_ext = file_name.substr(file_name.find_last_of(".") + 1); + if (file_ext == "vtp") { + vtk_xml_parser::load_time_varying_field_vtu(file_name, field_name, mesh); + } else if (file_ext == "vtu") { + vtk_xml_parser::load_time_varying_field_vtu(file_name, field_name, mesh); + } + #endif +} + //---------------- // read_vtu_pdata //---------------- @@ -1023,11 +1070,29 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array& lA, Array& lY, Array& lD, const std::string& fName); diff --git a/Code/Source/svFSI/vtk_xml_parser.cpp b/Code/Source/svFSI/vtk_xml_parser.cpp index 4d452903..3abd5eb6 100644 --- a/Code/Source/svFSI/vtk_xml_parser.cpp +++ b/Code/Source/svFSI/vtk_xml_parser.cpp @@ -35,12 +35,14 @@ #include "vtk_xml_parser.h" #include "Array.h" +#include "Array3.h" #include #include "vtkCellData.h" #include #include #include +#include #include #include #include @@ -51,6 +53,8 @@ #include #include #include +#include +#include namespace vtk_xml_parser { @@ -750,5 +754,76 @@ void load_vtu(const std::string& file_name, mshType& mesh) store_element_conn(vtk_ugrid, mesh); } + +/// @brief Read a time series field from a VTK .vtu file. +/// +/// Mesh variables set +/// mesh.Ys - time series field data (num_components, num_nodes, num_time_steps) +/// +// +void load_time_varying_field_vtu(const std::string file_name, const std::string field_name, mshType& mesh) +{ + #define n_debug_load_vtu + #ifdef debug_load_vtu + std::cout << "[load_vtu] " << std::endl; + std::cout << "[load_vtu] ===== vtk_xml_parser::load_time_varying_field_vtu ===== " << std::endl; + std::cout << "[load_vtu] file_name: " << file_name << std::endl; + #endif + + auto reader = vtkSmartPointer::New(); + reader->SetFileName(file_name.c_str()); + reader->Update(); + vtkSmartPointer vtk_ugrid = reader->GetOutput(); + vtkIdType num_nodes = vtk_ugrid->GetNumberOfPoints(); + int array_count = 0; + std::vector> array_names; + + if (num_nodes == 0) { + throw std::runtime_error("Failed reading the VTK file '" + file_name + "'."); + } + // Store all array names + for (int i = 0; i < vtk_ugrid->GetPointData()->GetNumberOfArrays(); i++) { + std::string array_name = vtk_ugrid->GetPointData()->GetArrayName(i); + size_t pos = array_name.find(field_name.c_str()); + if (pos != std::string::npos) { + auto not_digit = [](char c) { return !std::isdigit(c); }; + auto it = std::find_if(array_name.rbegin(), array_name.rend(), not_digit); + std::string time_step = std::string(it.base(), array_name.end()); + array_count++; + if (!time_step.empty()) { + array_names.push_back({array_name, std::stoi(time_step)}); + } else { + array_names.push_back({array_name, 0}); + } + } + } + // Check if there are any fields present in the VTK file + if (array_count == 0) { + throw std::runtime_error("No '" + field_name + "' data found in the VTK file '" + file_name + "'."); + } + + // Order all array names by time step + std::sort(array_names.begin(), array_names.end(), [](const std::pair& a, const std::pair& b) { + return a.second < b.second; + }); + // Get the expected number of state-variable components + int num_components = vtk_ugrid->GetPointData()->GetArray(array_names[0].first.c_str())->GetNumberOfComponents(); + mesh.Ys.resize(num_components, num_nodes, array_count); + + for (int i = 0; i < array_count; i++) { + auto array = vtk_ugrid->GetPointData()->GetArray(array_names[i].first.c_str()); + if (array == nullptr) { + throw std::runtime_error("No '" + array_names[i].first + "' data found in the VTK file '" + file_name + "'."); + } + if (array->GetNumberOfComponents() != num_components) { + throw std::runtime_error("The number of components in the field '" + array_names[i].first + "' is not equal to the number of components in the first field."); + } + for (int j = 0; j < num_nodes; j++) { + for (int k = 0; k < num_components; k++) { + mesh.Ys(k, j, i) = array->GetComponent(j, k); + } + } + } +} } // namespace vtk_utils diff --git a/Code/Source/svFSI/vtk_xml_parser.h b/Code/Source/svFSI/vtk_xml_parser.h index 54de5fc9..dec6bf06 100644 --- a/Code/Source/svFSI/vtk_xml_parser.h +++ b/Code/Source/svFSI/vtk_xml_parser.h @@ -55,6 +55,7 @@ void load_vtp(const std::string& file_name, mshType& mesh); void load_vtu(const std::string& file_name, mshType& mesh); +void load_time_varying_field_vtu(const std::string file_name, const std::string field_name, mshType& mesh); }; #endif diff --git a/tests/cases/fluid/precomputed_dye_AD/README.md b/tests/cases/fluid/precomputed_dye_AD/README.md new file mode 100644 index 00000000..03b079d1 --- /dev/null +++ b/tests/cases/fluid/precomputed_dye_AD/README.md @@ -0,0 +1,40 @@ + +# **Problem Description** + +Solve dye transportation with fluid flow in a cylindrical tube with zero neumann boundary condition at the outlet and steady flow at the inlet. The dye is passively transported by the flow through advection and diffusion. + +The input file `svFSI.inp` follows the master input file [`svFSI_master.inp`](./svFSI_master.inp) as a template. Some specific input options are discussed below: + +## Scalar transport equation +Unlike the test case in `04-fluid/02-dye_AD`, the advection-diffusion equation that governs the dye transportation is added to the input file without a proceeding fluid equation. +Here we assume that we already have a velocity field that has been precomputed by some prior physics simulation (denoted `precomputed_velocity.vtu`). The specified velocity field solution +is then taken to advect a scalar field. +This equation is built on top of the heat transfer equation, so some of the terminology in the input file follows those in the heat equation, e.g. "Conductivity", "Temperature". + +Some noticeable setting in the input file are: + +``` + true + precomputed_velocity.vtu + Velocity +``` + +This tell the solver that this is a one-way coupling study, and the dye is passively transported by the flow. + +``` + + Concentration + +``` + +This renames "Temperature" to "Concentration" for ease of interpretation. + +``` + + Dirichlet + Steady + 1.0 + +``` + +The dye is constantly released at the inlet. diff --git a/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-complete.exterior.vtp b/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-complete.exterior.vtp new file mode 100644 index 00000000..2b8a88e9 --- /dev/null +++ b/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-complete.exterior.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e3e23979d7ab682a8778f5d6a61e1d2623597d2e18a7185aa03730249d56c8f9 +size 239624 diff --git a/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-complete.mesh.vtu b/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..bb61e092 --- /dev/null +++ b/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6a31e8fa64cf3d6f14570182dfadd3a19e7f9b1b0f2e0c6280e8e30550ec724f +size 2084502 diff --git a/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_inlet.vtp b/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_inlet.vtp new file mode 100644 index 00000000..828d4a1f --- /dev/null +++ b/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_inlet.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:87e454caaf6c1a588f2027e540c9bf08673984b5304bdfe915105005c92d6b5f +size 15157 diff --git a/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_outlet.vtp b/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_outlet.vtp new file mode 100644 index 00000000..ee73d8d2 --- /dev/null +++ b/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_outlet.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:eda9ee852deed1bb7c79f4611e7f7d21f5ec0d43761f368765e0e0bf9c14d5e2 +size 14566 diff --git a/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_wall.vtp b/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_wall.vtp new file mode 100644 index 00000000..4e6e856d --- /dev/null +++ b/tests/cases/fluid/precomputed_dye_AD/mesh/mesh-surfaces/lumen_wall.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:42f000320f28af94fd4e17ce1907d9dd20b8739062a08ea55fbf3c1cd95ead8d +size 215849 diff --git a/tests/cases/fluid/precomputed_dye_AD/mesh/walls_combined.vtp b/tests/cases/fluid/precomputed_dye_AD/mesh/walls_combined.vtp new file mode 100644 index 00000000..4e6e856d --- /dev/null +++ b/tests/cases/fluid/precomputed_dye_AD/mesh/walls_combined.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:42f000320f28af94fd4e17ce1907d9dd20b8739062a08ea55fbf3c1cd95ead8d +size 215849 diff --git a/tests/cases/fluid/precomputed_dye_AD/precomputed_velocity.vtu b/tests/cases/fluid/precomputed_dye_AD/precomputed_velocity.vtu new file mode 100644 index 00000000..e9387cf2 --- /dev/null +++ b/tests/cases/fluid/precomputed_dye_AD/precomputed_velocity.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e4414f48b24a78632e77c1cff32fb3aca9d5b5d6a53465853ae2db8dde42f8e0 +size 4361565 diff --git a/tests/cases/fluid/precomputed_dye_AD/result_001.vtu b/tests/cases/fluid/precomputed_dye_AD/result_001.vtu new file mode 100644 index 00000000..4cfab084 --- /dev/null +++ b/tests/cases/fluid/precomputed_dye_AD/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ce606ed8358e18433cb2c02c70e5093005116aeeb6d9ea73ce90851607a60b36 +size 4324456 diff --git a/tests/cases/fluid/precomputed_dye_AD/svFSI.xml b/tests/cases/fluid/precomputed_dye_AD/svFSI.xml new file mode 100644 index 00000000..fbcc46c7 --- /dev/null +++ b/tests/cases/fluid/precomputed_dye_AD/svFSI.xml @@ -0,0 +1,87 @@ + + + + + 0 + 3 + 1 + 0.01 + true + precomputed_velocity.vtu + Velocity + 0.50 + STOP_SIM + + 1 + result + 1 + 0 + + 5 + 0 + + 1 + 0 + 1 + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/lumen_inlet.vtp + + + + mesh/mesh-surfaces/lumen_outlet.vtp + + + + mesh/mesh-surfaces/lumen_wall.vtp + + + + + + + false + 2 + 5 + 1e-6 + + 1e-8 + 0.0 + + + true + true + + + + true + + + + Concentration + + + + FSILS + 1e-6 + 100 + 50 + + + + Dirichlet + Steady + 1.0 + + + + + + + + diff --git a/tests/conftest.py b/tests/conftest.py index 8b6b3f81..4c4fd670 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -13,6 +13,7 @@ RTOL = { "Action_potential": 1.0e-10, "Cauchy_stress": 1.0e-4, + "Concentration": 1.0e-10, "Def_grad": 1.0e-10, "Divergence": 1.0e-9, "Displacement": 1.0e-10, diff --git a/tests/test_fluid.py b/tests/test_fluid.py index 641bb0d7..7317f571 100644 --- a/tests/test_fluid.py +++ b/tests/test_fluid.py @@ -23,6 +23,9 @@ def test_dye_AD(n_proc): test_folder = "dye_AD" run_with_reference(base_folder, test_folder, fields, n_proc) +def test_precomputed_dye_AD(n_proc): + test_folder = "precomputed_dye_AD" + run_with_reference(base_folder, test_folder, ['Velocity', 'Concentration'], n_proc) def test_newtonian(n_proc): test_folder = "newtonian"