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<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)
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<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
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<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)
 {
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<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);
 
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<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();
+  }
 }
 
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<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) {
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 6ec5c6df..f7d66007 100644
--- a/Code/Source/svFSI/main.cpp
+++ b/Code/Source/svFSI/main.cpp
@@ -57,6 +57,8 @@
 #include <stdlib.h>
 #include <iomanip>
 #include <iostream>
+#include <cmath>
+#include <fstream>
 
 /// @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<int>(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);
@@ -706,6 +795,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<double>& Ag, Array<double>& Yg, Array<do
     dmsg << "coef: " << coef[0] << " " << coef[1] << " " << coef[2] << " " << coef[3];
     #endif
 
-    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);
-      }
+    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++) {
+                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);
+            }
+        }
     }
   }
 
diff --git a/Code/Source/svFSI/read_files.cpp b/Code/Source/svFSI/read_files.cpp
index 8996f393..41312c00 100644
--- a/Code/Source/svFSI/read_files.cpp
+++ b/Code/Source/svFSI/read_files.cpp
@@ -1730,11 +1730,12 @@ void read_files(Simulation* simulation, const std::string& file_name)
     }     
 
     if (eq.phys == EquationType::phys_heatF) {   
-      auto& eq1_params = simulation->parameters.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 965248be..27e87910 100644
--- a/Code/Source/svFSI/vtk_xml.cpp
+++ b/Code/Source/svFSI/vtk_xml.cpp
@@ -602,6 +602,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
 //----------------
@@ -1032,11 +1079,29 @@ void write_vtus(Simulation* simulation, const Array<double>& lA, const Array<dou
           break;
 
           case OutputType::outGrp_Y:
-            for (int a = 0; a < msh.nNo; a++) {
-              int Ac = msh.gN(a);
-              for (int i = 0; i < l; i++) {
-                 d[iM].x(i+is,a) = lY(i+s,Ac);
-              }
+            if (eq.phys != EquationType::phys_heatF) {
+                for (int a = 0; a < msh.nNo; a++) {
+                    int Ac = msh.gN(a);
+                    for (int i = 0; i < l; i++) {
+                        d[iM].x(i + is, a) = lY(i + s, Ac);
+                    }
+                }
+            } else {
+                if (eq.output[iOut].name == "Velocity") {
+                    for (int a = 0; a < msh.nNo; a++) {
+                        int Ac = msh.gN(a);
+                        for (int i = 0; i < l; i++) {
+                            d[iM].x(i + is, a) = lY(i, Ac);
+                        }
+                    }
+                } else {
+                    for (int a = 0; a < msh.nNo; a++) {
+                        int Ac = msh.gN(a);
+                        for (int i = 0; i < l; i++) {
+                            d[iM].x(i + is, a) = lY(i + s, Ac);
+                        }
+                    }
+                }
             }
           break;
 
diff --git a/Code/Source/svFSI/vtk_xml.h b/Code/Source/svFSI/vtk_xml.h
index 63370b28..043cce85 100644
--- a/Code/Source/svFSI/vtk_xml.h
+++ b/Code/Source/svFSI/vtk_xml.h
@@ -47,6 +47,8 @@ void read_vtp_pdata(const std::string& fName, const std::string& kwrd, const int
 
 void read_vtu(const std::string& file_name, mshType& mesh);
 
+void read_precomputed_solution_vtu(const std::string& file_name, const std::string& field_name, mshType& mesh);
+
 void read_vtu_pdata(const std::string& fName, const std::string& kwrd, const int nsd, const int m, const int idx, mshType& mesh);
 
 void read_vtus(Simulation* simulation, Array<double>& lA, Array<double>& lY, Array<double>& lD, const std::string& fName);
diff --git a/Code/Source/svFSI/vtk_xml_parser.cpp b/Code/Source/svFSI/vtk_xml_parser.cpp
index bbe97d49..876eca71 100644
--- a/Code/Source/svFSI/vtk_xml_parser.cpp
+++ b/Code/Source/svFSI/vtk_xml_parser.cpp
@@ -35,6 +35,7 @@
 
 #include "vtk_xml_parser.h" 
 #include "Array.h" 
+#include "Array3.h"
 
 #include <vtkDoubleArray.h>
 #include "vtkCellData.h"
@@ -42,6 +43,7 @@
 #include <vtkGeometryFilter.h>
 #include <vtkIntArray.h>
 #include <vtkPointData.h>
+#include <vtkDataArray.h>
 #include <vtkPolyData.h>
 #include <vtkSmartPointer.h>
 #include <vtkUnsignedCharArray.h>
@@ -52,6 +54,8 @@
 #include <string>
 #include <map>
 #include <vector>
+#include <algorithm>
+#include <cctype>
 
 namespace vtk_xml_parser {
 
@@ -828,6 +832,7 @@ void load_vtu(const std::string& file_name, mshType& mesh)
   store_element_conn(vtk_ugrid, mesh);
 }
 
+
 /// @brief Store a surface mesh read from a VTK .vtu file into a Face object.
 //
 void load_vtu(const std::string& file_name, faceType& face)
@@ -874,5 +879,78 @@ void load_vtu(const std::string& file_name, faceType& face)
   store_element_ids(vtk_ugrid, face);
 }
 
+=======
+
+/// @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<vtkXMLUnstructuredGridReader>::New();
+    reader->SetFileName(file_name.c_str());
+    reader->Update();
+    vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid = reader->GetOutput();
+    vtkIdType num_nodes = vtk_ugrid->GetNumberOfPoints();
+    int array_count = 0;
+    std::vector<std::pair<std::string, int>> 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<std::string, int>& a, const std::pair<std::string, int>& 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 5e8587af..fa89d710 100644
--- a/Code/Source/svFSI/vtk_xml_parser.h
+++ b/Code/Source/svFSI/vtk_xml_parser.h
@@ -57,6 +57,7 @@ void load_vtu(const std::string& file_name, mshType& mesh);
 
 void load_vtu(const std::string& file_name, faceType& face);
 
+void load_time_varying_field_vtu(const std::string file_name, const std::string field_name, mshType& mesh);
 
 };
 
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:
+
+```
+   <Use_precomputed_solution> true </Use_precomputed_solution>
+   <Precomputed_solution_file_path> precomputed_velocity.vtu </Precomputed_solution_file_path>
+   <Precomputed_solution_field_name> Velocity </Precomputed_solution_field_name>
+```
+
+This tell the solver that this is a one-way coupling study, and the dye is passively transported by the flow.
+
+```
+   <Output type="Alias" >
+     <Temperature> Concentration </Temperature>
+   </Output>
+```
+
+This renames  "Temperature" to "Concentration" for ease of interpretation.
+
+```
+   <Add_BC name="lumen_inlet" >
+      <Type> Dirichlet </Type>
+      <Time_dependence> Steady </Time_dependence>
+      <Value> 1.0 </Value>
+   </Add_BC>
+```
+
+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 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<svFSIFile version="0.1">
+
+<GeneralSimulationParameters>
+  <Continue_previous_simulation> 0 </Continue_previous_simulation>
+  <Number_of_spatial_dimensions> 3 </Number_of_spatial_dimensions> 
+  <Number_of_time_steps> 1 </Number_of_time_steps> 
+  <Time_step_size> 0.01 </Time_step_size> 
+  <Use_precomputed_solution> true </Use_precomputed_solution>
+  <Precomputed_solution_file_path> precomputed_velocity.vtu </Precomputed_solution_file_path>
+  <Precomputed_solution_field_name> Velocity </Precomputed_solution_field_name>
+  <Spectral_radius_of_infinite_time_step> 0.50 </Spectral_radius_of_infinite_time_step> 
+  <Searched_file_name_to_trigger_stop> STOP_SIM </Searched_file_name_to_trigger_stop> 
+
+  <Save_results_to_VTK_format> 1 </Save_results_to_VTK_format> 
+  <Name_prefix_of_saved_VTK_files> result </Name_prefix_of_saved_VTK_files> 
+  <Increment_in_saving_VTK_files> 1 </Increment_in_saving_VTK_files> 
+  <Start_saving_after_time_step> 0 </Start_saving_after_time_step> 
+
+  <Increment_in_saving_restart_files> 5 </Increment_in_saving_restart_files> 
+  <Convert_BIN_to_VTK_format> 0 </Convert_BIN_to_VTK_format> 
+
+  <Verbose> 1 </Verbose> 
+  <Warning> 0 </Warning> 
+  <Debug> 1 </Debug> 
+</GeneralSimulationParameters>
+
+<Add_mesh name="msh" > 
+
+  <Mesh_file_path> mesh/mesh-complete.mesh.vtu </Mesh_file_path>
+
+  <Add_face name="lumen_inlet">
+      <Face_file_path> mesh/mesh-surfaces/lumen_inlet.vtp </Face_file_path>
+  </Add_face>
+
+  <Add_face name="lumen_outlet">
+      <Face_file_path> mesh/mesh-surfaces/lumen_outlet.vtp </Face_file_path>
+  </Add_face>
+
+  <Add_face name="lumen_wall">
+      <Face_file_path> mesh/mesh-surfaces/lumen_wall.vtp </Face_file_path>
+  </Add_face>
+
+</Add_mesh>
+
+
+<Add_equation type="scalarTransport" >
+   <Coupled> false </Coupled>
+   <Min_iterations> 2 </Min_iterations>
+   <Max_iterations> 5 </Max_iterations>
+   <Tolerance> 1e-6 </Tolerance>
+
+   <Conductivity> 1e-8 </Conductivity>
+   <Source_term> 0.0  </Source_term>
+
+   <Output type="Spatial">
+     <Velocity> true </Velocity>
+     <Temperature> true </Temperature>
+   </Output>
+
+   <Output type="Volume_integral" >
+     <Temperature> true </Temperature>
+   </Output>
+
+   <Output type="Alias" >
+     <Temperature> Concentration </Temperature>
+   </Output>
+
+   <LS type="GMRES" >
+      <Preconditioner> FSILS </Preconditioner> 
+      <Tolerance> 1e-6 </Tolerance>
+      <Max_iterations> 100 </Max_iterations>
+      <Krylov_space_dimension> 50 </Krylov_space_dimension>
+   </LS>
+
+   <Add_BC name="lumen_inlet" >
+      <Type> Dirichlet </Type>
+      <Time_dependence> Steady </Time_dependence>
+      <Value> 1.0 </Value>
+   </Add_BC>
+
+</Add_equation>
+
+
+</svFSIFile>
+
+
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 b86b503c..4dd38b48 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"