adding distribution functionality for Array3 and testing AD integration

Code/Source/svFSI/all_fun.cpp

@@ -1018,6 +1018,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();
+ r = U.ncols();
+ n = U.nslices();
+ 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, r);
+ Vector<double> tmpU(m * com_mod.gtnNo * r);
+
+ if (cm.mas(cm_mod)) {
+ for (int a = 0; a < com_mod.gtnNo; a++) {
+ int s = m * a;
+ for (int i = 0; i < U.nslices(); i++) {
+ int e = i * m * (com_mod.gtnNo + 1);
+ for (int j = 0; j < U.ncols(); 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 * (r + 1);
+ 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

@@ -1423,6 +1423,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 +2024,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/main.cpp

@@ -243,6 +243,26 @@ void iterate_solution(Simulation* simulation)
 
  set_bc::set_bc_dir(com_mod, An, Yn, Dn);
 
+ /*
+ if (com_mod.usePrecomp) {
+ #ifdef debug_iterate_solution
+ dmsg << "Use precomputed values ..." << std::endl;
+ #endif
+ double cT = cTS * dt;
+ double precompDt = com_mod.precompDt;
+ int n1 = static_cast<int>(cT / precompDt) + 1;
+ int n2 = n1 + 1;
+ double alpha = (cT - (n1 - 1) * precompDt) / precompDt;
+ for (int l = 0; l < com_mod.nMsh; l++) {
+ auto lM = com_mod.msh[l];
+ for (int i = 0; i < nsd; i++) {
+ for (int j = 0; j < tnNo; j++) {
+ Yn(i, j) = (1.0 - alpha) * lM.Ys(i,j,n1) + alpha * lM.Ys(i, j, n2);
+ }
+ }
+ }
+ }
+ */
  // Inner loop for iteration
  //
  int inner_count = 1;
@@ -706,12 +726,13 @@ int main(int argc, char *argv[])
  #endif
  read_files(simulation, file_name);
 
+ std::cout << "Degrees of Freedom [end read_files]: " << simulation->com_mod.tDof << std::endl;
  // Distribute data to processors.
  #ifdef debug_main
  dmsg << "Distribute data to processors " << " ... ";
  #endif
  distribute(simulation);
-
+ std::cout << "Degrees of Freedom [end distribute]: " << simulation->com_mod.tDof << std::endl;
  // Initialize simulation data.
  //
  Vector<double> init_time(3);
@@ -720,7 +741,7 @@ int main(int argc, char *argv[])
  dmsg << "Initialize " << " ... ";
  #endif
  initialize(simulation, init_time);
-
+ std::cout << "Degrees of Freedom [end initialize]: " << simulation->com_mod.tDof << std::endl;
  #ifdef debug_main
  for (int iM = 0; iM < simulation->com_mod.nMsh; iM++) {
  dmsg << "---------- iM " << iM;

Code/Source/svFSI/read_files.cpp

@@ -1730,11 +1730,15 @@ 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 (general_params.use_precomputed_solution.value()) {
+ std::cout << "Degrees of Freedom: " << com_mod.tDof << std::endl;
+ }
  } 
 
  if (eq.phys == EquationType::phys_mesh) { 
@@ -1812,6 +1816,7 @@ void read_files(Simulation* simulation, const std::string& file_name)
  #ifdef debug_read_files
  dmsg << "Done" << " ";
  #endif
+ std::cout << "Degrees of Freedom [end read_eq]: " << com_mod.tDof << std::endl;
 }
 
 //--------------------------------