Precomputed velocities #202

Code/Source/svFSI/ComMod.h

@@ -937,6 +937,10 @@ class mshType
  /// davep double Nxx(:,:,:)
  Array3<double> Nxx;
 
+ /// @brief Solution field (displacement, velocity, pressure, etc.) for a known, potentially
+ /// time-varying, quantity of interest across a mesh
+ Array3<double> 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)

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()

Code/Source/svFSI/Parameters.h

@@ -1207,9 +1207,11 @@ class GeneralSimulationParameters : public ParameterLists
  Parameter<bool> start_averaging_from_zero;
  Parameter<bool> verbose;
  Parameter<bool> warning;
+ Parameter<bool> use_precomputed_solution;
 
  Parameter<double> spectral_radius_of_infinite_time_step;
  Parameter<double> time_step_size;
+ Parameter<double> precomputed_time_step_size;
 
  Parameter<int> increment_in_saving_restart_files;
  Parameter<int> increment_in_saving_vtk_files;
@@ -1223,7 +1225,9 @@ class GeneralSimulationParameters : public ParameterLists
  Parameter<std::string> restart_file_name; 
  Parameter<std::string> searched_file_name_to_trigger_stop; 
  Parameter<std::string> save_results_in_folder; 
- Parameter<std::string> simulation_initialization_file_path; 
+ Parameter<std::string> simulation_initialization_file_path;
+ Parameter<std::string> precomputed_solution_file_path;
+ Parameter<std::string> precomputed_solution_field_name;
 };
 
 /// @brief The FaceParameters class is used to store parameters for the

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();

Code/Source/svFSI/all_fun.cpp

@@ -1095,6 +1095,67 @@
 }
 
 
+Array3<double>
+local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Array3<double>& U){
+ if (com_mod.ltg.size() == 0) {
+ throw std::runtime_error("ltg is not set yet");
+ }
+
+ Array3<double> 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<double> 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<double> 
 mkc(const ComMod& com_mod, Vector<double>& U)
 {

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<double>& x, const Array<double>&Nxi);
 
  Vector<int> local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Vector<int>& u);
+
  Array<double> local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Array<double>& u);
 
+ Array3<double> local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Array3<double>& u);
+
  Vector<double> mkc(const ComMod& com_mod, Vector<double>& U);
  Array<double> mkc(const ComMod& com_mod, Array<double>& U);
 

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<int>& gmtl, int nP, Vector<float>& wgt)
 {
@@ -2017,6 +2026,27 @@ void part_msh(Simulation* simulation, int iM, mshType& lM, Vector<int>& 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<double> 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();
+ }
 }
 

Code/Source/svFSI/initialize.cpp

@@ -434,16 +434,39 @@
  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<double>::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 @@
  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) {

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);
  }