Iterate both heat and FEA solvers for thermoelastic simulations

Common/src/CConfig.cpp

@@ -4003,12 +4003,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 
   /*--- Set the number of external iterations to 1 for the steady state problem ---*/
 
-  if (Kind_Solver == MAIN_SOLVER::FEM_ELASTICITY) {
-    nMGLevels = 0;
-    if (Kind_Struct_Solver == STRUCT_DEFORMATION::SMALL){
-      MinLogResidual = log10(Linear_Solver_Error);
-    }
-  }
+  if (Kind_Solver == MAIN_SOLVER::FEM_ELASTICITY) nMGLevels = 0;
 
   Radiation = (Kind_Radiation != RADIATION_MODEL::NONE);
 

SU2_CFD/src/iteration/CFEAIteration.cpp

@@ -40,140 +40,158 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom
 
   const bool nonlinear = (config[val_iZone]->GetGeometricConditions() == STRUCT_DEFORMATION::LARGE);
   const bool linear = (config[val_iZone]->GetGeometricConditions() == STRUCT_DEFORMATION::SMALL);
+  const bool heat = config[val_iZone]->GetWeakly_Coupled_Heat();
   const bool disc_adj_fem = config[val_iZone]->GetDiscrete_Adjoint();
 
   /*--- Loads applied in steps (not used for discrete adjoint). ---*/
   const bool incremental_load = config[val_iZone]->GetIncrementalLoad() && !disc_adj_fem;
 
+  /*--- We need to restore the current inner iter in discrete adjoint cases because file output depends on it. ---*/
+  const auto CurIter = config[val_iZone]->GetInnerIter();
+
+  CGeometry** geo = geometry[val_iZone][val_iInst];
   CIntegration* feaIntegration = integration[val_iZone][val_iInst][FEA_SOL];
   CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL];
 
-  /*--- Add heat solver integration step. ---*/
-  if (config[val_iZone]->GetWeakly_Coupled_Heat()) {
+  const auto nPoint = geo[MESH_0]->GetnPoint();
+  const auto nDim = geo[MESH_0]->GetnDim();
+
+  su2activematrix Coords;
+  if (nonlinear && heat && Coords.empty()) {
+    /*--- Save the mesh coordinates to solve the heat equations in the deformed configuration
+     * and then restore the coordinates for the FEA solver. We could use CoordsOld of CPoint but
+     * that would require accounting for a lot of physics-specific logic in a geometry class. ---*/
+    Coords = geo[MESH_0]->nodes->GetCoord();
+  }
+
+  auto IterateHeat = [&]() {
+    /*--- Update the FVM mesh for the heat solver based on the structural deformations. ---*/
+    if (nonlinear) {
+      SU2_OMP_PARALLEL_(for schedule(static))
+      for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) {
+        for (auto iDim = 0u; iDim < nDim; ++iDim) {
+          geo[MESH_0]->nodes->SetCoord(iPoint, iDim, Coords(iPoint, iDim) + feaSolver->GetNodes()->GetSolution(iPoint, iDim));
+        }
+      }
+      END_SU2_OMP_PARALLEL
+
+      CGeometry::UpdateGeometry(geo, config[val_iZone]);
+    }
+
     config[val_iZone]->SetGlobalParam(MAIN_SOLVER::HEAT_EQUATION, RUNTIME_HEAT_SYS);
     integration[val_iZone][val_iInst][HEAT_SOL]->SingleGrid_Iteration(geometry, solver, numerics, config,
                                                                       RUNTIME_HEAT_SYS, val_iZone, val_iInst);
-  }
 
