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

@@ -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 @@ local(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, Array<double
 }
 
 
+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 @@ void initialize(Simulation* simulation, Vector<double>& 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<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 @@ 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) {

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