Add improvements, bug fixes, + notes from j-signorelli/su2-precice-pywrapper (#34)

* Nondimensional mode implementation
* Fix compilation w/ OMP bug
* Specify problem dimension as argument
* Clarify meaning of Monitor
* Save time and iter for implicit coupling
* Fix CHT initial read/write (was misleading) (read_data --> write_data)
* Add disclaimer/note on restarting simulations

Author: j-signorelli

README.md

@@ -166,7 +166,11 @@ mpirun -n 8 python3 SU2_preCICE_CHT.py -f SU2_config_file.cfg --parallel
 
 **NOTE**: As of SU2 v7.5.1: Deforming `MARKER_EULER`'s are buggy when simulations are run in parallel, leading to unexpected results. More information can be found at this discussion here: https://github.com/su2code/SU2/discussions/1931.
 
-### Further notes
+## Important note on restarts
+
+This code **has not been tested** for restarts using initializations *from* SU2. Any restarted simulations should have SU2 be the first participant and receive initialization data. It is possible that, if SU2 must send initialization data, that it is incorrect (it may use default values in the config file, or just be zeros if the data hasn't been computed until after/during a first iteration). Admittedly, this is from a lack of understanding of the specifics of how SU2 operates and there may not be a trivial work-around.
+
+## Further notes
 
 Result files (vtu) generated from SU2 might be incompatible with your ParaView version. For example, ParaView 5.11.2 on Ubuntu 22.04 is known to fail with SU2 7.5.1 result files, but ParaView 5.12 works.
 

replacement_files/python_wrapper_structure.cpp

@@ -47,7 +47,41 @@ void CDriver::PythonInterface_Preprocessing(CConfig **config, CGeometry ****geom
       if (rank == MASTER_NODE) cout << "Setting customized boundary conditions for zone " << iZone << endl;
       for (iMesh = 0; iMesh <= config[iZone]->GetnMGLevels(); iMesh++) {
         geometry[iZone][INST_0][iMesh]->SetCustomBoundary(config[iZone]);
+        // preCICE: Update input file initial custom boundary temperature/HF to be nondimensional.
+        //          Was fixed in SU2 v8 here: https://github.com/su2code/SU2/pull/2078
+        //------------------------------------------------------------------------------------------------
+        // Basically just copied + pasted (terrible coding but as a quick/least invasive fix) CGeometry::SetCustomBoundary
+        unsigned short iMarker;
+        unsigned long iVertex;
+        string Marker_Tag;
+
+        for(iMarker=0; iMarker < geometry[iZone][INST_0][iMesh]->GetnMarker(); iMarker++){
+          Marker_Tag = config[iZone]->GetMarker_All_TagBound(iMarker);
+          if(config[iZone]->GetMarker_All_PyCustom(iMarker)){
+            switch(config[iZone]->GetMarker_All_KindBC(iMarker)){
+              case HEAT_FLUX:
+                for(iVertex=0; iVertex < geometry[iZone][INST_0][iMesh]->GetnVertex(iMarker); iVertex++){
+                  geometry[iZone][INST_0][iMesh]->SetCustomBoundaryHeatFlux(iMarker, iVertex, config[iZone]->GetWall_HeatFlux(Marker_Tag)/config[iZone]->GetHeat_Flux_Ref());
+                }
+                break;
+              case ISOTHERMAL:
+                for(iVertex=0; iVertex < geometry[iZone][INST_0][iMesh]->GetnVertex(iMarker); iVertex++){
+                  geometry[iZone][INST_0][iMesh]->SetCustomBoundaryTemperature(iMarker, iVertex, config[iZone]->GetIsothermal_Temperature(Marker_Tag)/config[iZone]->GetTemperature_Ref());
+                }
+                break;
+              case INLET_FLOW:
+                // This case is handled in the solver class.
+                break;
+              default:
+                cout << "WARNING: Marker " << Marker_Tag << " is not customizable. Using default behavior." << endl;
+                break;
+            }
+          }
+        }
+          //-------------------------------------------------------------------------------------------------
       }
+
+
       geometry[iZone][INST_0][MESH_0]->UpdateCustomBoundaryConditions(geometry[iZone][INST_0], config[iZone]);
 
       if ((config[iZone]->GetKind_Solver() == MAIN_SOLVER::EULER) ||
@@ -491,7 +525,7 @@ void CDriver::FinalizeMESH_SOL() {
   const bool secondOrder = config_container[ZONE_0]->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND;
   const su2double invTimeStep = 1.0 / config_container[ZONE_0]->GetDelta_UnstTimeND();
 
-  SU2_OMP_FOR_STAT(omp_chunk_size)
+  //SU2_OMP_FOR_STAT(omp_chunk_size) - omp_chunk_size part of CMeshSolver
   for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
 
     /*--- Coordinates of the current point at n+1, n, & n-1 time levels. ---*/
@@ -513,7 +547,7 @@ void CDriver::FinalizeMESH_SOL() {
       geometry_container[ZONE_0][INST_0][MESH_0]->nodes->SetGridVel(iPoint, iDim, GridVel);
     }
   }
-  END_SU2_OMP_FOR
+  //END_SU2_OMP_FOR
 
   for (auto iMGlevel = 1u; iMGlevel <= config_container[ZONE_0]->GetnMGLevels(); iMGlevel++)
     geometry_container[ZONE_0][INST_0][iMGlevel]->SetRestricted_GridVelocity(geometry_container[ZONE_0][INST_0][iMGlevel-1]);
@@ -635,13 +669,15 @@ passivedouble CDriver::GetVertexTemperature(unsigned short iMarker, unsigned lon
     vertexWallTemp = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetTemperature(iPoint);
   }
 
-  return SU2_TYPE::GetValue(vertexWallTemp);
+  //preCICE: re-dimensionalize before returning
+  return SU2_TYPE::GetValue(vertexWallTemp * config_container[ZONE_0]->GetTemperature_Ref());
 
 }
 
 void CDriver::SetVertexTemperature(unsigned short iMarker, unsigned long iVertex, passivedouble val_WallTemp){
 
-  geometry_container[ZONE_0][INST_0][MESH_0]->SetCustomBoundaryTemperature(iMarker, iVertex, val_WallTemp);
+  // preCICE: non-dimensionalize before setting
+  geometry_container[ZONE_0][INST_0][MESH_0]->SetCustomBoundaryTemperature(iMarker, iVertex, val_WallTemp / config_container[ZONE_0]->GetTemperature_Ref());
 }
 
 vector<passivedouble> CDriver::GetVertexHeatFluxes(unsigned short iMarker, unsigned long iVertex) const {
@@ -671,9 +707,10 @@ vector<passivedouble> CDriver::GetVertexHeatFluxes(unsigned short iMarker, unsig
     }
   }
 
-  HeatFluxPassive[0] = SU2_TYPE::GetValue(HeatFlux[0]);
-  HeatFluxPassive[1] = SU2_TYPE::GetValue(HeatFlux[1]);
-  HeatFluxPassive[2] = SU2_TYPE::GetValue(HeatFlux[2]);
+  //preCICE: re-dimensionalize before returning
+  HeatFluxPassive[0] = SU2_TYPE::GetValue(HeatFlux[0] * config_container[ZONE_0]->GetHeat_Flux_Ref());
+  HeatFluxPassive[1] = SU2_TYPE::GetValue(HeatFlux[1] * config_container[ZONE_0]->GetHeat_Flux_Ref());
+  HeatFluxPassive[2] = SU2_TYPE::GetValue(HeatFlux[2] * config_container[ZONE_0]->GetHeat_Flux_Ref());
 
   return HeatFluxPassive;
 }
@@ -718,12 +755,14 @@ passivedouble CDriver::GetVertexNormalHeatFlux(unsigned short iMarker, unsigned
     vertexWallHeatFlux = -thermal_conductivity*dTdn;
   }
 
-  return SU2_TYPE::GetValue(vertexWallHeatFlux);
+   //preCICE: re-dimensionalize before returning
+  return SU2_TYPE::GetValue(vertexWallHeatFlux * config_container[ZONE_0]->GetHeat_Flux_Ref());
 }
 
 void CDriver::SetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVertex, passivedouble val_WallHeatFlux){
 
-  geometry_container[ZONE_0][INST_0][MESH_0]->SetCustomBoundaryHeatFlux(iMarker, iVertex, val_WallHeatFlux);
+  // preCICE: non-dimensionalize before setting
+  geometry_container[ZONE_0][INST_0][MESH_0]->SetCustomBoundaryHeatFlux(iMarker, iVertex, val_WallHeatFlux / config_container[ZONE_0]->GetHeat_Flux_Ref());
 }
 
 passivedouble CDriver::GetThermalConductivity(unsigned short iMarker, unsigned long iVertex) const {

run/SU2_preCICE_CHT.py

@@ -55,9 +55,11 @@ def main():
   parser.add_option("-c", "--precice-config", dest="precice_config", help="Specify preCICE config file", default="../precice-config.xml")
   parser.add_option("-m", "--precice-mesh", dest="precice_mesh", help="Specify the preCICE mesh name", default="Fluid-Mesh")
   parser.add_option("-r", "--precice-reverse", action="store_true", dest="precice_reverse", help="Include flag to have SU2 write temperature, read heat flux", default=False)
-
+  
+  # Dimension
+  parser.add_option("-d", "--dimension", dest="nDim", help="Dimension of fluid domain", type="int", default=3)
+  
   (options, args) = parser.parse_args()
-  options.nDim = int(2) # Specify dimension here
   options.nZone = int(1) # Specify number of zones here (1)
 
   # Import mpi4py for parallel run
@@ -177,9 +179,9 @@ def main():
   if (interface.is_action_required(precice.action_write_initial_data())):
 
     for i, iVertex in enumerate(iVertices_CHTMarker_PHYS):
-      read_data[i] = GetInitialFxn(CHTMarkerID, iVertex)
+      write_data[i] = GetInitialFxn(CHTMarkerID, iVertex)
 
-    interface.write_block_scalar_data(write_data_id, vertex_ids, read_data)
+    interface.write_block_scalar_data(write_data_id, vertex_ids, write_data)
     interface.mark_action_fulfilled(precice.action_write_initial_data())
 
   interface.initialize_data()
@@ -194,13 +196,16 @@ def main():
   if options.with_MPI == True:
     comm.Barrier()
 
-
+  precice_saved_time = 0
+  precice_saved_iter = 0
   while (interface.is_coupling_ongoing()):
 
     # Implicit coupling
     if (interface.is_action_required(precice.action_write_iteration_checkpoint())):
       # Save the state
       SU2Driver.SaveOldState()
+      precice_saved_time = time
+      precice_saved_iter = TimeIter
       interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
 
     if (interface.is_read_data_available()):
@@ -234,7 +239,7 @@ def main():
     # Update the solver for the next time iteration
     SU2Driver.Update()
 
-    # Monitor the solver and output solution to file if required
+    # Monitor the solver
     stopCalc = SU2Driver.Monitor(TimeIter)
 
     if (interface.is_write_data_required(deltaT)):
@@ -253,6 +258,8 @@ def main():
     if (interface.is_action_required(precice.action_read_iteration_checkpoint())):
       # Reload old state
       SU2Driver.ReloadOldState()
+      time = precice_saved_time
+      TimeIter = precice_saved_iter
       interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
     else: # Output and increment as usual
       SU2Driver.Output(TimeIter)

run/SU2_preCICE_FSI.py

@@ -55,8 +55,10 @@ def main():
     parser.add_option("-c", "--precice-config", dest="precice_config", help="Specify preCICE config file", default="../precice-config.xml")
     parser.add_option("-m", "--precice-mesh", dest="precice_mesh", help="Specify the preCICE mesh name", default="Fluid-Mesh")
 
+    # Dimension
+    parser.add_option("-d", "--dimension", dest="nDim", help="Dimension of fluid domain", type="int", default=3)
+  
     (options, args) = parser.parse_args()
-    options.nDim  = int(2)
     options.nZone = int(1)
 
     # Import mpi4py for parallel run
@@ -180,12 +182,16 @@ def main():
     if options.with_MPI == True:
         comm.Barrier()
 
+    precice_saved_time = 0
+    precice_saved_iter = 0
     while (interface.is_coupling_ongoing()):#(TimeIter < nTimeIter):
 
         # Implicit coupling
         if (interface.is_action_required(precice.action_write_iteration_checkpoint())):
             # Save the state
             SU2Driver.SaveOldState()
+            precice_saved_time = time
+            precice_saved_iter = TimeIter
             interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
 
         if (interface.is_read_data_available()):
@@ -220,7 +226,7 @@ def main():
         # Update the solver for the next time iteration
         SU2Driver.Update()
 
-        # Monitor the solver and output solution to file if required
+        # Monitor the solver
         stopCalc = SU2Driver.Monitor(TimeIter)
 
 
@@ -241,6 +247,8 @@ def main():
         if (interface.is_action_required(precice.action_read_iteration_checkpoint())):
             # Reload old state
             SU2Driver.ReloadOldState()
+            time = precice_saved_time
+            TimeIter = precice_saved_iter
             interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
         else: # Output and increment as usual
             SU2Driver.Output(TimeIter)