-  /*--- FEA equations ---*/
-  config[val_iZone]->SetGlobalParam(MAIN_SOLVER::FEM_ELASTICITY, RUNTIME_FEA_SYS);
+    /*--- Restore the coordinates for the FEA solver. ---*/
+    if (nonlinear) {
+      SU2_OMP_PARALLEL_(for schedule(static))
+      for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) {
+        for (auto iDim = 0u; iDim < nDim; ++iDim) {
+          geo[MESH_0]->nodes->SetCoord(iPoint, iDim, Coords(iPoint, iDim));
+        }
+      }
+      END_SU2_OMP_PARALLEL
+    }
+  };
 
-  if (linear) {
-    /*--- Run the (one) iteration ---*/
+  auto Iterate = [&](unsigned long IntIter) {
+    config[val_iZone]->SetInnerIter(IntIter);
 
-    config[val_iZone]->SetInnerIter(0);
+    /*--- Heat transfer equations. ---*/
+    if (heat) IterateHeat();
 
+    /*--- FEA equations. ---*/
+    config[val_iZone]->SetGlobalParam(MAIN_SOLVER::FEM_ELASTICITY, RUNTIME_FEA_SYS);
     feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, val_iInst);
 
     if (!disc_adj_fem) {
-      Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox,
-              val_iZone, INST_0);
-
-      /*--- Set the convergence monitor to true, to prevent the solver to stop in intermediate FSI subiterations ---*/
-      output->SetConvergence(true);
+      StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement,
+                         FFDBox, val_iZone, INST_0);
     }
+  };
 
-  } else if (nonlinear && !incremental_load) {
-    /*--- THIS IS THE DIRECT APPROACH (NO INCREMENTAL LOAD APPLIED) ---*/
-
-    /*--- Keep the current inner iter, we need to restore it in discrete adjoint cases
-     * because file output depends on it. ---*/
-    const auto CurIter = config[val_iZone]->GetInnerIter();
 
+  if (linear || (nonlinear && !incremental_load)) {
+    /*--- THIS IS THE DIRECT APPROACH (NO INCREMENTAL LOAD APPLIED) ---*/
     /*--- Newton-Raphson subiterations ---*/
 
     for (IntIter = 0; IntIter < config[val_iZone]->GetnInner_Iter(); IntIter++) {
-      config[val_iZone]->SetInnerIter(IntIter);
-
-      feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, val_iInst);
+      Iterate(IntIter);
 
       /*--- Limit to only one iteration for the discrete adjoint recording, restore inner iter (see above) ---*/
       if (disc_adj_fem) {
         config[val_iZone]->SetInnerIter(CurIter);
         break;
       }
-      StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement,
-                         FFDBox, val_iZone, INST_0);
-
-      if (StopCalc && (IntIter > 0)) break;
+      /*--- Linear elasticity without thermal effects only needs one iteration. ---*/
+      if (linear && !heat) {
+        output->SetConvergence(true);
+        break;
+      }
+      /*--- Normal stopping criteria. ---*/
+      if (StopCalc && IntIter > 0) break;
     }
+    return;
+  }
 
-  } else {
-    /*--- THIS IS THE INCREMENTAL LOAD APPROACH (only makes sense for nonlinear) ---*/
-
-    /*--- Assume the initial load increment as 1.0 ---*/
+  /*--- THIS IS THE INCREMENTAL LOAD APPROACH (only makes sense for nonlinear) ---*/
 
-    feaSolver->SetLoad_Increment(0, 1.0);
-    feaSolver->SetForceCoeff(1.0);
+  /*--- Assume the initial load increment as 1.0 ---*/
+  feaSolver->SetLoad_Increment(0, 1.0);
+  feaSolver->SetForceCoeff(1.0);
 
-    /*--- Run two nonlinear iterations to check if incremental loading can be skipped ---*/
+  /*--- Run two nonlinear iterations to check if incremental loading can be skipped ---*/
+  for (IntIter = 0; IntIter < 2; ++IntIter) {
+    Iterate(IntIter);
+  }
 
-    auto Iterate = [&](unsigned long IntIter) {
-      config[val_iZone]->SetInnerIter(IntIter);
-      feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, val_iInst);
+  /*--- Early return if we already meet the convergence criteria. ---*/
+  if (StopCalc) return;
 
-      StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement,
-                         FFDBox, val_iZone, INST_0);
-    };
+  /*--- Check user-defined criteria to either increment loads or continue with NR iterations. ---*/
+  bool meetCriteria = true;
+  for (int i = 0; i < 3; ++i) {
+    meetCriteria &= (log10(feaSolver->GetRes_FEM(i)) < config[val_iZone]->GetIncLoad_Criteria(i));
+  }
 
-    for (IntIter = 0; IntIter < 2; ++IntIter) {
+  /*--- If the criteria is met, i.e. the load is not too large, continue the regular calculation. ---*/
+  if (meetCriteria) {
+    for (IntIter = 2; IntIter < config[val_iZone]->GetnInner_Iter(); IntIter++) {
       Iterate(IntIter);
+      if (StopCalc) break;
     }
+    return;
+  }
 
-    /*--- Early return if we already meet the convergence criteria. ---*/
-    if (StopCalc) return;
-
-    /*--- Check user-defined criteria to either increment loads or continue with NR iterations. ---*/
-
-    bool meetCriteria = true;
-    for (int i = 0; i < 3; ++i)
-      meetCriteria &= (log10(feaSolver->GetRes_FEM(i)) < config[val_iZone]->GetIncLoad_Criteria(i));
-
-    /*--- If the criteria is met, i.e. the load is not too large, continue the regular calculation. ---*/
-
-    if (meetCriteria) {
-      /*--- Newton-Raphson subiterations ---*/
-
-      for (IntIter = 2; IntIter < config[val_iZone]->GetnInner_Iter(); IntIter++) {
-        Iterate(IntIter);
-        if (StopCalc) break;
-      }
-
-    }
-
-    /*--- If the criteria is not met, a whole set of subiterations for the different loads must be done. ---*/
-
-    else {
-      /*--- Restore solution to initial. Because we ramp the load from zero, in multizone cases it is not
-       * adequate to take "old values" as those will be for maximum loading on the previous outer iteration. ---*/
+  /*--- If the criteria is not met, a whole set of subiterations for the different loads must be done.
+   * Restore solution to initial. Because we ramp the load from zero, in multizone cases it is not
+   * adequate to take "old values" as those will be for maximum loading on the previous outer iteration. ---*/
 
-      feaSolver->SetInitialCondition(geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone],
-                                     TimeIter);
+  feaSolver->SetInitialCondition(geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone],
+                                 TimeIter);
 
-      /*--- For the number of increments ---*/
-      for (auto iIncrement = 1ul; iIncrement <= nIncrements; iIncrement++) {
-        /*--- Set the load increment and the initial condition, and output the
-         *    parameters of UTOL, RTOL, ETOL for the previous iteration ---*/
+  /*--- For the number of increments ---*/
+  for (auto iIncrement = 1ul; iIncrement <= nIncrements; iIncrement++) {
+    /*--- Set the load increment and the initial condition, and output the
+     * parameters of UTOL, RTOL, ETOL for the previous iteration ---*/
 
-        su2double loadIncrement = su2double(iIncrement) / nIncrements;
-        feaSolver->SetLoad_Increment(iIncrement, loadIncrement);
+    su2double loadIncrement = su2double(iIncrement) / nIncrements;
+    feaSolver->SetLoad_Increment(iIncrement, loadIncrement);
 
-        /*--- Set the convergence monitor to false, to force the solver to converge every subiteration ---*/
-        output->SetConvergence(false);
+    /*--- Set the convergence monitor to false, to force the solver to converge every subiteration. ---*/
+    output->SetConvergence(false);
 
-        if (rank == MASTER_NODE) cout << "\nIncremental load: increment " << iIncrement << endl;
+    if (rank == MASTER_NODE) cout << "\nIncremental load: increment " << iIncrement << endl;
 
-        /*--- Newton-Raphson subiterations ---*/
+    /*--- Newton-Raphson subiterations ---*/
 
-        for (IntIter = 0; IntIter < config[val_iZone]->GetnInner_Iter(); IntIter++) {
-          Iterate(IntIter);
-          if (StopCalc && (IntIter > 0)) break;
-        }
-      }
-      /*--- Just to be sure, set default increment settings. ---*/
-      feaSolver->SetLoad_Increment(0, 1.0);
+    for (IntIter = 0; IntIter < config[val_iZone]->GetnInner_Iter(); IntIter++) {
+      Iterate(IntIter);
+      if (StopCalc && IntIter > 0) break;
     }
   }
+  /*--- Just to be sure, set default increment settings. ---*/
+  feaSolver->SetLoad_Increment(0, 1.0);
 }
 
 void CFEAIteration::Update(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver,

TestCases/fea_fsi/ThermalBeam_3d/configBeamNonlinear_3d.cfg

@@ -0,0 +1,56 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% SU2 configuration file                                                       %
+% Case description: 3D beam with thermal expansion                             %
+% File Version 8.4.0 "Harrier"                                                 %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+SOLVER= ELASTICITY
+MATH_PROBLEM= DIRECT
+GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS
+MATERIAL_MODEL= NEO_HOOKEAN
+MESH_FILENAME= ../StatBeam_3d/meshBeam_3d.su2
+ELASTICITY_MODULUS= 3E7
+POISSON_RATIO= 0.3
+MATERIAL_THERMAL_EXPANSION_COEFF= 2e-5
+MATERIAL_REFERENCE_TEMPERATURE= 288.15
+MATERIAL_DENSITY= 7854
+MARKER_CLAMPED= ( left, right )
+MARKER_PRESSURE= ( lower, 0, symleft, 0, symright, 0 )
+MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0 )
+LINEAR_SOLVER= FGCRODR
+LINEAR_SOLVER_PREC= ILU
+LINEAR_SOLVER_ERROR= 1E-4
+LINEAR_SOLVER_ITER= 1000
+LINEAR_SOLVER_RESTART_FREQUENCY= 50
+LINEAR_SOLVER_RESTART_DEFLATION= 10
+MESH_FORMAT= SU2
+TABULAR_FORMAT= CSV
+CONV_FILENAME= history_beam
+VOLUME_FILENAME= beam
+RESTART_FILENAME= restart_beam
+SOLUTION_FILENAME= restart_beam
+OUTPUT_WRT_FREQ= 1
+
+CONV_FIELD= REL_RMS_RTOL, REL_RMS_TEMPERATURE
+CONV_STARTITER= 0
+CONV_RESIDUAL_MINVAL= -5
+INNER_ITER= 10
+
+% Coupling with heat solver.
+WEAKLY_COUPLED_HEAT_EQUATION= YES
+FREESTREAM_TEMPERATURE= 300
+SPECIFIC_HEAT_CP= 460
+THERMAL_CONDUCTIVITY_CONSTANT= 45
+
+% Copy markers to allow specifying boundary conditions for both solvers on the
+% same markers.
+MARKER_CREATE_COPY= ( left, left_heat, right, right_heat )
+MARKER_ISOTHERMAL= ( left_heat, 400, right_heat, 300 )
+
+NUM_METHOD_GRAD= GREEN_GAUSS
+TIME_DISCRE_HEAT= EULER_IMPLICIT
+CFL_NUMBER= 1e8
+
+MARKER_MONITORING= ( left_heat )
+SCREEN_OUTPUT= INNER_ITER, RMS_RES, LINSOL, VMS, TOTAL_HEATFLUX
+

TestCases/fea_fsi/ThermalBeam_3d/configBeam_3d.cfg

@@ -9,26 +9,30 @@ MATH_PROBLEM= DIRECT
 GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS
 MATERIAL_MODEL= LINEAR_ELASTIC
 MESH_FILENAME= ../StatBeam_3d/meshBeam_3d.su2
-ELASTICITY_MODULUS=3E7
-POISSON_RATIO=0.3
+ELASTICITY_MODULUS= 3E7
+POISSON_RATIO= 0.3
 MATERIAL_THERMAL_EXPANSION_COEFF= 2e-5
 MATERIAL_REFERENCE_TEMPERATURE= 288.15
-MATERIAL_DENSITY=7854
+MATERIAL_DENSITY= 7854
 MARKER_CLAMPED= ( left, right )
 MARKER_PRESSURE= ( lower, 0, symleft, 0, symright, 0 )
 MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0 )
 LINEAR_SOLVER= CONJUGATE_GRADIENT
 LINEAR_SOLVER_PREC= ILU
-LINEAR_SOLVER_ERROR= 1E-8
-LINEAR_SOLVER_ITER= 1000
+LINEAR_SOLVER_ERROR= 1E-4
+LINEAR_SOLVER_ITER= 500
 MESH_FORMAT= SU2
 TABULAR_FORMAT= CSV
 CONV_FILENAME= history_beam
 VOLUME_FILENAME= beam
 RESTART_FILENAME= restart_beam
 SOLUTION_FILENAME= restart_beam
 OUTPUT_WRT_FREQ= 1
-INNER_ITER= 1
+
+CONV_FIELD= REL_RMS_DISP_X, REL_RMS_DISP_Y, REL_RMS_DISP_Z, REL_RMS_TEMPERATURE
+CONV_STARTITER= 0
+CONV_RESIDUAL_MINVAL= -4
+INNER_ITER= 10
 
 % Coupling with heat solver.
 WEAKLY_COUPLED_HEAT_EQUATION= YES

TestCases/parallel_regression.py

@@ -1276,12 +1276,20 @@ def main():
     thermal_beam_3d = TestCase('thermal_beam_3d')
     thermal_beam_3d.cfg_dir = "fea_fsi/ThermalBeam_3d"
     thermal_beam_3d.cfg_file = "configBeam_3d.cfg"
-    thermal_beam_3d.test_iter = 0
-    thermal_beam_3d.test_vals = [-6.140220, -5.842734, -5.972391, -8.091358, 262, -8.246755, 81, -8.298569, 135620, 144.65]
-    thermal_beam_3d.test_vals_aarch64 = [-6.131860, -5.845164, -5.964084, -8.091358, 262, -8.244325, 81, -8.298569, 135620, 144.65]
+    thermal_beam_3d.test_iter = 4
+    thermal_beam_3d.test_vals = [-8.137623, -7.840833, -7.903285, -13.978110, 217, -4.093241, 39, -4.072614, 1.3676e+05, 75.0]
     thermal_beam_3d.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f")
     test_list.append(thermal_beam_3d)
 
+    # Static beam, 3d with coupled temperature, nonlinear elasticity
+    thermal_beam_nl_3d = TestCase('thermal_beam_nl_3d')
+    thermal_beam_nl_3d.cfg_dir = "fea_fsi/ThermalBeam_3d"
+    thermal_beam_nl_3d.cfg_file = "configBeamNonlinear_3d.cfg"
+    thermal_beam_nl_3d.test_iter = 8
+    thermal_beam_nl_3d.test_vals = [-7.564309, -2.992893, -12.242503, -14.068322, 57, -4.017665, 24, -4.204804, 138710, 75.233]
+    thermal_beam_nl_3d.command = TestCase.Command(exec = "parallel_computation_fsi.py", param = "-f")
+    test_list.append(thermal_beam_nl_3d)
+
     # Rotating cylinder, 3d
     rotating_cylinder_fea           = TestCase('rotating_cylinder_fea')
     rotating_cylinder_fea.cfg_dir   = "fea_fsi/rotating_cylinder"