From d5b1be59c585db3553b2332c1f50d9f17082c63b Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Wed, 27 Jan 2021 18:40:45 +0100 Subject: [PATCH 01/32] Temporary fix to restart waiting for complete PR --- SU2_CFD/src/solvers/CEulerSolver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 1629d078f7e6..76c2c05622c2 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -9287,7 +9287,7 @@ void CEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig /*--- For dynamic meshes, read in and store the grid coordinates and grid velocities for each node. ---*/ - if (dynamic_grid && val_update_geo) { + if (dynamic_grid && val_update_geo && !config->GetDeform_Mesh()) { /*--- Read in the next 2 or 3 variables which are the grid velocities ---*/ /*--- If we are restarting the solution from a previously computed static calculation (no grid movement) ---*/ @@ -9313,7 +9313,7 @@ void CEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig /*--- For static FSI problems, grid_movement is 0 but we need to read in and store the grid coordinates for each node (but not the grid velocities, as there are none). ---*/ - if (static_fsi && val_update_geo) { + if (static_fsi && val_update_geo && !config->GetDeform_Mesh()) { /*--- Rewind the index to retrieve the Coords. ---*/ index = counter*Restart_Vars[1]; Coord = &Restart_Data[index]; From 614d009710076cee88b52815e2cfdaf4943b6351 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 8 Feb 2021 20:48:45 +0100 Subject: [PATCH 02/32] Tentative removing switch from python --- SU2_PY/FSI_tools/FSI_config.py | 57 ++++++----- SU2_PY/FSI_tools/__init__.py | 1 - SU2_PY/FSI_tools/switch.py | 52 ---------- SU2_PY/SU2_Nastran/pysu2_nastran.py | 145 ++++++++++++++-------------- 4 files changed, 97 insertions(+), 158 deletions(-) delete mode 100644 SU2_PY/FSI_tools/switch.py diff --git a/SU2_PY/FSI_tools/FSI_config.py b/SU2_PY/FSI_tools/FSI_config.py index f9bba78037d5..0172e8fabfea 100644 --- a/SU2_PY/FSI_tools/FSI_config.py +++ b/SU2_PY/FSI_tools/FSI_config.py @@ -39,8 +39,10 @@ # Imports # ---------------------------------------------------------------------- -import os, sys, shutil, copy -from FSI_tools.switch import switch +import os +import sys +import shutil +import copy # ---------------------------------------------------------------------- # FSI Configuration Class @@ -86,40 +88,35 @@ def readConfig(self): this_param = line[0].strip() this_value = line[1].strip() - for case in switch(this_param): #integer values - if case("NDIM") : pass - if case("RESTART_ITER") : pass - if case("TIME_TRESHOLD") : pass - if case("NB_FSI_ITER") : - self._ConfigContent[this_param] = int(this_value) - break + if (this_param == "NDIM") || \ + (this_param == "RESTART_ITER") || \ + (this_param == "TIME_TRESHOLD") || \ + (this_param == "NB_FSI_ITER") : + self._ConfigContent[this_param] = int(this_value) #float values - if case("RBF_RADIUS") : pass - if case("AITKEN_PARAM") : pass - if case("UNST_TIMESTEP") : pass - if case("UNST_TIME") : pass - if case("FSI_TOLERANCE") : - self._ConfigContent[this_param] = float(this_value) - break + elif (this_param == "RBF_RADIUS") || \ + (this_param == "AITKEN_PARAM") || \ + (this_param == "UNST_TIMESTEP") || \ + (this_param == "UNST_TIME") || \ + (this_param == "FSI_TOLERANCE") : + self._ConfigContent[this_param] = float(this_value) #string values - if case("CFD_CONFIG_FILE_NAME") : pass - if case("CSD_SOLVER") : pass - if case("CSD_CONFIG_FILE_NAME") : pass - if case("RESTART_SOL") : pass - if case("MATCHING_MESH") : pass - if case("MESH_INTERP_METHOD") : pass - if case("DISP_PRED") : pass - if case("AITKEN_RELAX") : pass - if case("TIME_MARCHING") : - self._ConfigContent[this_param] = this_value - break + elif (this_param == "CFD_CONFIG_FILE_NAME") || \ + (this_param == "CSD_SOLVER") || \ + (this_param == "CSD_CONFIG_FILE_NAME") || \ + (this_param == "RESTART_SOL") || \ + (this_param == "MATCHING_MESH") || \ + (this_param == "MESH_INTERP_METHOD") || \ + (this_param == "DISP_PRED") || \ + (this_param == "AITKEN_RELAX") || \ + (this_param == "TIME_MARCHING") : + self._ConfigContent[this_param] = this_value - if case(): - print(this_param + " is an invalid option !") - break + else : + print(this_param + " is an invalid option !") def applyDefaults(self): if self._ConfigContent["CSD_SOLVER"] == "IMPOSED": diff --git a/SU2_PY/FSI_tools/__init__.py b/SU2_PY/FSI_tools/__init__.py index e843f12bbe49..ec8ad1b64101 100644 --- a/SU2_PY/FSI_tools/__init__.py +++ b/SU2_PY/FSI_tools/__init__.py @@ -1,3 +1,2 @@ from FSI_tools.FSIInterface import Interface -from FSI_tools.switch import switch from FSI_tools.FSI_config import FSIConfig diff --git a/SU2_PY/FSI_tools/switch.py b/SU2_PY/FSI_tools/switch.py deleted file mode 100644 index b42eaf6e9ddd..000000000000 --- a/SU2_PY/FSI_tools/switch.py +++ /dev/null @@ -1,52 +0,0 @@ -# ------------------------------------------------------------------- -# Switch Class -# ------------------------------------------------------------------- -# source: Brian Beck, PSF License, ActiveState Code -# http://code.activestate.com/recipes/410692/ - -class switch(object): - """ Readable switch construction - - Example: - - c = 'z' - for case in switch(c): - if case('a'): pass # only necessary if the rest of the suite is empty - if case('b'): pass - # ... - if case('y'): pass - if case('z'): - print("c is lowercase!") - break - if case('A'): pass - # ... - if case('Z'): - print("c is uppercase!") - break - if case(): # default - print("I dunno what c was!") - - source: Brian Beck, PSF License, ActiveState Code - http://code.activestate.com/recipes/410692/ - """ - - def __init__(self, value): - self.value = value - self.fall = False - - def __iter__(self): - """Return the match method once, then stop""" - yield self.match - raise StopIteration - - def match(self, *args): - """Indicate whether or not to enter a case suite""" - if self.fall or not args: - return True - elif self.value in args: - self.fall = True - return True - else: - return False - -#: class switch() diff --git a/SU2_PY/SU2_Nastran/pysu2_nastran.py b/SU2_PY/SU2_Nastran/pysu2_nastran.py index d3a82503448f..53226db3bc1c 100644 --- a/SU2_PY/SU2_Nastran/pysu2_nastran.py +++ b/SU2_PY/SU2_Nastran/pysu2_nastran.py @@ -29,12 +29,13 @@ # Imports # ---------------------------------------------------------------------- -import os, shutil, copy +import os +import shutil +import copy import numpy as np import scipy as sp import scipy.linalg as linalg from math import * -from FSI_tools.switch import switch # ---------------------------------------------------------------------- # Config class @@ -45,60 +46,56 @@ class ImposedMotionFunction: def __init__(self,time0,tipo,parameters): self.time0 = time0 self.tipo = tipo - for case in switch(self.tipo): - if case("SINUSOIDAL"): - self.bias = parameters[0] - self.amplitude = parameters[1] - self.frequency = parameters[2] - break - if case("BLENDED_STEP"): - self.kmax = parameters[0] - self.vinf = parameters[1] - self.lref = parameters[2] - self.amplitude = parameters[3] - self.tmax = 2*pi/self.kmax*self.lref/self.vinf - self.omega0 = 1/2*self.kmax - break - if case(): - raise Exception('Imposed function {} not found, please implement it in pysu2_nastran.py'.format(self.tipo)) - break + if self.tipo == "SINUSOIDAL": + self.bias = parameters[0] + self.amplitude = parameters[1] + self.frequency = parameters[2] + + elif self.tipo == "BLENDED_STEP": + self.kmax = parameters[0] + self.vinf = parameters[1] + self.lref = parameters[2] + self.amplitude = parameters[3] + self.tmax = 2*pi/self.kmax*self.lref/self.vinf + self.omega0 = 1/2*self.kmax + + else: + raise Exception('Imposed function {} not found, please implement it in pysu2_nastran.py'.format(self.tipo)) + def GetDispl(self,time): time = time - self.time0 - for case in switch(self.tipo): - if case("SINUSOIDAL"): - return self.bias+self.amplitude*sin(2*pi*self.frequency*time) - break - if case("BLENDED_STEP"): - if time < self.tmax: - return self.amplitude/2.0*(1.0-cos(self.omega0*time*self.vinf/self.lref)) - return self.amplitude - break + if self.tipo == "SINUSOIDAL": + return self.bias+self.amplitude*sin(2*pi*self.frequency*time) + + if self.tipo == "BLENDED_STEP": + if time < self.tmax: + return self.amplitude/2.0*(1.0-cos(self.omega0*time*self.vinf/self.lref)) + return self.amplitude + def GetVel(self,time): time = time - self.time0 - for case in switch(self.tipo): - if case("SINUSOIDAL"): - return self.amplitude*cos(2*pi*self.frequency*time)*2*pi*self.frequency - break - if case("BLENDED_STEP"): - if time < self.tmax: - return self.amplitude/2.0*sin(self.omega0*time*self.vinf/self.lref)*(self.omega0*self.vinf/self.lref) - return 0.0 - break + + if self.tipo == "SINUSOIDAL": + return self.amplitude*cos(2*pi*self.frequency*time)*2*pi*self.frequency + + if self.tipo == "BLENDED_STEP": + if time < self.tmax: + return self.amplitude/2.0*sin(self.omega0*time*self.vinf/self.lref)*(self.omega0*self.vinf/self.lref) + return 0.0 def GetAcc(self,time): time = time - self.time0 - for case in switch(self.tipo): - if case("SINUSOIDAL"): - return -self.amplitude*sin(2*pi*self.frequency*time)*(2*pi*self.frequency)**2 - break - if case("BLENDED_STEP"): - if time < self.tmax: - return self.amplitude/2.0*cos(self.omega0*time*self.vinf/self.lref)*(self.omega0*self.vinf/self.lref)**2 - return 0.0 - break + + if self.tipo == "SINUSOIDAL": + return -self.amplitude*sin(2*pi*self.frequency*time)*(2*pi*self.frequency)**2 + + if self.tipo == "BLENDED_STEP": + if time < self.tmax: + return self.amplitude/2.0*cos(self.omega0*time*self.vinf/self.lref)*(self.omega0*self.vinf/self.lref)**2 + return 0.0 class RefSystem: @@ -335,39 +332,37 @@ def __readConfig(self): this_param = line[0].strip() this_value = line[1].strip() - for case in switch(this_param): - #integer values - if case("NMODES") : pass - if case("RESTART_ITER") : - self.Config[this_param] = int(this_value) - break + #integer values + if (this_param == "NMODES") || \ + (this_param == "RESTART_ITER": + self.Config[this_param] = int(this_value) - #float values - if case("DELTA_T") : pass - if case("MODAL_DAMPING") : pass - if case("RHO") : - self.Config[this_param] = float(this_value) - break - #string values - if case("TIME_MARCHING") : pass - if case("MESH_FILE") : pass - if case("PUNCH_FILE") : pass - if case("RESTART_SOL") : pass - if case("MOVING_MARKER") : - self.Config[this_param] = this_value - break + #float values + elif (this_param == "DELTA_T") || \ + (this_param == "MODAL_DAMPING") || \ + (this_param == "RHO"): + self.Config[this_param] = float(this_value) - #lists values - if case("INITIAL_MODES"): pass - if case("IMPOSED_MODES"): pass - if case("IMPOSED_PARAMETERS"): - self.Config[this_param] = eval(this_value) - break - if case(): - raise Exception('{} is an invalid option !'.format(this_param)) - break + #string values + elif (this_param == "TIME_MARCHING") || \ + (this_param == "MESH_FILE") || \ + (this_param == "PUNCH_FILE") || \ + (this_param == "RESTART_SOL") || \ + (this_param == "MOVING_MARKER"): + self.Config[this_param] = this_value + + + #lists values + elif (this_param == "INITIAL_MODES") || \ + (this_param == "IMPOSED_MODES") || \ + (this_param == "IMPOSED_PARAMETERS"): + self.Config[this_param] = eval(this_value) + + + else: + raise Exception('{} is an invalid option !'.format(this_param)) From 8a5bfa0996b9ec9405e3583d546b637089facd32 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 8 Feb 2021 21:02:11 +0100 Subject: [PATCH 03/32] First fix to removed switch --- SU2_PY/meson.build | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/SU2_PY/meson.build b/SU2_PY/meson.build index 1deb38390c6c..38136b25db98 100644 --- a/SU2_PY/meson.build +++ b/SU2_PY/meson.build @@ -66,8 +66,7 @@ install_data(['SU2/util/bunch.py', install_data(['FSI_tools/__init__.py', 'FSI_tools/FSIInterface.py', - 'FSI_tools/FSI_config.py', - 'FSI_tools/switch.py'], + 'FSI_tools/FSI_config.py'], install_dir: join_paths(get_option('bindir'), 'FSI_tools')) install_data(['SU2_Nastran/__init__.py', From 750d4df70a158ad47b1a421693bf2dfcf9213b00 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 8 Feb 2021 21:05:44 +0100 Subject: [PATCH 04/32] fixed or --- SU2_PY/FSI_tools/FSI_config.py | 30 ++++++++++++++--------------- SU2_PY/SU2_Nastran/pysu2_nastran.py | 18 ++++++++--------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/SU2_PY/FSI_tools/FSI_config.py b/SU2_PY/FSI_tools/FSI_config.py index 0172e8fabfea..ec8e3c55f34d 100644 --- a/SU2_PY/FSI_tools/FSI_config.py +++ b/SU2_PY/FSI_tools/FSI_config.py @@ -89,29 +89,29 @@ def readConfig(self): this_value = line[1].strip() #integer values - if (this_param == "NDIM") || \ - (this_param == "RESTART_ITER") || \ - (this_param == "TIME_TRESHOLD") || \ + if (this_param == "NDIM") or \ + (this_param == "RESTART_ITER") or \ + (this_param == "TIME_TRESHOLD") or \ (this_param == "NB_FSI_ITER") : self._ConfigContent[this_param] = int(this_value) #float values - elif (this_param == "RBF_RADIUS") || \ - (this_param == "AITKEN_PARAM") || \ - (this_param == "UNST_TIMESTEP") || \ - (this_param == "UNST_TIME") || \ + elif (this_param == "RBF_RADIUS") or \ + (this_param == "AITKEN_PARAM") or \ + (this_param == "UNST_TIMESTEP") or \ + (this_param == "UNST_TIME") or \ (this_param == "FSI_TOLERANCE") : self._ConfigContent[this_param] = float(this_value) #string values - elif (this_param == "CFD_CONFIG_FILE_NAME") || \ - (this_param == "CSD_SOLVER") || \ - (this_param == "CSD_CONFIG_FILE_NAME") || \ - (this_param == "RESTART_SOL") || \ - (this_param == "MATCHING_MESH") || \ - (this_param == "MESH_INTERP_METHOD") || \ - (this_param == "DISP_PRED") || \ - (this_param == "AITKEN_RELAX") || \ + elif (this_param == "CFD_CONFIG_FILE_NAME") or \ + (this_param == "CSD_SOLVER") or \ + (this_param == "CSD_CONFIG_FILE_NAME") or \ + (this_param == "RESTART_SOL") or \ + (this_param == "MATCHING_MESH") or \ + (this_param == "MESH_INTERP_METHOD") or \ + (this_param == "DISP_PRED") or \ + (this_param == "AITKEN_RELAX") or \ (this_param == "TIME_MARCHING") : self._ConfigContent[this_param] = this_value diff --git a/SU2_PY/SU2_Nastran/pysu2_nastran.py b/SU2_PY/SU2_Nastran/pysu2_nastran.py index 53226db3bc1c..d7b431d71945 100644 --- a/SU2_PY/SU2_Nastran/pysu2_nastran.py +++ b/SU2_PY/SU2_Nastran/pysu2_nastran.py @@ -333,30 +333,30 @@ def __readConfig(self): this_value = line[1].strip() #integer values - if (this_param == "NMODES") || \ + if (this_param == "NMODES") or \ (this_param == "RESTART_ITER": self.Config[this_param] = int(this_value) #float values - elif (this_param == "DELTA_T") || \ - (this_param == "MODAL_DAMPING") || \ + elif (this_param == "DELTA_T") or \ + (this_param == "MODAL_DAMPING") or \ (this_param == "RHO"): self.Config[this_param] = float(this_value) #string values - elif (this_param == "TIME_MARCHING") || \ - (this_param == "MESH_FILE") || \ - (this_param == "PUNCH_FILE") || \ - (this_param == "RESTART_SOL") || \ + elif (this_param == "TIME_MARCHING") or \ + (this_param == "MESH_FILE") or \ + (this_param == "PUNCH_FILE") or \ + (this_param == "RESTART_SOL") or \ (this_param == "MOVING_MARKER"): self.Config[this_param] = this_value #lists values - elif (this_param == "INITIAL_MODES") || \ - (this_param == "IMPOSED_MODES") || \ + elif (this_param == "INITIAL_MODES") or \ + (this_param == "IMPOSED_MODES") or \ (this_param == "IMPOSED_PARAMETERS"): self.Config[this_param] = eval(this_value) From a285aed6765b7c477e591eb7184cb835e13be593 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 8 Feb 2021 21:07:00 +0100 Subject: [PATCH 05/32] Small issue with parenthesis --- SU2_PY/SU2_Nastran/pysu2_nastran.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_PY/SU2_Nastran/pysu2_nastran.py b/SU2_PY/SU2_Nastran/pysu2_nastran.py index d7b431d71945..00e3ba9a5a58 100644 --- a/SU2_PY/SU2_Nastran/pysu2_nastran.py +++ b/SU2_PY/SU2_Nastran/pysu2_nastran.py @@ -334,7 +334,7 @@ def __readConfig(self): #integer values if (this_param == "NMODES") or \ - (this_param == "RESTART_ITER": + (this_param == "RESTART_ITER"): self.Config[this_param] = int(this_value) From fbcfdc29c901a1a04f4c9eb9805b6ec27f2e9e72 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 8 Feb 2021 21:13:35 +0100 Subject: [PATCH 06/32] Uniformed indentation --- SU2_PY/SU2_Nastran/pysu2_nastran.py | 83 ++++++++++++++--------------- 1 file changed, 41 insertions(+), 42 deletions(-) diff --git a/SU2_PY/SU2_Nastran/pysu2_nastran.py b/SU2_PY/SU2_Nastran/pysu2_nastran.py index 00e3ba9a5a58..28a88c3f32eb 100644 --- a/SU2_PY/SU2_Nastran/pysu2_nastran.py +++ b/SU2_PY/SU2_Nastran/pysu2_nastran.py @@ -43,59 +43,58 @@ class ImposedMotionFunction: - def __init__(self,time0,tipo,parameters): - self.time0 = time0 - self.tipo = tipo - if self.tipo == "SINUSOIDAL": - self.bias = parameters[0] - self.amplitude = parameters[1] - self.frequency = parameters[2] - - elif self.tipo == "BLENDED_STEP": - self.kmax = parameters[0] - self.vinf = parameters[1] - self.lref = parameters[2] - self.amplitude = parameters[3] - self.tmax = 2*pi/self.kmax*self.lref/self.vinf - self.omega0 = 1/2*self.kmax - - else: - raise Exception('Imposed function {} not found, please implement it in pysu2_nastran.py'.format(self.tipo)) + def __init__(self,time0,tipo,parameters): + self.time0 = time0 + self.tipo = tipo + if self.tipo == "SINUSOIDAL": + self.bias = parameters[0] + self.amplitude = parameters[1] + self.frequency = parameters[2] + + elif self.tipo == "BLENDED_STEP": + self.kmax = parameters[0] + self.vinf = parameters[1] + self.lref = parameters[2] + self.amplitude = parameters[3] + self.tmax = 2*pi/self.kmax*self.lref/self.vinf + self.omega0 = 1/2*self.kmax + else: + raise Exception('Imposed function {} not found, please implement it in pysu2_nastran.py'.format(self.tipo)) - def GetDispl(self,time): - time = time - self.time0 - if self.tipo == "SINUSOIDAL": - return self.bias+self.amplitude*sin(2*pi*self.frequency*time) + def GetDispl(self,time): + time = time - self.time0 + if self.tipo == "SINUSOIDAL": + return self.bias+self.amplitude*sin(2*pi*self.frequency*time) - if self.tipo == "BLENDED_STEP": - if time < self.tmax: - return self.amplitude/2.0*(1.0-cos(self.omega0*time*self.vinf/self.lref)) - return self.amplitude + if self.tipo == "BLENDED_STEP": + if time < self.tmax: + return self.amplitude/2.0*(1.0-cos(self.omega0*time*self.vinf/self.lref)) + return self.amplitude - def GetVel(self,time): - time = time - self.time0 + def GetVel(self,time): + time = time - self.time0 - if self.tipo == "SINUSOIDAL": - return self.amplitude*cos(2*pi*self.frequency*time)*2*pi*self.frequency + if self.tipo == "SINUSOIDAL": + return self.amplitude*cos(2*pi*self.frequency*time)*2*pi*self.frequency - if self.tipo == "BLENDED_STEP": - if time < self.tmax: - return self.amplitude/2.0*sin(self.omega0*time*self.vinf/self.lref)*(self.omega0*self.vinf/self.lref) - return 0.0 + if self.tipo == "BLENDED_STEP": + if time < self.tmax: + return self.amplitude/2.0*sin(self.omega0*time*self.vinf/self.lref)*(self.omega0*self.vinf/self.lref) + return 0.0 - def GetAcc(self,time): - time = time - self.time0 + def GetAcc(self,time): + time = time - self.time0 - if self.tipo == "SINUSOIDAL": - return -self.amplitude*sin(2*pi*self.frequency*time)*(2*pi*self.frequency)**2 + if self.tipo == "SINUSOIDAL": + return -self.amplitude*sin(2*pi*self.frequency*time)*(2*pi*self.frequency)**2 - if self.tipo == "BLENDED_STEP": - if time < self.tmax: - return self.amplitude/2.0*cos(self.omega0*time*self.vinf/self.lref)*(self.omega0*self.vinf/self.lref)**2 - return 0.0 + if self.tipo == "BLENDED_STEP": + if time < self.tmax: + return self.amplitude/2.0*cos(self.omega0*time*self.vinf/self.lref)*(self.omega0*self.vinf/self.lref)**2 + return 0.0 class RefSystem: From 0dd42ad22c6bb5bec0c8f0dca6d27739071f82cf Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Fri, 12 Feb 2021 17:12:59 +0100 Subject: [PATCH 07/32] Update geo in the mesh deformation was not used --- SU2_CFD/src/drivers/CDriver.cpp | 2 +- SU2_CFD/src/solvers/CMeshSolver.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 7ee5c3571587..e68f59b2ae44 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1334,7 +1334,7 @@ void CDriver::Solver_Restart(CSolver ***solver, CGeometry **geometry, if ((restart || restart_flow) && config->GetDeform_Mesh() && update_geo){ /*--- Always restart with the last state ---*/ val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; - solver[MESH_0][MESH_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); + solver[MESH_0][MESH_SOL]->LoadRestart(geometry, solver, config, val_iter); } /*--- Exit if a restart was requested for a solver that is not available. ---*/ diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 5704f6c7bbd1..3874f65fd281 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -713,7 +713,7 @@ void CMeshSolver::SetDualTime_Mesh(void){ nodes->Set_Solution_time_n(); } -void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int val_iter, bool val_update_geo) { +void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int val_iter) { /*--- Read the restart data from either an ASCII or binary SU2 file. ---*/ From 2f6635a0b711a8c4493cfa0c4e0597daf8e257c8 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Fri, 12 Feb 2021 17:22:23 +0100 Subject: [PATCH 08/32] Updategeo removed from hpp files also --- SU2_CFD/include/solvers/CMeshSolver.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/SU2_CFD/include/solvers/CMeshSolver.hpp b/SU2_CFD/include/solvers/CMeshSolver.hpp index 55c14d7aa6f0..e1e597323b58 100644 --- a/SU2_CFD/include/solvers/CMeshSolver.hpp +++ b/SU2_CFD/include/solvers/CMeshSolver.hpp @@ -148,13 +148,11 @@ class CMeshSolver final : public CFEASolver { * \param[in] solver - Container vector with all of the solvers. * \param[in] config - Definition of the particular problem. * \param[in] val_iter - Current external iteration number. - * \param[in] val_update_geo - Flag for updating coords and grid velocity. */ void LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, - int val_iter, - bool val_update_geo) override; + int val_iter) override; /*! * \brief Load the geometries at the previous time states n and nM1. From 29b9d4c480cfeb52eaf62dbdb1ccf987c96f1305 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 13 Feb 2021 10:48:54 +0000 Subject: [PATCH 09/32] cleanup, try to make vertex tractions more general --- .../include/solvers/CFVMFlowSolverBase.inl | 19 +- SU2_CFD/include/solvers/CSolver.hpp | 10 +- .../src/iteration/CDiscAdjFluidIteration.cpp | 6 +- SU2_CFD/src/python_wrapper_structure.cpp | 241 +++++------------- SU2_CFD/src/solvers/CSolver.cpp | 168 +++++------- 5 files changed, 147 insertions(+), 297 deletions(-) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index f17b2bb21777..865b6b216a02 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -225,11 +225,20 @@ void CFVMFlowSolverBase::Allocate(const CConfig& config) { } } - /*--- Only initialize when there is a Marker_Fluid_Load defined - *--- (this avoids overhead in all other cases while a more permanent structure is being developed) ---*/ - if ((config.GetnMarker_Fluid_Load() > 0) && (MGLevel == MESH_0)) { - Alloc3D(nMarker, nVertex, nDim, VertexTraction); - if (config.GetDiscrete_Adjoint()) Alloc3D(nMarker, nVertex, nDim, VertexTractionAdjoint); + if (MGLevel == MESH_0) { + VertexTraction.resize(nMarker); + for (iMarker = 0; iMarker < nMarker; iMarker++) { + if (config.GetSolid_Wall(iMarker)) + VertexTraction[iMarker].resize(nVertex[iMarker], nDim) = su2double(0.0); + } + + if (config.GetDiscrete_Adjoint()) { + VertexTractionAdjoint.resize(nMarker); + for (iMarker = 0; iMarker < nMarker; iMarker++) { + if (config.GetSolid_Wall(iMarker)) + VertexTractionAdjoint[iMarker].resize(nVertex[iMarker], nDim) = su2double(0.0); + } + } } /*--- Initialize the BGS residuals in FSI problems. ---*/ diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 062711cb47b6..021ec21956e4 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -133,8 +133,8 @@ class CSolver { bool dynamic_grid; /*!< \brief Flag that determines whether the grid is dynamic (moving or deforming + grid velocities). */ - su2double ***VertexTraction; /*- Temporary, this will be moved to a new postprocessing structure once in place -*/ - su2double ***VertexTractionAdjoint; /*- Also temporary -*/ + vector VertexTraction; /*- Temporary, this will be moved to a new postprocessing structure once in place -*/ + vector VertexTractionAdjoint; /*- Also temporary -*/ string SolverName; /*!< \brief Store the name of the solver for output purposes. */ @@ -4339,7 +4339,7 @@ class CSolver { * \param[in] geometry - Geometrical definition. * \param[in] config - Definition of the particular problem. */ - void ComputeVertexTractions(CGeometry *geometry, CConfig *config); + void ComputeVertexTractions(CGeometry *geometry, const CConfig *config); /*! * \brief Set the adjoints of the vertex tractions. @@ -4356,7 +4356,7 @@ class CSolver { * \param[in] geometry - Geometrical definition. * \param[in] config - Definition of the particular problem. */ - void RegisterVertexTractions(CGeometry *geometry, CConfig *config); + void RegisterVertexTractions(CGeometry *geometry, const CConfig *config); /*! * \brief Store the adjoints of the vertex tractions. @@ -4377,7 +4377,7 @@ class CSolver { * \param[in] geometry - Geometrical definition. * \param[in] config - Definition of the particular problem. */ - void SetVertexTractionsAdjoint(CGeometry *geometry, CConfig *config); + void SetVertexTractionsAdjoint(CGeometry *geometry, const CConfig *config); /*! * \brief Get minimun volume in the mesh diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 9153084ae11c..94bfab7dc34f 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -370,7 +370,6 @@ void CDiscAdjFluidIteration::InitializeAdjoint(CSolver***** solver, CGeometry*** unsigned short iZone, unsigned short iInst) { bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); bool heat = config[iZone]->GetWeakly_Coupled_Heat(); - bool interface_boundary = (config[iZone]->GetnMarker_Fluid_Load() > 0); /*--- Initialize the adjoints the conservative variables ---*/ @@ -390,7 +389,7 @@ void CDiscAdjFluidIteration::InitializeAdjoint(CSolver***** solver, CGeometry*** solver[iZone][iInst][MESH_0][ADJRAD_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]); } - if (interface_boundary) { + if (config[iZone]->GetFluidProblem()) { solver[iZone][iInst][MESH_0][FLOW_SOL]->SetVertexTractionsAdjoint(geometry[iZone][iInst][MESH_0], config[iZone]); } } @@ -510,7 +509,6 @@ void CDiscAdjFluidIteration::RegisterOutput(CSolver***** solver, CGeometry**** g COutput* output, unsigned short iZone, unsigned short iInst) { bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); bool heat = config[iZone]->GetWeakly_Coupled_Heat(); - bool interface_boundary = (config[iZone]->GetnMarker_Fluid_Load() > 0); /*--- Register conservative variables as output of the iteration ---*/ @@ -526,7 +524,7 @@ void CDiscAdjFluidIteration::RegisterOutput(CSolver***** solver, CGeometry**** g if (config[iZone]->AddRadiation()) { solver[iZone][iInst][MESH_0][ADJRAD_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); } - if (interface_boundary) { + if (config[iZone]->GetFluidProblem()) { solver[iZone][iInst][MESH_0][FLOW_SOL]->RegisterVertexTractions(geometry[iZone][iInst][MESH_0], config[iZone]); } } diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index fde9093bc583..95cbf77ec6d7 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -79,19 +79,10 @@ passivedouble CDriver::Get_Drag() { unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); - su2double CDrag, RefDensity, RefArea, RefVel2, factor, val_Drag; - - /*--- Export free-stream density and reference area ---*/ - RefDensity = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetDensity_Inf(); - RefArea = config_container[val_iZone]->GetRefArea(); - - /*--- Calculate free-stream velocity (squared) ---*/ - RefVel2 = 0.0; - for(unsigned short iDim = 0; iDim < nDim; iDim++) - RefVel2 += pow(solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetVelocity_Inf(iDim),2); + su2double CDrag, factor, val_Drag; /*--- Calculate drag force based on drag coefficient ---*/ - factor = 0.5*RefDensity*RefArea*RefVel2; + factor = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetAeroCoeffsReferenceForce(); CDrag = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CD(); val_Drag = CDrag*factor; @@ -103,19 +94,10 @@ passivedouble CDriver::Get_Lift() { unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); - su2double CLift, RefDensity, RefArea, RefVel2, factor, val_Lift; - - /*--- Export free-stream density and reference area ---*/ - RefDensity = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetDensity_Inf(); - RefArea = config_container[val_iZone]->GetRefArea(); - - /*--- Calculate free-stream velocity (squared) ---*/ - RefVel2 = 0.0; - for(unsigned short iDim = 0; iDim < nDim; iDim++) - RefVel2 += pow(solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetVelocity_Inf(iDim),2); + su2double CLift, factor, val_Lift; /*--- Calculate drag force based on drag coefficient ---*/ - factor = 0.5*RefDensity*RefArea*RefVel2; + factor = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetAeroCoeffsReferenceForce(); CLift = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CL(); val_Lift = CLift*factor; @@ -127,100 +109,73 @@ passivedouble CDriver::Get_Mx(){ unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); - su2double CMx, RefDensity, RefArea, RefLengthCoeff, RefVel2, factor, val_Mx; + su2double CMx, RefLengthCoeff, factor, val_Mx; - /*--- Export free-stream density and reference area ---*/ - RefDensity = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetDensity_Inf(); - RefArea = config_container[val_iZone]->GetRefArea(); RefLengthCoeff = config_container[val_iZone]->GetRefLength(); - /*--- Calculate free-stream velocity (squared) ---*/ - RefVel2 = 0.0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) - RefVel2 += pow(solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetVelocity_Inf(iDim),2); - /*--- Calculate moment around x-axis based on coefficients ---*/ - factor = 0.5*RefDensity*RefArea*RefVel2; + factor = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetAeroCoeffsReferenceForce(); CMx = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CMx(); val_Mx = CMx*factor*RefLengthCoeff; return SU2_TYPE::GetValue(val_Mx); - } passivedouble CDriver::Get_My(){ unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); - su2double CMy, RefDensity, RefArea, RefLengthCoeff, RefVel2, factor, val_My; + su2double CMy, RefLengthCoeff, factor, val_My; - /*--- Export free-stream density and reference area ---*/ - RefDensity = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetDensity_Inf(); - RefArea = config_container[val_iZone]->GetRefArea(); RefLengthCoeff = config_container[val_iZone]->GetRefLength(); - /*--- Calculate free-stream velocity (squared) ---*/ - RefVel2 = 0.0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) - RefVel2 += pow(solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetVelocity_Inf(iDim),2); - /*--- Calculate moment around x-axis based on coefficients ---*/ - factor = 0.5*RefDensity*RefArea*RefVel2; + factor = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetAeroCoeffsReferenceForce(); CMy = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CMy(); val_My = CMy*factor*RefLengthCoeff; return SU2_TYPE::GetValue(val_My); - } passivedouble CDriver::Get_Mz() { unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); - su2double CMz, RefDensity, RefArea, RefLengthCoeff, RefVel2, factor, val_Mz; + su2double CMz, RefLengthCoeff, factor, val_Mz; - /*--- Export free-stream density and reference area ---*/ - RefDensity = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetDensity_Inf(); - RefArea = config_container[val_iZone]->GetRefArea(); RefLengthCoeff = config_container[val_iZone]->GetRefLength(); - /*--- Calculate free-stream velocity (squared) ---*/ - RefVel2 = 0.0; - for(unsigned short iDim = 0; iDim < nDim; iDim++) - RefVel2 += pow(solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetVelocity_Inf(iDim),2); - /*--- Calculate moment around z-axis based on coefficients ---*/ - factor = 0.5*RefDensity*RefArea*RefVel2; + factor = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetAeroCoeffsReferenceForce(); CMz = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CMz(); val_Mz = CMz*factor*RefLengthCoeff; return SU2_TYPE::GetValue(val_Mz); - } passivedouble CDriver::Get_DragCoeff() { - unsigned short val_iZone = ZONE_0; - unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); - su2double CDrag; + unsigned short val_iZone = ZONE_0; + unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); + su2double CDrag; - CDrag = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CD(); + CDrag = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CD(); - return SU2_TYPE::GetValue(CDrag); + return SU2_TYPE::GetValue(CDrag); } passivedouble CDriver::Get_LiftCoeff() { - unsigned short val_iZone = ZONE_0; - unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); - su2double CLift; + unsigned short val_iZone = ZONE_0; + unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); + su2double CLift; - CLift = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CL(); + CLift = solver_container[val_iZone][INST_0][FinestMesh][FLOW_SOL]->GetTotal_CL(); - return SU2_TYPE::GetValue(CLift); + return SU2_TYPE::GetValue(CLift); } unsigned short CDriver::GetMovingMarker() { @@ -305,55 +260,39 @@ passivedouble CDriver::GetUnsteady_TimeStep(){ passivedouble CDriver::GetVertexCoordX(unsigned short iMarker, unsigned long iVertex) { - su2double* Coord; - unsigned long iPoint; - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); + auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); + auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); return SU2_TYPE::GetValue(Coord[0]); - } passivedouble CDriver::GetVertexCoordY(unsigned short iMarker, unsigned long iVertex) { - su2double* Coord; - unsigned long iPoint; - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); + auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); + auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); return SU2_TYPE::GetValue(Coord[1]); } passivedouble CDriver::GetVertexCoordZ(unsigned short iMarker, unsigned long iVertex) { - su2double* Coord; - unsigned long iPoint; - - if(nDim == 3) { - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); - return SU2_TYPE::GetValue(Coord[2]); - } - else { - return 0.0; - } + if(nDim == 2) return 0.0; + auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); + auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); + return SU2_TYPE::GetValue(Coord[2]); } bool CDriver::ComputeVertexForces(unsigned short iMarker, unsigned long iVertex) { unsigned long iPoint; - unsigned short iDim, jDim; - su2double *Normal, AreaSquare, Area; - bool halo; + unsigned short iDim; + su2double *Normal, Area; unsigned short FinestMesh = config_container[ZONE_0]->GetFinestMesh(); /*--- Check the kind of fluid problem ---*/ bool compressible = (config_container[ZONE_0]->GetKind_Regime() == COMPRESSIBLE); bool incompressible = (config_container[ZONE_0]->GetKind_Regime() == INCOMPRESSIBLE); - bool viscous_flow = ((config_container[ZONE_0]->GetKind_Solver() == NAVIER_STOKES) || - (config_container[ZONE_0]->GetKind_Solver() == RANS) ); + bool viscous_flow = config_container[ZONE_0]->GetViscous(); /*--- Parameters for the calculations ---*/ // Pn: Pressure @@ -367,68 +306,50 @@ bool CDriver::ComputeVertexForces(unsigned short iMarker, unsigned long iVertex) iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); /*--- It is necessary to distinguish the halo nodes from the others, since they introduce non physical forces. ---*/ - if(geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetDomain(iPoint)) { - /*--- Get the normal at the vertex: this normal goes inside the fluid domain. ---*/ - Normal = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNormal(); - AreaSquare = 0.0; - for(iDim = 0; iDim < nDim; iDim++) { - AreaSquare += Normal[iDim]*Normal[iDim]; - } - Area = sqrt(AreaSquare); + if (!geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetDomain(iPoint)) + return true; - /*--- Get the values of pressure and viscosity ---*/ - Pn = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetPressure(iPoint); - if (viscous_flow) { - Viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); - } + /*--- Get the normal at the vertex: this normal goes inside the fluid domain. ---*/ + Normal = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNormal(); + Area = GeometryToolbox::Norm(nDim, Normal); - /*--- Calculate the inviscid (pressure) part of tn in the fluid nodes (force units) ---*/ - for (iDim = 0; iDim < nDim; iDim++) { - PyWrapNodalForce[iDim] = -(Pn-Pinf)*Normal[iDim]; //NB : norm(Normal) = Area - } + /*--- Get the values of pressure and viscosity ---*/ + Pn = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetPressure(iPoint); + if (viscous_flow) { + Viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); + } - /*--- Calculate the viscous (shear stress) part of tn in the fluid nodes (force units) ---*/ - if ((incompressible || compressible) && viscous_flow) { - CNumerics::ComputeStressTensor(nDim, Tau, - solver_container[ZONE_0][INST_0][FinestMesh][FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint)+1, Viscosity); - for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - PyWrapNodalForce[iDim] += Tau[iDim][jDim]*Normal[jDim]; - } - } - } + /*--- Calculate the inviscid (pressure) part of tn in the fluid nodes (force units) ---*/ + for (iDim = 0; iDim < nDim; iDim++) { + PyWrapNodalForce[iDim] = -(Pn-Pinf)*Normal[iDim]; //NB : norm(Normal) = Area + } - //Divide by local are in case of force density communication. - for(iDim = 0; iDim < nDim; iDim++) { - PyWrapNodalForceDensity[iDim] = PyWrapNodalForce[iDim]/Area; + /*--- Calculate the viscous (shear stress) part of tn in the fluid nodes (force units) ---*/ + if ((incompressible || compressible) && viscous_flow) { + CNumerics::ComputeStressTensor(nDim, Tau, + solver_container[ZONE_0][INST_0][FinestMesh][FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint)+1, Viscosity); + for (iDim = 0; iDim < nDim; iDim++) { + PyWrapNodalForce[iDim] += GeometryToolbox::DotProduct(nDim, Tau[iDim], Normal); } - - halo = false; - } - else { - halo = true; } - return halo; + //Divide by local are in case of force density communication. + for(iDim = 0; iDim < nDim; iDim++) + PyWrapNodalForceDensity[iDim] = PyWrapNodalForce[iDim]/Area; + return false; } passivedouble CDriver::GetVertexForceX(unsigned short iMarker, unsigned long iVertex) { - return SU2_TYPE::GetValue(PyWrapNodalForce[0]); - } passivedouble CDriver::GetVertexForceY(unsigned short iMarker, unsigned long iVertex) { - return SU2_TYPE::GetValue(PyWrapNodalForce[1]); - } passivedouble CDriver::GetVertexForceZ(unsigned short iMarker, unsigned long iVertex) { - return SU2_TYPE::GetValue(PyWrapNodalForce[2]); - } passivedouble CDriver::GetVertexForceDensityX(unsigned short iMarker, unsigned long iVertex) { @@ -445,34 +366,24 @@ passivedouble CDriver::GetVertexForceDensityZ(unsigned short iMarker, unsigned l void CDriver::SetVertexCoordX(unsigned short iMarker, unsigned long iVertex, passivedouble newPosX) { - unsigned long iPoint; - su2double *Coord; - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); + auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); + auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); PyWrapVarCoord[0] = newPosX - Coord[0]; - } void CDriver::SetVertexCoordY(unsigned short iMarker, unsigned long iVertex, passivedouble newPosY) { - unsigned long iPoint; - su2double *Coord; - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); + auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); + auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); PyWrapVarCoord[1] = newPosY - Coord[1]; } void CDriver::SetVertexCoordZ(unsigned short iMarker, unsigned long iVertex, passivedouble newPosZ) { - unsigned long iPoint; - su2double *Coord; - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); + auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); + auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); if(nDim > 2) { PyWrapVarCoord[2] = newPosZ - Coord[2]; @@ -486,8 +397,8 @@ passivedouble CDriver::SetVertexVarCoord(unsigned short iMarker, unsigned long i su2double nodalVarCoordNorm; - geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->SetVarCoord(PyWrapVarCoord); - nodalVarCoordNorm = sqrt((PyWrapVarCoord[0])*(PyWrapVarCoord[0]) + (PyWrapVarCoord[1])*(PyWrapVarCoord[1]) + (PyWrapVarCoord[2])*(PyWrapVarCoord[2])); + geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->SetVarCoord(PyWrapVarCoord); + nodalVarCoordNorm = sqrt((PyWrapVarCoord[0])*(PyWrapVarCoord[0]) + (PyWrapVarCoord[1])*(PyWrapVarCoord[1]) + (PyWrapVarCoord[2])*(PyWrapVarCoord[2])); return SU2_TYPE::GetValue(nodalVarCoordNorm); @@ -652,8 +563,6 @@ vector CDriver::GetVertexUnitNormal(unsigned short iMarker, unsig ret_Normal_passive[2] = SU2_TYPE::GetValue(ret_Normal[2]); return ret_Normal_passive; - - } vector CDriver::GetAllBoundaryMarkersTag(){ @@ -957,9 +866,7 @@ void CFluidDriver::BoundaryConditionsUpdate(){ int rank = MASTER_NODE; unsigned short iZone; -#ifdef HAVE_MPI - MPI_Comm_rank(SU2_MPI::GetComm(), &rank); -#endif + SU2_MPI::Comm_rank(SU2_MPI::GetComm(), &rank); if(rank == MASTER_NODE) cout << "Updating boundary conditions." << endl; for(iZone = 0; iZone < nZone; iZone++){ @@ -1070,11 +977,6 @@ vector CDriver::GetFEA_Velocity(unsigned short iMarker, unsigned else Velocity[2] = 0.0; } - else{ - Velocity[0] = 0.0; - Velocity[1] = 0.0; - Velocity[2] = 0.0; - } Velocity_passive[0] = SU2_TYPE::GetValue(Velocity[0]); Velocity_passive[1] = SU2_TYPE::GetValue(Velocity[1]); @@ -1101,11 +1003,6 @@ vector CDriver::GetFEA_Velocity_n(unsigned short iMarker, unsigne else Velocity_n[2] = 0.0; } - else{ - Velocity_n[0] = 0.0; - Velocity_n[1] = 0.0; - Velocity_n[2] = 0.0; - } Velocity_n_passive[0] = SU2_TYPE::GetValue(Velocity_n[0]); Velocity_n_passive[1] = SU2_TYPE::GetValue(Velocity_n[1]); @@ -1156,11 +1053,6 @@ vector CDriver::GetFlowLoad(unsigned short iMarker, unsigned long else FlowLoad[2] = 0.0; } - else{ - FlowLoad[0] = 0.0; - FlowLoad[1] = 0.0; - FlowLoad[2] = 0.0; - } FlowLoad_passive[0] = SU2_TYPE::GetValue(FlowLoad[0]); FlowLoad_passive[1] = SU2_TYPE::GetValue(FlowLoad[1]); @@ -1217,11 +1109,6 @@ vector CDriver::GetVertex_UndeformedCoord(unsigned short iMarker, else MeshCoord[2] = 0.0; } - else{ - MeshCoord[0] = 0.0; - MeshCoord[1] = 0.0; - MeshCoord[2] = 0.0; - } MeshCoord_passive[0] = SU2_TYPE::GetValue(MeshCoord[0]); MeshCoord_passive[1] = SU2_TYPE::GetValue(MeshCoord[1]); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index f608da291387..a29ab35b1fee 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -119,10 +119,6 @@ CSolver::CSolver(bool mesh_deform_mode) : System(mesh_deform_mode) { /*--- Flags for the dynamic grid (rigid movement or unsteady deformation). ---*/ dynamic_grid = false; - /*--- Container to store the vertex tractions. ---*/ - VertexTraction = nullptr; - VertexTractionAdjoint = nullptr; - /*--- Auxiliary data needed for CFL adaption. ---*/ Old_Func = 0; @@ -137,7 +133,6 @@ CSolver::CSolver(bool mesh_deform_mode) : System(mesh_deform_mode) { CSolver::~CSolver(void) { unsigned short iVar; - unsigned long iMarker, iVertex; /*--- Public variables, may be accessible outside ---*/ @@ -222,24 +217,6 @@ CSolver::~CSolver(void) { delete [] Jacobian_jj; } - if (VertexTraction != nullptr) { - for (iMarker = 0; iMarker < nMarker; iMarker++) { - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) - delete [] VertexTraction[iMarker][iVertex]; - delete [] VertexTraction[iMarker]; - } - delete [] VertexTraction; - } - - if (VertexTractionAdjoint != nullptr) { - for (iMarker = 0; iMarker < nMarker; iMarker++) { - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) - delete [] VertexTractionAdjoint[iMarker][iVertex]; - delete [] VertexTractionAdjoint[iMarker]; - } - delete [] VertexTractionAdjoint; - } - delete [] nVertex; delete [] Restart_Vars; @@ -4014,24 +3991,17 @@ void CSolver::LoadInletProfile(CGeometry **geometry, } -void CSolver::ComputeVertexTractions(CGeometry *geometry, CConfig *config){ +void CSolver::ComputeVertexTractions(CGeometry *geometry, const CConfig *config){ /*--- Compute the constant factor to dimensionalize pressure and shear stress. ---*/ - su2double *Velocity_ND, *Velocity_Real; + const su2double *Velocity_ND, *Velocity_Real; su2double Density_ND, Density_Real, Velocity2_Real, Velocity2_ND; su2double factor; - unsigned short iDim, jDim; + unsigned short iDim; // Check whether the problem is viscous - bool viscous_flow = ((config->GetKind_Solver() == NAVIER_STOKES) || - (config->GetKind_Solver() == INC_NAVIER_STOKES) || - (config->GetKind_Solver() == RANS) || - (config->GetKind_Solver() == INC_RANS) || - (config->GetKind_Solver() == DISC_ADJ_NAVIER_STOKES) || - (config->GetKind_Solver() == DISC_ADJ_INC_NAVIER_STOKES) || - (config->GetKind_Solver() == DISC_ADJ_INC_RANS) || - (config->GetKind_Solver() == DISC_ADJ_RANS)); + bool viscous_flow = config->GetViscous(); // Parameters for the calculations su2double Pn = 0.0; @@ -4039,7 +4009,7 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, CConfig *config){ unsigned short iMarker; unsigned long iVertex, iPoint; - su2double const *iNormal; + const su2double* iNormal; su2double Pressure_Inf = config->GetPressure_FreeStreamND(); @@ -4049,59 +4019,52 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, CConfig *config){ Velocity_ND = config->GetVelocity_FreeStreamND(); Density_ND = config->GetDensity_FreeStreamND(); - Velocity2_Real = 0.0; - Velocity2_ND = 0.0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) { - Velocity2_Real += Velocity_Real[iDim]*Velocity_Real[iDim]; - Velocity2_ND += Velocity_ND[iDim]*Velocity_ND[iDim]; - } + Velocity2_Real = GeometryToolbox::SquaredNorm(nDim, Velocity_Real); + Velocity2_ND = GeometryToolbox::SquaredNorm(nDim, Velocity_ND); factor = Density_Real * Velocity2_Real / ( Density_ND * Velocity2_ND ); for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - /*--- If this is defined as an interface marker ---*/ - if (config->GetMarker_All_Fluid_Load(iMarker) == YES) { + /*--- If this is defined as a wall ---*/ + if (!config->GetSolid_Wall(iMarker)) continue; - // Loop over the vertices - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + // Loop over the vertices + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - // Recover the point index - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - // Get the normal at the vertex: this normal goes inside the fluid domain. - iNormal = geometry->vertex[iMarker][iVertex]->GetNormal(); + // Recover the point index + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + // Get the normal at the vertex: this normal goes inside the fluid domain. + iNormal = geometry->vertex[iMarker][iVertex]->GetNormal(); - /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { + /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ + if (geometry->nodes->GetDomain(iPoint)) { - // Retrieve the values of pressure - Pn = base_nodes->GetPressure(iPoint); + // Retrieve the values of pressure + Pn = base_nodes->GetPressure(iPoint); - // Calculate tn in the fluid nodes for the inviscid term --> Units of force (non-dimensional). - for (iDim = 0; iDim < nDim; iDim++) - auxForce[iDim] = -(Pn-Pressure_Inf)*iNormal[iDim]; + // Calculate tn in the fluid nodes for the inviscid term --> Units of force (non-dimensional). + for (iDim = 0; iDim < nDim; iDim++) + auxForce[iDim] = -(Pn-Pressure_Inf)*iNormal[iDim]; - // Calculate tn in the fluid nodes for the viscous term - if (viscous_flow) { - su2double Viscosity = base_nodes->GetLaminarViscosity(iPoint); - su2double Tau[3][3]; - CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetGradient_Primitive(iPoint)+1, Viscosity); - for (iDim = 0; iDim < nDim; iDim++) { - for (jDim = 0 ; jDim < nDim; jDim++) { - auxForce[iDim] += Tau[iDim][jDim]*iNormal[jDim]; - } - } - } - - // Redimensionalize the forces + // Calculate tn in the fluid nodes for the viscous term + if (viscous_flow) { + su2double Viscosity = base_nodes->GetLaminarViscosity(iPoint); + su2double Tau[3][3]; + CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetGradient_Primitive(iPoint)+1, Viscosity); for (iDim = 0; iDim < nDim; iDim++) { - VertexTraction[iMarker][iVertex][iDim] = factor * auxForce[iDim]; + auxForce[iDim] += GeometryToolbox::DotProduct(nDim, Tau[iDim], iNormal); } } - else{ - for (iDim = 0; iDim < nDim; iDim++) { - VertexTraction[iMarker][iVertex][iDim] = 0.0; - } + + // Redimensionalize the forces + for (iDim = 0; iDim < nDim; iDim++) { + VertexTraction[iMarker][iVertex][iDim] = factor * auxForce[iDim]; + } + } + else{ + for (iDim = 0; iDim < nDim; iDim++) { + VertexTraction[iMarker][iVertex][iDim] = 0.0; } } } @@ -4109,7 +4072,7 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, CConfig *config){ } -void CSolver::RegisterVertexTractions(CGeometry *geometry, CConfig *config){ +void CSolver::RegisterVertexTractions(CGeometry *geometry, const CConfig *config){ unsigned short iMarker, iDim; unsigned long iVertex, iPoint; @@ -4117,31 +4080,28 @@ void CSolver::RegisterVertexTractions(CGeometry *geometry, CConfig *config){ /*--- Loop over all the markers ---*/ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - /*--- If this is defined as an interface marker ---*/ - if (config->GetMarker_All_Fluid_Load(iMarker) == YES) { + /*--- If this is defined as a wall ---*/ + if (!config->GetSolid_Wall(iMarker)) continue; - /*--- Loop over the vertices ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { + /*--- Loop over the vertices ---*/ + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - /*--- Recover the point index ---*/ - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + /*--- Recover the point index ---*/ + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Register the vertex traction as output ---*/ - for (iDim = 0; iDim < nDim; iDim++) { - AD::RegisterOutput(VertexTraction[iMarker][iVertex][iDim]); - } + /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ + if (!geometry->nodes->GetDomain(iPoint)) continue; - } + /*--- Register the vertex traction as output ---*/ + for (iDim = 0; iDim < nDim; iDim++) { + AD::RegisterOutput(VertexTraction[iMarker][iVertex][iDim]); } } } } -void CSolver::SetVertexTractionsAdjoint(CGeometry *geometry, CConfig *config){ +void CSolver::SetVertexTractionsAdjoint(CGeometry *geometry, const CConfig *config){ unsigned short iMarker, iDim; unsigned long iVertex, iPoint; @@ -4149,26 +4109,22 @@ void CSolver::SetVertexTractionsAdjoint(CGeometry *geometry, CConfig *config){ /*--- Loop over all the markers ---*/ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - /*--- If this is defined as an interface marker ---*/ - if (config->GetMarker_All_Fluid_Load(iMarker) == YES) { - - /*--- Loop over the vertices ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - - /*--- Recover the point index ---*/ - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + /*--- If this is defined as a wall ---*/ + if (!config->GetSolid_Wall(iMarker)) continue; - /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { + /*--- Loop over the vertices ---*/ + for (iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++) { - /*--- Set the adjoint of the vertex traction from the value received ---*/ - for (iDim = 0; iDim < nDim; iDim++) { + /*--- Recover the point index ---*/ + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - SU2_TYPE::SetDerivative(VertexTraction[iMarker][iVertex][iDim], - SU2_TYPE::GetValue(VertexTractionAdjoint[iMarker][iVertex][iDim])); - } + /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ + if (!geometry->nodes->GetDomain(iPoint)) continue; - } + /*--- Set the adjoint of the vertex traction from the value received ---*/ + for (iDim = 0; iDim < nDim; iDim++) { + SU2_TYPE::SetDerivative(VertexTraction[iMarker][iVertex][iDim], + SU2_TYPE::GetValue(VertexTractionAdjoint[iMarker][iVertex][iDim])); } } } From 0f9342fd6d81fed5557e19a806f1223f2b6dbae9 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Sat, 13 Feb 2021 15:13:50 +0100 Subject: [PATCH 10/32] Several modifications to clean up the FSI interface --- SU2_PY/FSI_tools/FSIInterface.py | 280 ++++++++++++++-------------- SU2_PY/SU2_Nastran/pysu2_nastran.py | 6 +- 2 files changed, 146 insertions(+), 140 deletions(-) diff --git a/SU2_PY/FSI_tools/FSIInterface.py b/SU2_PY/FSI_tools/FSIInterface.py index d3d10cc58191..62438de8388b 100644 --- a/SU2_PY/FSI_tools/FSIInterface.py +++ b/SU2_PY/FSI_tools/FSIInterface.py @@ -114,7 +114,9 @@ def __init__(self, FSI_config, FluidSolver, SolidSolver, have_MPI): self.localFluidInterface_array_Y_init = None self.localFluidInterface_array_Z_init = None - self.haloNodesPositionsInit = {} #initial position of the halo nodes (fluid side only) + self.localSolidInterface_array_X_init = None #initial solid interface position on each partition (used for mesh mapping) + self.localSolidInterface_array_Y_init = None + self.localSolidInterface_array_Z_init = None self.solidInterface_array_DispX = None #solid interface displacement self.solidInterface_array_DispY = None @@ -311,6 +313,7 @@ def connect(self, FSI_config, FluidSolver, SolidSolver): # Same thing for the solid part self.nLocalSolidInterfaceHaloNode = 0 + # TODO when the solid solver will run in parallel, add here the calculation of halo nodes self.nLocalSolidInterfacePhysicalNodes = self.nLocalSolidInterfaceNodes - self.nLocalSolidInterfaceHaloNode if self.have_MPI: self.SolidHaloNodeList = self.comm.allgather(self.SolidHaloNodeList) @@ -319,36 +322,36 @@ def connect(self, FSI_config, FluidSolver, SolidSolver): # --- Calculate the total number of nodes (with and without halo) at the fluid interface (sum over all the partitions) and broadcast the number accross all processors --- - sendBuffHalo = np.array(int(self.nLocalFluidInterfaceNodes)) + sendBuffTotal = np.array(int(self.nLocalFluidInterfaceNodes)) sendBuffPhysical = np.array(int(self.nLocalFluidInterfacePhysicalNodes)) - rcvBuffHalo = np.zeros(1, dtype=int) + rcvBuffTotal = np.zeros(1, dtype=int) rcvBuffPhysical = np.zeros(1, dtype=int) if self.have_MPI: self.comm.barrier() - self.comm.Allreduce(sendBuffHalo,rcvBuffHalo,op=self.MPI.SUM) + self.comm.Allreduce(sendBuffTotal,rcvBuffTotal,op=self.MPI.SUM) self.comm.Allreduce(sendBuffPhysical,rcvBuffPhysical,op=self.MPI.SUM) - self.nFluidInterfaceNodes = rcvBuffHalo[0] + self.nFluidInterfaceNodes = rcvBuffTotal[0] self.nFluidInterfacePhysicalNodes = rcvBuffPhysical[0] else: - self.nFluidInterfaceNodes = np.copy(sendBuffHalo) + self.nFluidInterfaceNodes = np.copy(sendBuffTotal) self.nFluidInterfacePhysicalNodes = np.copy(sendBuffPhysical) - del sendBuffHalo, rcvBuffHalo, sendBuffPhysical, rcvBuffPhysical + del sendBuffTotal, rcvBuffTotal, sendBuffPhysical, rcvBuffPhysical # Same thing for the solid part - sendBuffHalo = np.array(int(self.nLocalSolidInterfaceNodes)) + sendBuffTotal = np.array(int(self.nLocalSolidInterfaceNodes)) sendBuffPhysical = np.array(int(self.nLocalSolidInterfacePhysicalNodes)) - rcvBuffHalo = np.zeros(1, dtype=int) + rcvBuffTotal = np.zeros(1, dtype=int) rcvBuffPhysical = np.zeros(1, dtype=int) if self.have_MPI: self.comm.barrier() - self.comm.Allreduce(sendBuffHalo,rcvBuffHalo,op=self.MPI.SUM) + self.comm.Allreduce(sendBuffTotal,rcvBuffTotal,op=self.MPI.SUM) self.comm.Allreduce(sendBuffPhysical,rcvBuffPhysical,op=self.MPI.SUM) - self.nSolidInterfaceNodes = rcvBuffHalo[0] + self.nSolidInterfaceNodes = rcvBuffTotal[0] self.nSolidInterfacePhysicalNodes = rcvBuffPhysical[0] else: - self.nSolidInterfaceNodes = np.copy(sendBuffHalo) + self.nSolidInterfaceNodes = np.copy(sendBuffTotal) self.nSolidInterfacePhysicalNodes = np.copy(sendBuffPhysical) - del sendBuffHalo, rcvBuffHalo, sendBuffPhysical, rcvBuffPhysical + del sendBuffTotal, rcvBuffTotal, sendBuffPhysical, rcvBuffPhysical # --- Store the number of physical interface nodes on each processor and allgather the information --- self.fluidPhysicalInterfaceNodesDistribution = np.zeros(MPIsize, dtype=int) @@ -362,7 +365,7 @@ def connect(self, FSI_config, FluidSolver, SolidSolver): # Same thing for the solid part self.solidPhysicalInterfaceNodesDistribution = np.zeros(MPIsize, dtype=int) if self.have_MPI: - sendBuffPhysical = np.array(int(self.nLocalSolidInterfaceNodes)) + sendBuffPhysical = np.array(int(self.nLocalSolidInterfacePhysicalNodes)) self.comm.Allgather(sendBuffPhysical,self.solidPhysicalInterfaceNodesDistribution) del sendBuffPhysical else: @@ -392,7 +395,7 @@ def connect(self, FSI_config, FluidSolver, SolidSolver): globalIndexStart = 0 for iProc in range(myid): globalIndexStart += self.solidPhysicalInterfaceNodesDistribution[iProc] - globalIndexStop = globalIndexStart + self.nLocalSolidInterfaceNodes-1 + globalIndexStop = globalIndexStart + self.nLocalSolidInterfacePhysicalNodes-1 else: globalIndexStart = 0 globalIndexStop = 0 @@ -562,13 +565,11 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config): # Note that the fluid solver is separated in more processors outside the python script # thus when, from a core, we request for the vertices on the interface, we only obtain # those in that node - GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex) + GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex) #TODO obtain here the undeformed mesh posx = FluidSolver.GetVertexCoordX(self.fluidInterfaceIdentifier, iVertex) posy = FluidSolver.GetVertexCoordY(self.fluidInterfaceIdentifier, iVertex) posz = FluidSolver.GetVertexCoordZ(self.fluidInterfaceIdentifier, iVertex) - if GlobalIndex in self.FluidHaloNodeList[myid].keys(): - self.haloNodesPositionsInit[GlobalIndex] = (posx, posy, posz) - else: + if GlobalIndex not in self.FluidHaloNodeList[myid].keys(): fluidIndexing_temp[GlobalIndex] = self.__getGlobalIndex('fluid', myid, localIndex) self.localFluidInterface_array_X_init[localIndex] = posx self.localFluidInterface_array_Y_init[localIndex] = posy @@ -587,19 +588,17 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config): # --- Get the solid interface from solid solver on each partition --- localIndex = 0 solidIndexing_temp = {} - self.localSolidInterface_array_X = np.zeros(self.nLocalSolidInterfaceNodes) - self.localSolidInterface_array_Y = np.zeros(self.nLocalSolidInterfaceNodes) - self.localSolidInterface_array_Z = np.zeros(self.nLocalSolidInterfaceNodes) + self.localSolidInterface_array_X_init = np.zeros(self.nLocalSolidInterfaceNodes) + self.localSolidInterface_array_Y_init = np.zeros(self.nLocalSolidInterfaceNodes) + self.localSolidInterface_array_Z_init = np.zeros(self.nLocalSolidInterfaceNodes) for iVertex in range(self.nLocalSolidInterfaceNodes): GlobalIndex = SolidSolver.getInterfaceNodeGlobalIndex(self.solidInterfaceIdentifier, iVertex) - posx, posy, posz = SolidSolver.getInterfaceNodePos(self.solidInterfaceIdentifier, iVertex) - if GlobalIndex in self.SolidHaloNodeList[myid].keys(): - pass - else: + posx, posy, posz = SolidSolver.getInterfaceNodePosInit(self.solidInterfaceIdentifier, iVertex) + if GlobalIndex not in self.SolidHaloNodeList[myid].keys(): solidIndexing_temp[GlobalIndex] = self.__getGlobalIndex('solid', myid, localIndex) - self.localSolidInterface_array_X[localIndex] = posx - self.localSolidInterface_array_Y[localIndex] = posy - self.localSolidInterface_array_Z[localIndex] = posz + self.localSolidInterface_array_X_init[localIndex] = posx + self.localSolidInterface_array_Y_init[localIndex] = posy + self.localSolidInterface_array_Z_init[localIndex] = posz localIndex += 1 if self.have_MPI: solidIndexing_temp = self.comm.allgather(solidIndexing_temp) @@ -682,17 +681,17 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config): if myid == iProc: for jProc in self.solidInterfaceProcessors: if jProc != iProc: - self.comm.ssend(self.localSolidInterface_array_X, dest=jProc, tag=1) - self.comm.ssend(self.localSolidInterface_array_Y, dest=jProc, tag=2) - self.comm.ssend(self.localSolidInterface_array_Z, dest=jProc, tag=3) + self.comm.ssend(self.localSolidInterface_array_X_init, dest=jProc, tag=1) + self.comm.ssend(self.localSolidInterface_array_Y_init, dest=jProc, tag=2) + self.comm.ssend(self.localSolidInterface_array_Z_init, dest=jProc, tag=3) else: sizeOfBuff = self.solidPhysicalInterfaceNodesDistribution[iProc] solidInterfaceBuffRcv_X = np.zeros(sizeOfBuff) solidInterfaceBuffRcv_Y = np.zeros(sizeOfBuff) solidInterfaceBuffRcv_Z = np.zeros(sizeOfBuff) - solidInterfaceBuffRcv_X = self.localSolidInterface_array_X - solidInterfaceBuffRcv_Y = self.localSolidInterface_array_Y - solidInterfaceBuffRcv_Z = self.localSolidInterface_array_Z + solidInterfaceBuffRcv_X = self.localSolidInterface_array_X_init + solidInterfaceBuffRcv_Y = self.localSolidInterface_array_Y_init + solidInterfaceBuffRcv_Z = self.localSolidInterface_array_Z_init if myid in self.solidInterfaceProcessors: if myid != iProc: sizeOfBuff = self.solidPhysicalInterfaceNodesDistribution[iProc] @@ -708,9 +707,9 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config): self.TPSMeshMapping_A(solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc) else: if FSI_config['MESH_INTERP_METHOD'] == 'RBF': - self.RBFMeshMapping_A(self.localSolidInterface_array_X, self.localSolidInterface_array_Y, self.localSolidInterface_array_Z, 0, self.RBF_rad) + self.RBFMeshMapping_A(self.localSolidInterface_array_X_init, self.localSolidInterface_array_Y_init, self.localSolidInterface_array_Z_init, 0, self.RBF_rad) else: - self.TPSMeshMapping_A(self.localSolidInterface_array_X, self.localSolidInterface_array_Y, self.localSolidInterface_array_Z, 0) + self.TPSMeshMapping_A(self.localSolidInterface_array_X_init, self.localSolidInterface_array_Y_init, self.localSolidInterface_array_Z_init, 0) self.MappingMatrixA.assemblyBegin() self.MappingMatrixA.assemblyEnd() self.MappingMatrixA_T.assemblyBegin() @@ -724,17 +723,17 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config): if myid == iProc: for jProc in self.fluidInterfaceProcessors: if jProc != iProc: - self.comm.ssend(self.localSolidInterface_array_X, dest=jProc, tag=1) - self.comm.ssend(self.localSolidInterface_array_Y, dest=jProc, tag=2) - self.comm.ssend(self.localSolidInterface_array_Z, dest=jProc, tag=3) + self.comm.ssend(self.localSolidInterface_array_X_init, dest=jProc, tag=1) + self.comm.ssend(self.localSolidInterface_array_Y_init, dest=jProc, tag=2) + self.comm.ssend(self.localSolidInterface_array_Z_init, dest=jProc, tag=3) else: sizeOfBuff = self.solidPhysicalInterfaceNodesDistribution[iProc] solidInterfaceBuffRcv_X = np.zeros(sizeOfBuff) solidInterfaceBuffRcv_Y = np.zeros(sizeOfBuff) solidInterfaceBuffRcv_Z = np.zeros(sizeOfBuff) - solidInterfaceBuffRcv_X = self.localSolidInterface_array_X - solidInterfaceBuffRcv_Y = self.localSolidInterface_array_Y - solidInterfaceBuffRcv_Z = self.localSolidInterface_array_Z + solidInterfaceBuffRcv_X = self.localSolidInterface_array_X_init + solidInterfaceBuffRcv_Y = self.localSolidInterface_array_Y_init + solidInterfaceBuffRcv_Z = self.localSolidInterface_array_Z_init if myid in self.fluidInterfaceProcessors: if myid != iProc: sizeOfBuff = self.solidPhysicalInterfaceNodesDistribution[iProc] @@ -756,13 +755,13 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config): else: if FSI_config['MATCHING_MESH'] == 'NO': if FSI_config['MESH_INTERP_METHOD'] == 'RBF': - self.RBFMeshMapping_B(self.localSolidInterface_array_X, self.localSolidInterface_array_Y, self.localSolidInterface_array_Z, 0, self.RBF_rad) + self.RBFMeshMapping_B(self.localSolidInterface_array_X_init, self.localSolidInterface_array_Y_init, self.localSolidInterface_array_Z_init, 0, self.RBF_rad) elif FSI_config['MESH_INTERP_METHOD'] == 'TPS' : - self.TPSMeshMapping_B(self.localSolidInterface_array_X, self.localSolidInterface_array_Y, self.localSolidInterface_array_Z, 0) + self.TPSMeshMapping_B(self.localSolidInterface_array_X_init, self.localSolidInterface_array_Y_init, self.localSolidInterface_array_Z_init, 0) else: - self.NearestNeighboorMeshMapping(self.localSolidInterface_array_X, self.localSolidInterface_array_Y, self.localSolidInterface_array_Z, 0) + self.NearestNeighboorMeshMapping(self.localSolidInterface_array_X_init, self.localSolidInterface_array_Y_init, self.localSolidInterface_array_Z_init, 0) else: - self.matchingMeshMapping(self.localSolidInterface_array_X, self.localSolidInterface_array_Y, self.localSolidInterface_array_Z, 0) + self.matchingMeshMapping(self.localSolidInterface_array_X_init, self.localSolidInterface_array_Y_init, self.localSolidInterface_array_Z_init, 0) if FSI_config['MATCHING_MESH'] == 'NO' and (FSI_config['MESH_INTERP_METHOD'] == 'RBF' or FSI_config['MESH_INTERP_METHOD'] == 'TPS'): self.MappingMatrixB.assemblyBegin() @@ -779,9 +778,9 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config): self.MPIBarrier() - del self.localSolidInterface_array_X - del self.localSolidInterface_array_Y - del self.localSolidInterface_array_Z + del self.localSolidInterface_array_X_init + del self.localSolidInterface_array_Y_init + del self.localSolidInterface_array_Z_init del self.localFluidInterface_array_X_init del self.localFluidInterface_array_Y_init del self.localFluidInterface_array_Z_init @@ -906,10 +905,10 @@ def RBFMeshMapping_A(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, sol else : SolidSpatialTree.add(jVertex, (posX, posY, posZ)) - for iVertexSolid in range(self.nLocalSolidInterfaceNodes): - posX = self.localSolidInterface_array_X[iVertexSolid] - posY = self.localSolidInterface_array_Y[iVertexSolid] - posZ = self.localSolidInterface_array_Z[iVertexSolid] + for iVertexSolid in range(self.nLocalSolidInterfacePhysicalNodes): + posX = self.localSolidInterface_array_X_init[iVertexSolid] + posY = self.localSolidInterface_array_Y_init[iVertexSolid] + posZ = self.localSolidInterface_array_Z_init[iVertexSolid] NodeA = np.array([posX, posY, posZ]) iGlobalVertexSolid = self.__getGlobalIndex('solid', myid, iVertexSolid) if self.nDim == 2: @@ -1003,10 +1002,10 @@ def TPSMeshMapping_A(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, sol nSolidNodes = solidInterfaceBuffRcv_X.shape[0] - for iVertexSolid in range(self.nLocalSolidInterfaceNodes): - posX = self.localSolidInterface_array_X[iVertexSolid] - posY = self.localSolidInterface_array_Y[iVertexSolid] - posZ = self.localSolidInterface_array_Z[iVertexSolid] + for iVertexSolid in range(self.nLocalSolidInterfacePhysicalNodes): + posX = self.localSolidInterface_array_X_init[iVertexSolid] + posY = self.localSolidInterface_array_Y_init[iVertexSolid] + posZ = self.localSolidInterface_array_Z_init[iVertexSolid] NodeA = np.array([posX, posY, posZ]) iGlobalVertexSolid = self.__getGlobalIndex('solid', myid, iVertexSolid) for jVertexSolid in range(nSolidNodes): @@ -1135,9 +1134,7 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config): KSP_solver.getPC().setType('jacobi') KSP_solver.setOperators(self.MappingMatrixA) KSP_solver.setFromOptions() - #print(KSP_solver.getInitialGuessNonzero()) KSP_solver.setInitialGuessNonzero(True) - #print(KSP_solver.getInitialGuessNonzero()) KSP_solver.solve(self.solidInterface_array_DispX, gamma_array_DispX) KSP_solver.solve(self.solidInterface_array_DispY, gamma_array_DispY) if self.nDim==3: @@ -1204,9 +1201,6 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config): del sendBuffNumber, rcvBuffNumber - #print("DEBUG MESSAGE From proc {}, counts = {}".format(myid, counts)) - #print("DEBUG MESSAGE From proc {}, displ = {}".format(myid, displ)) - self.comm.Gatherv(self.fluidInterface_array_DispX.getArray(), [self.fluidInterface_array_DispX_recon, counts, displ, self.MPI.DOUBLE], root=self.rootProcess) self.comm.Gatherv(self.fluidInterface_array_DispY.getArray(), [self.fluidInterface_array_DispY_recon, counts, displ, self.MPI.DOUBLE], root=self.rootProcess) self.comm.Gatherv(self.fluidInterface_array_DispZ.getArray(), [self.fluidInterface_array_DispZ_recon, counts, displ, self.MPI.DOUBLE], root=self.rootProcess) @@ -1227,7 +1221,7 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config): self.localFluidInterface_array_DispX = np.copy(sendBuff_X) self.localFluidInterface_array_DispY = np.copy(sendBuff_Y) self.localFluidInterface_array_DispZ = np.copy(sendBuff_Z) - if iProc != self.rootProcess: + else: self.comm.ssend(sendBuff_X, dest=iProc, tag = 1) self.comm.ssend(sendBuff_Y, dest=iProc, tag = 2) self.comm.ssend(sendBuff_Z, dest=iProc, tag = 3) @@ -1383,17 +1377,17 @@ def interpolateFluidLoadsOnSolidMesh(self, FSI_config): self.comm.ssend(sendBuff_Y, dest=iProc, tag = 2) self.comm.ssend(sendBuff_Z, dest=iProc, tag = 3) else: - self.localSolidLoads_array_X = np.zeros(self.nLocalSolidInterfaceNodes) - self.localSolidLoads_array_Y = np.zeros(self.nLocalSolidInterfaceNodes) - self.localSolidLoads_array_Z = np.zeros(self.nLocalSolidInterfaceNodes) + self.localSolidLoads_array_X = np.zeros(self.nLocalSolidInterfacePhysicalNodes) + self.localSolidLoads_array_Y = np.zeros(self.nLocalSolidInterfacePhysicalNodes) + self.localSolidLoads_array_Z = np.zeros(self.nLocalSolidInterfacePhysicalNodes) self.localSolidLoads_array_X = sendBuff_X self.localSolidLoads_array_Y = sendBuff_Y self.localSolidLoads_array_Z = sendBuff_Z if myid in self.solidInterfaceProcessors: if myid != self.rootProcess: - self.localSolidLoads_array_X = np.zeros(self.nLocalSolidInterfaceNodes) - self.localSolidLoads_array_Y = np.zeros(self.nLocalSolidInterfaceNodes) - self.localSolidLoads_array_Z = np.zeros(self.nLocalSolidInterfaceNodes) + self.localSolidLoads_array_X = np.zeros(self.nLocalSolidInterfacePhysicalNodes) + self.localSolidLoads_array_Y = np.zeros(self.nLocalSolidInterfacePhysicalNodes) + self.localSolidLoads_array_Z = np.zeros(self.nLocalSolidInterfacePhysicalNodes) self.localSolidLoads_array_X = self.comm.recv(source=self.rootProcess, tag = 1) self.localSolidLoads_array_Y = self.comm.recv(source=self.rootProcess, tag = 2) self.localSolidLoads_array_Z = self.comm.recv(source=self.rootProcess, tag = 3) @@ -1408,8 +1402,31 @@ def interpolateFluidLoadsOnSolidMesh(self, FSI_config): self.localSolidLoads_array_Y = self.solidLoads_array_Y.getArray().copy() self.localSolidLoads_array_Z = self.solidLoads_array_Z.getArray().copy() - # Special treatment for the halo nodes on the fluid interface - # TODO when we will use parallel solid solver !! + # Special treatment for the halo nodes on the solid interface + self.haloNodesLoads = {} + sendBuff = {} + if self.have_MPI: + if myid == self.rootProcess: + for iProc in self.solidInterfaceProcessors: + sendBuff = {} + for key in self.SolidHaloNodeList[iProc].keys(): + globalIndex = self.solidIndexing[key] + DispX = self.solidLoads_array_X_recon[globalIndex] + DispY = self.solidLoads_array_Y_recon[globalIndex] + DispZ = self.solidLoads_array_Z_recon[globalIndex] + sendBuff[key] = (DispX, DispY, DispZ) + if iProc == self.rootProcess: + self.haloNodesLoads = sendBuff + else: + self.comm.ssend(sendBuff, dest = iProc, tag=4) + if myid in self.solidInterfaceProcessors: + if myid != self.rootProcess: + self.haloNodesLoads = self.comm.recv(source = self.rootProcess, tag = 4) + self.comm.barrier() + del self.solidLoads_array_X_recon + del self.solidLoads_array_Y_recon + del self.solidLoads_array_Z_recon + del sendBuff def getSolidInterfaceDisplacement(self, SolidSolver): """ @@ -1425,9 +1442,7 @@ def getSolidInterfaceDisplacement(self, SolidSolver): localIndex = 0 for iVertex in range(self.nLocalSolidInterfaceNodes): GlobalIndex = SolidSolver.getInterfaceNodeGlobalIndex(self.solidInterfaceIdentifier, iVertex) - if GlobalIndex in self.SolidHaloNodeList[myid].keys(): - pass - else: + if GlobalIndex not in self.SolidHaloNodeList[myid].keys(): newDispx, newDispy, newDispz = SolidSolver.getInterfaceNodeDisp(self.solidInterfaceIdentifier, iVertex) iGlobalVertex = self.__getGlobalIndex('solid', myid, localIndex) self.solidInterface_array_DispX.setValues([iGlobalVertex],newDispx) @@ -1451,10 +1466,8 @@ def getFluidInterfaceNodalForce(self, FSI_config, FluidSolver): else: myid = 0 + GlobalIndex = int() localIndex = 0 - FX = 0.0 - FY = 0.0 - FZ = 0.0 # --- Get the fluid interface loads from the fluid solver and directly fill the corresponding PETSc vector --- for iVertex in range(self.nLocalFluidInterfaceNodes): @@ -1465,16 +1478,8 @@ def getFluidInterfaceNodalForce(self, FSI_config, FluidSolver): self.fluidLoads_array_X.setValues([iGlobalVertex], loadX) self.fluidLoads_array_Y.setValues([iGlobalVertex], loadY) self.fluidLoads_array_Z.setValues([iGlobalVertex], loadZ) - FX += loadX - FY += loadY - FZ += loadZ localIndex += 1 - if self.have_MPI: - FX = self.comm.allreduce(FX) - FY = self.comm.allreduce(FY) - FZ = self.comm.allreduce(FZ) - self.fluidLoads_array_X.assemblyBegin() self.fluidLoads_array_X.assemblyEnd() self.fluidLoads_array_Y.assemblyBegin() @@ -1482,10 +1487,6 @@ def getFluidInterfaceNodalForce(self, FSI_config, FluidSolver): self.fluidLoads_array_Z.assemblyBegin() self.fluidLoads_array_Z.assemblyEnd() - FX_b = self.fluidLoads_array_X.sum() - FY_b = self.fluidLoads_array_Y.sum() - FZ_b = self.fluidLoads_array_Z.sum() - def setFluidInterfaceVarCoord(self, FluidSolver): """ @@ -1515,7 +1516,7 @@ def setFluidInterfaceVarCoord(self, FluidSolver): def setSolidInterfaceLoads(self, SolidSolver, FSI_config): """ Communicates the new solid interface loads to the solid solver. - In case of rigid body motion, calculates the new resultant forces (lift, drag, ...). + Calculates the new resultant forces (lift, drag, ...). """ if self.have_MPI: myid = self.comm.Get_rank() @@ -1534,7 +1535,7 @@ def setSolidInterfaceLoads(self, SolidSolver, FSI_config): FFY = self.fluidLoads_array_Y.sum() FFZ = self.fluidLoads_array_Z.sum() - for iVertex in range(self.nLocalSolidInterfaceNodes): + for iVertex in range(self.nLocalSolidInterfacePhysicalNodes): FX += self.localSolidLoads_array_X[iVertex] FY += self.localSolidLoads_array_Y[iVertex] FZ += self.localSolidLoads_array_Z[iVertex] @@ -1551,17 +1552,16 @@ def setSolidInterfaceLoads(self, SolidSolver, FSI_config): # --- Send the new solid interface loads to the solid solver (on each partition, halo nodes included) --- GlobalIndex = int() localIndex = 0 - if myid in self.solidInterfaceProcessors: - for iVertex in range(self.nLocalSolidInterfaceNodes): - GlobalIndex = SolidSolver.getInterfaceNodeGlobalIndex(self.solidInterfaceIdentifier, iVertex) - if GlobalIndex in self.SolidHaloNodeList[myid].keys(): - pass - else: - Fx = self.localSolidLoads_array_X[localIndex] - Fy = self.localSolidLoads_array_Y[localIndex] - Fz = self.localSolidLoads_array_Z[localIndex] - SolidSolver.applyload(iVertex, Fx, Fy, Fz) - localIndex += 1 + for iVertex in range(self.nLocalSolidInterfaceNodes): + GlobalIndex = SolidSolver.getInterfaceNodeGlobalIndex(self.solidInterfaceIdentifier, iVertex) + if GlobalIndex in self.SolidHaloNodeList[myid].keys(): + pass #TODO here, when the solid solver will run in parallel, we will need to pass the halo loads + else: + Fx = self.localSolidLoads_array_X[localIndex] + Fy = self.localSolidLoads_array_Y[localIndex] + Fz = self.localSolidLoads_array_Z[localIndex] + SolidSolver.applyload(iVertex, Fx, Fy, Fz) + localIndex += 1 def computeSolidInterfaceResidual(self, SolidSolver): """ @@ -1593,14 +1593,16 @@ def computeSolidInterfaceResidual(self, SolidSolver): predDisp_array_X.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF) predDisp_array_Y.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF) predDisp_array_Z.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF) + predDisp_array_X.set(0.0) + predDisp_array_Y.set(0.0) + predDisp_array_Z.set(0.0) - if myid in self.solidSolverProcessors: - for iVertex in range(self.nLocalSolidInterfaceNodes): - predDispx, predDispy, predDispz = SolidSolver.getInterfaceNodeDisp(self.solidInterfaceIdentifier, iVertex) - iGlobalVertex = self.__getGlobalIndex('solid', myid, iVertex) - predDisp_array_X.setValues([iGlobalVertex], predDispx) - predDisp_array_Y.setValues([iGlobalVertex], predDispy) - predDisp_array_Z.setValues([iGlobalVertex], predDispz) + for iVertex in range(self.nLocalSolidInterfaceNodes): + predDispx, predDispy, predDispz = SolidSolver.getInterfaceNodeDisp(self.solidInterfaceIdentifier, iVertex) + iGlobalVertex = self.__getGlobalIndex('solid', myid, iVertex) + predDisp_array_X.setValues([iGlobalVertex], predDispx) + predDisp_array_Y.setValues([iGlobalVertex], predDispy) + predDisp_array_Z.setValues([iGlobalVertex], predDispz) predDisp_array_X.assemblyBegin() predDisp_array_X.assemblyEnd() @@ -1680,7 +1682,7 @@ def setAitkenCoefficient(self, FSI_config): deltaResx_array_Y.setType('seq') deltaResx_array_Z = PETSc.Vec().create() deltaResx_array_Z.setType('seq') - deltaResx_array_X.setSizes(self.nSolidInterfacePhysicalNodes) + deltaResx_array_X.setSizes(self.nSolidInterfacePhysicalNodes) #TODO I think here we should add self.d_RBF, check the sizes deltaResx_array_X.set(0.0) deltaResx_array_Y.setSizes(self.nSolidInterfacePhysicalNodes) deltaResx_array_Y.set(0.0) @@ -1788,9 +1790,7 @@ def displacementPredictor(self, FSI_config , SolidSolver, deltaT): localIndex = 0 for iVertex in range(self.nLocalSolidInterfaceNodes): GlobalIndex = SolidSolver.getInterfaceNodeGlobalIndex(self.solidInterfaceIdentifier, iVertex) - if GlobalIndex in self.SolidHaloNodeList[myid].keys(): - pass - else: + if GlobalIndex not in self.SolidHaloNodeList[myid].keys(): iGlobalVertex = self.__getGlobalIndex('solid', myid, localIndex) velx, vely, velz = SolidSolver.getInterfaceNodeVel(self.solidInterfaceIdentifier, iVertex) velxNm1, velyNm1, velzNm1 = SolidSolver.getInterfaceNodeVelNm1(self.solidInterfaceIdentifier, iVertex) @@ -1931,7 +1931,8 @@ def UnsteadyFSI(self,FSI_config, FluidSolver, SolidSolver): # then to push it back once more to compute the solution at the next time level # Also this is required because in the fluid iteration preprocessor, if we do not update # and step to the next time level, there is a flag "fsi" that will initialise the flow - FluidSolver.Update() + if myid in self.fluidSolverProcessors: + FluidSolver.Update() if myid in self.solidSolverProcessors: SolidSolver.updateSolution() #If no restart @@ -1943,7 +1944,8 @@ def UnsteadyFSI(self,FSI_config, FluidSolver, SolidSolver): self.interpolateSolidPositionOnFluidMesh(FSI_config) self.setFluidInterfaceVarCoord(FluidSolver) self.MPIPrint('\nPerforming static mesh deformation (ALE) of initial mesh...\n') - FluidSolver.SetInitialMesh() # if there is an initial deformation in the solid, it has to be communicated to the fluid solver + if myid in self.fluidSolverProcessors: + FluidSolver.SetInitialMesh() # if there is an initial deformation in the solid, it has to be communicated to the fluid solver self.MPIPrint('\nFSI initial conditions are set') self.MPIPrint('Beginning time integration\n') @@ -1969,18 +1971,20 @@ def UnsteadyFSI(self,FSI_config, FluidSolver, SolidSolver): self.interpolateSolidPositionOnFluidMesh(FSI_config) self.MPIPrint('\nPerforming dynamic mesh deformation (ALE)...\n') self.setFluidInterfaceVarCoord(FluidSolver) - if self.FSIIter == 0: - FluidSolver.Preprocess(TimeIter) # set some parameters before temporal fluid iteration and dynamic mesh update - else: - FluidSolver.DynamicMeshUpdate(TimeIter) + if myid in self.fluidSolverProcessors: + if self.FSIIter == 0: + FluidSolver.Preprocess(TimeIter) # set some parameters before temporal fluid iteration and dynamic mesh update + else: + FluidSolver.DynamicMeshUpdate(TimeIter) # --- Fluid solver call for FSI subiteration --- # self.MPIPrint('\nLaunching fluid solver for one single dual-time iteration...') self.MPIBarrier() - FluidSolver.ResetConvergence() - FluidSolver.Run() - self.MPIBarrier() - FluidSolver.Postprocess() - self.MPIBarrier() + if myid in self.fluidSolverProcessors: + FluidSolver.ResetConvergence() + FluidSolver.Run() + self.MPIBarrier() + FluidSolver.Postprocess() + self.MPIBarrier() # --- Surface fluid loads interpolation and communication --- # if not self.ImposedMotion: @@ -2020,9 +2024,10 @@ def UnsteadyFSI(self,FSI_config, FluidSolver, SolidSolver): self.writeFSIHistory(TimeIter, time, varCoordNorm, FSIConv) # --- Update, monitor and output the fluid solution before the next time step ---# - FluidSolver.Update() - FluidSolver.Monitor(TimeIter) - FluidSolver.Output(TimeIter) + if myid in self.fluidSolverProcessors: + FluidSolver.Update() + FluidSolver.Monitor(TimeIter) + FluidSolver.Output(TimeIter) if TimeIter >= TimeIterTreshold: if myid in self.solidSolverProcessors: @@ -2084,15 +2089,16 @@ def SteadyFSI(self, FSI_config,FluidSolver, SolidSolver): self.MPIPrint('\nLaunching fluid solver for a steady computation...') # --- Fluid solver call for FSI subiteration ---# - FluidSolver.ResetConvergence() #This is setting to zero the convergence in the integrator, important to reset it - # The mesh will be deformed in the context of the preprocessor, there is no need to set the initial - # mesh pushing back the solution to avoid spurious velocities, as the velocity is not computed at all - self.MPIPrint('\nPerforming static mesh deformation...\n') - FluidSolver.Preprocess(0)# This will attempt to always set the initial condition, but there is a flag on the unsteady computation that will avoid it - FluidSolver.Run() - FluidSolver.Postprocess() - FluidSolver.Monitor(0) #This is actually not needed, it only saves the fact that the fluid solver converged innerly or reached max iterations - FluidSolver.Output(0) + if myid in self.fluidSolverProcessors: + FluidSolver.ResetConvergence() #This is setting to zero the convergence in the integrator, important to reset it. + # The mesh will be deformed in the context of the preprocessor, there is no need to set the initial + # mesh pushing back the solution to avoid spurious velocities, as the velocity is not computed at all + self.MPIPrint('\nPerforming static mesh deformation...\n') + FluidSolver.Preprocess(0)# This will attempt to always set the initial condition, but there is a flag on the unsteady computation that will avoid it + FluidSolver.Run() #TODO check how the preprocess work if fsi is false + FluidSolver.Postprocess() + FluidSolver.Monitor(0) #This is actually not needed, it only saves the fact that the fluid solver converged innerly or reached max iterations + FluidSolver.Output(0) # --- Surface fluid loads interpolation and communication ---# if not self.ImposedMotion: diff --git a/SU2_PY/SU2_Nastran/pysu2_nastran.py b/SU2_PY/SU2_Nastran/pysu2_nastran.py index 28a88c3f32eb..cf28f42a5ebb 100644 --- a/SU2_PY/SU2_Nastran/pysu2_nastran.py +++ b/SU2_PY/SU2_Nastran/pysu2_nastran.py @@ -948,11 +948,11 @@ def getInterfaceNodeGlobalIndex(self, markerID, iVertex): return self.markers[markerID][iVertex] - def getInterfaceNodePos(self, markerID, iVertex): + def getInterfaceNodePosInit(self, markerID, iVertex): iPoint = self.markers[markerID][iVertex] - Coord = self.node[iPoint].GetCoord() - return Coord + Coord0 = self.node[iPoint].GetCoord0() + return Coord0 def getInterfaceNodeDisp(self, markerID, iVertex): From d8c1778e79938017f363e947b83087d7ec8baa3b Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Sat, 13 Feb 2021 15:16:48 +0100 Subject: [PATCH 11/32] Revert "Update geo in the mesh deformation was not used" This reverts commit 0dd42ad22c6bb5bec0c8f0dca6d27739071f82cf. --- SU2_CFD/src/drivers/CDriver.cpp | 2 +- SU2_CFD/src/solvers/CMeshSolver.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index e68f59b2ae44..7ee5c3571587 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1334,7 +1334,7 @@ void CDriver::Solver_Restart(CSolver ***solver, CGeometry **geometry, if ((restart || restart_flow) && config->GetDeform_Mesh() && update_geo){ /*--- Always restart with the last state ---*/ val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; - solver[MESH_0][MESH_SOL]->LoadRestart(geometry, solver, config, val_iter); + solver[MESH_0][MESH_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); } /*--- Exit if a restart was requested for a solver that is not available. ---*/ diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 3874f65fd281..5704f6c7bbd1 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -713,7 +713,7 @@ void CMeshSolver::SetDualTime_Mesh(void){ nodes->Set_Solution_time_n(); } -void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int val_iter) { +void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int val_iter, bool val_update_geo) { /*--- Read the restart data from either an ASCII or binary SU2 file. ---*/ From 4b9b8e625ae6de8c939ef4e8d8e0ff60764e6d5e Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Sat, 13 Feb 2021 15:17:04 +0100 Subject: [PATCH 12/32] Revert "Updategeo removed from hpp files also" This reverts commit 2f6635a0b711a8c4493cfa0c4e0597daf8e257c8. --- SU2_CFD/include/solvers/CMeshSolver.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/include/solvers/CMeshSolver.hpp b/SU2_CFD/include/solvers/CMeshSolver.hpp index e1e597323b58..55c14d7aa6f0 100644 --- a/SU2_CFD/include/solvers/CMeshSolver.hpp +++ b/SU2_CFD/include/solvers/CMeshSolver.hpp @@ -148,11 +148,13 @@ class CMeshSolver final : public CFEASolver { * \param[in] solver - Container vector with all of the solvers. * \param[in] config - Definition of the particular problem. * \param[in] val_iter - Current external iteration number. + * \param[in] val_update_geo - Flag for updating coords and grid velocity. */ void LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, - int val_iter) override; + int val_iter, + bool val_update_geo) override; /*! * \brief Load the geometries at the previous time states n and nM1. From 25ee985941bedfa657cf7fe56b5f75348958fa26 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Sat, 13 Feb 2021 15:36:32 +0100 Subject: [PATCH 13/32] Only extract forces when required --- SU2_PY/FSI_tools/FSIInterface.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/SU2_PY/FSI_tools/FSIInterface.py b/SU2_PY/FSI_tools/FSIInterface.py index 62438de8388b..3bb224c8c3e0 100644 --- a/SU2_PY/FSI_tools/FSIInterface.py +++ b/SU2_PY/FSI_tools/FSIInterface.py @@ -1987,13 +1987,12 @@ def UnsteadyFSI(self,FSI_config, FluidSolver, SolidSolver): self.MPIBarrier() # --- Surface fluid loads interpolation and communication --- # - if not self.ImposedMotion: - self.MPIPrint('\nProcessing interface fluid loads...\n') - self.MPIBarrier() - self.getFluidInterfaceNodalForce(FSI_config, FluidSolver) - self.MPIBarrier() if TimeIter > TimeIterTreshold: if not self.ImposedMotion: + self.MPIPrint('\nProcessing interface fluid loads...\n') + self.MPIBarrier() + self.getFluidInterfaceNodalForce(FSI_config, FluidSolver) + self.MPIBarrier() self.interpolateFluidLoadsOnSolidMesh(FSI_config) self.setSolidInterfaceLoads(SolidSolver, FSI_config) From c0c42abe9d409b836601dadb77d03bfc37e95e8d Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Sat, 13 Feb 2021 15:59:23 +0100 Subject: [PATCH 14/32] Fixed size --- SU2_PY/FSI_tools/FSIInterface.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SU2_PY/FSI_tools/FSIInterface.py b/SU2_PY/FSI_tools/FSIInterface.py index 3bb224c8c3e0..aa5423980143 100644 --- a/SU2_PY/FSI_tools/FSIInterface.py +++ b/SU2_PY/FSI_tools/FSIInterface.py @@ -1682,11 +1682,11 @@ def setAitkenCoefficient(self, FSI_config): deltaResx_array_Y.setType('seq') deltaResx_array_Z = PETSc.Vec().create() deltaResx_array_Z.setType('seq') - deltaResx_array_X.setSizes(self.nSolidInterfacePhysicalNodes) #TODO I think here we should add self.d_RBF, check the sizes + deltaResx_array_X.setSizes(self.nSolidInterfacePhysicalNodes + self.d_RBF) deltaResx_array_X.set(0.0) - deltaResx_array_Y.setSizes(self.nSolidInterfacePhysicalNodes) + deltaResx_array_Y.setSizes(self.nSolidInterfacePhysicalNodes + self.d_RBF) deltaResx_array_Y.set(0.0) - deltaResx_array_Z.setSizes(self.nSolidInterfacePhysicalNodes) + deltaResx_array_Z.setSizes(self.nSolidInterfacePhysicalNodes + self.d_RBF) deltaResx_array_Z.set(0.0) # --- Compute the dynamic Aitken coefficient --- From d9ab5923760f752d602cfeb442ca95fc44e0e74d Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Sat, 13 Feb 2021 19:53:11 +0100 Subject: [PATCH 15/32] Removing old setInitialMesh --- SU2_CFD/include/drivers/CDriver.hpp | 15 ---------- SU2_CFD/src/python_wrapper_structure.cpp | 38 ------------------------ 2 files changed, 53 deletions(-) diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index 82bedff1ddc9..ff0570e4c9b6 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -360,11 +360,6 @@ class CDriver { */ virtual void DynamicMeshUpdate(unsigned short val_iZone, unsigned long TimeIter) { } - /*! - * \brief Perform a static mesh deformation, without considering grid velocity. - */ - virtual void StaticMeshUpdate() { } - /*! * \brief Perform a mesh deformation as initial condition. */ @@ -962,16 +957,6 @@ class CFluidDriver : public CDriver { */ void DynamicMeshUpdate(unsigned long TimeIter) override; - /*! - * \brief Perform a static mesh deformation, without considering grid velocity (multiple zone). - */ - void StaticMeshUpdate() override; - - /*! - * \brief Perform a mesh deformation as initial condition (multiple zone). - */ - void SetInitialMesh() override; - /*! * \brief Process the boundary conditions and update the multigrid structure. */ diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 95cbf77ec6d7..770b961899e1 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -759,44 +759,6 @@ void CDriver::ResetConvergence() { } -void CFluidDriver::StaticMeshUpdate() { - - int rank = MASTER_NODE; - -#ifdef HAVE_MPI - MPI_Comm_rank(SU2_MPI::GetComm(), &rank); -#endif - - for(iZone = 0; iZone < nZone; iZone++) { - if(rank == MASTER_NODE) cout << " Deforming the volume grid." << endl; - grid_movement[iZone][INST_0]->SetVolume_Deformation(geometry_container[iZone][INST_0][MESH_0], config_container[iZone], true); - - if(rank == MASTER_NODE) cout << "No grid velocity to be computde : static grid deformation." << endl; - - if(rank == MASTER_NODE) cout << " Updating multigrid structure." << endl; - grid_movement[iZone][INST_0]->UpdateMultiGrid(geometry_container[iZone][INST_0], config_container[iZone]); - } -} - -void CFluidDriver::SetInitialMesh() { - - StaticMeshUpdate(); - - /*--- Propagate the initial deformation to the past ---*/ - //if (!restart) { - for(iZone = 0; iZone < nZone; iZone++) { - for (iMesh = 0; iMesh <= config_container[iZone]->GetnMGLevels(); iMesh++) { - //solver_container[iZone][iMesh][FLOW_SOL]->nodes->Set_Solution_time_n(iPoint); - //solver_container[iZone][iMesh][FLOW_SOL]->nodes->Set_Solution_time_n1(iPoint); - geometry_container[iZone][INST_0][iMesh]->nodes->SetVolume_n(); - geometry_container[iZone][INST_0][iMesh]->nodes->SetVolume_nM1(); - geometry_container[iZone][INST_0][iMesh]->nodes->SetCoord_n(); - geometry_container[iZone][INST_0][iMesh]->nodes->SetCoord_n1(); - } - } - //} -} - void CSinglezoneDriver::SetInitialMesh() { DynamicMeshUpdate(0); From ccd7ac479c5b926635f52d75676d515b0874360b Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Sat, 13 Feb 2021 22:57:51 +0100 Subject: [PATCH 16/32] Compute the interface mapping starting from the undeformed mesh --- SU2_CFD/include/drivers/CDriver.hpp | 8 ++++++++ SU2_CFD/src/python_wrapper_structure.cpp | 10 ++++++++++ SU2_PY/FSI_tools/FSIInterface.py | 8 +++----- 3 files changed, 21 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index ff0570e4c9b6..5177f12eb9c7 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -466,6 +466,14 @@ class CDriver { */ unsigned long GetVertexGlobalIndex(unsigned short iMarker, unsigned long iVertex); + /*! + * \brief Get undeformed coordinates from the mesh solver. + * \param[in] iMarker - Marker identifier. + * \param[in] iVertex - Vertex identifier. + * \return x,y,z coordinates of the vertex. + */ + vector GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex); + /*! * \brief Get the x coordinate of a vertex on a specified marker. * \param[in] iMarker - Marker identifier. diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 770b961899e1..022c5353adbc 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -258,6 +258,16 @@ passivedouble CDriver::GetUnsteady_TimeStep(){ return SU2_TYPE::GetValue(config_container[ZONE_0]->GetTime_Step()); } +vector CDriver::GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex) { + + su2double coord[3] = {0.0}; + + auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); + for (auto iDim = 0 ; iDim < nDim ; iDim++){ + coord[iDim] = solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord(iPoint,iDim); + } +} + passivedouble CDriver::GetVertexCoordX(unsigned short iMarker, unsigned long iVertex) { auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); diff --git a/SU2_PY/FSI_tools/FSIInterface.py b/SU2_PY/FSI_tools/FSIInterface.py index aa5423980143..702ad4d78b0a 100644 --- a/SU2_PY/FSI_tools/FSIInterface.py +++ b/SU2_PY/FSI_tools/FSIInterface.py @@ -565,10 +565,8 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config): # Note that the fluid solver is separated in more processors outside the python script # thus when, from a core, we request for the vertices on the interface, we only obtain # those in that node - GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex) #TODO obtain here the undeformed mesh - posx = FluidSolver.GetVertexCoordX(self.fluidInterfaceIdentifier, iVertex) - posy = FluidSolver.GetVertexCoordY(self.fluidInterfaceIdentifier, iVertex) - posz = FluidSolver.GetVertexCoordZ(self.fluidInterfaceIdentifier, iVertex) + GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex) + posx, posy, posz = FluidSolver.GetInitialMeshCoord(self.fluidInterfaceIdentifier, iVertex) if GlobalIndex not in self.FluidHaloNodeList[myid].keys(): fluidIndexing_temp[GlobalIndex] = self.__getGlobalIndex('fluid', myid, localIndex) self.localFluidInterface_array_X_init[localIndex] = posx @@ -2094,7 +2092,7 @@ def SteadyFSI(self, FSI_config,FluidSolver, SolidSolver): # mesh pushing back the solution to avoid spurious velocities, as the velocity is not computed at all self.MPIPrint('\nPerforming static mesh deformation...\n') FluidSolver.Preprocess(0)# This will attempt to always set the initial condition, but there is a flag on the unsteady computation that will avoid it - FluidSolver.Run() #TODO check how the preprocess work if fsi is false + FluidSolver.Run() FluidSolver.Postprocess() FluidSolver.Monitor(0) #This is actually not needed, it only saves the fact that the fluid solver converged innerly or reached max iterations FluidSolver.Output(0) From a0089c0c6ca258331095fe7dbe4f172186f8ea4f Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Sat, 13 Feb 2021 23:06:16 +0100 Subject: [PATCH 17/32] Loads computed for all the solid markers --- SU2_CFD/include/drivers/CDriver.hpp | 2 +- SU2_CFD/src/python_wrapper_structure.cpp | 13 ++++++++++--- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index 5177f12eb9c7..a6665799ae57 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -472,7 +472,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return x,y,z coordinates of the vertex. */ - vector GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex); + vector GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex); /*! * \brief Get the x coordinate of a vertex on a specified marker. diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 022c5353adbc..4382d020c2e6 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -258,14 +258,21 @@ passivedouble CDriver::GetUnsteady_TimeStep(){ return SU2_TYPE::GetValue(config_container[ZONE_0]->GetTime_Step()); } -vector CDriver::GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex) { +vector CDriver::GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex) { - su2double coord[3] = {0.0}; + vector coord(3,0.0); + vector coord_passive(3, 0.0); auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); for (auto iDim = 0 ; iDim < nDim ; iDim++){ coord[iDim] = solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord(iPoint,iDim); } + + coord_passive[0] = SU2_TYPE::GetValue(coord[0]); + coord_passive[1] = SU2_TYPE::GetValue(coord[1]); + coord_passive[2] = SU2_TYPE::GetValue(coord[2]); + + return coord_passive; } passivedouble CDriver::GetVertexCoordX(unsigned short iMarker, unsigned long iVertex) { @@ -1017,7 +1024,7 @@ vector CDriver::GetFlowLoad(unsigned short iMarker, unsigned long CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]; CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; - if (config_container[ZONE_0]->GetMarker_All_Fluid_Load(iMarker) == YES) { + if (config_container[ZONE_0]->GetSolid_Wall(iMarker)) { FlowLoad[0] = solver->GetVertexTractions(iMarker, iVertex, 0); FlowLoad[1] = solver->GetVertexTractions(iMarker, iVertex, 1); if (geometry->GetnDim() == 3) From a54f6b6037d703519518157eb99f2ea87b3a5516 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Sun, 14 Feb 2021 16:08:47 +0100 Subject: [PATCH 18/32] Tentative modification of test case --- SU2_CFD/include/drivers/CDriver.hpp | 6 - SU2_CFD/src/python_wrapper_structure.cpp | 24 --- .../flatPlate_rigidMotion_Conf.cfg | 186 +----------------- .../launch_flatPlate_rigidMotion.py | 59 +----- 4 files changed, 17 insertions(+), 258 deletions(-) diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index a6665799ae57..c698a85a14ee 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -412,12 +412,6 @@ class CDriver { */ passivedouble Get_LiftCoeff(); - /*! - * \brief Get the moving marker identifier. - * \return Moving marker identifier. - */ - unsigned short GetMovingMarker(); - /*! * \brief Get the number of vertices (halo nodes included) from a specified marker. * \param[in] iMarker - Marker identifier. diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 4382d020c2e6..ff61589d1704 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -178,30 +178,6 @@ passivedouble CDriver::Get_LiftCoeff() { return SU2_TYPE::GetValue(CLift); } -unsigned short CDriver::GetMovingMarker() { - - unsigned short IDtoSend,iMarker, jMarker, Moving; - string Marker_Tag, Moving_Tag; - - IDtoSend = 0; - for (iMarker = 0; iMarker < config_container[ZONE_0]->GetnMarker_All(); iMarker++) { - Moving = config_container[ZONE_0]->GetMarker_All_Moving(iMarker); - if (Moving == YES) { - for (jMarker = 0; jMarkerGetnMarker_Moving(); jMarker++) { - Moving_Tag = config_container[ZONE_0]->GetMarker_Moving_TagBound(jMarker); - Marker_Tag = config_container[ZONE_0]->GetMarker_All_TagBound(iMarker); - if (Marker_Tag == Moving_Tag) { - IDtoSend = iMarker; - break; - } - } - } - } - - return IDtoSend; - -} - unsigned long CDriver::GetNumberVertices(unsigned short iMarker){ return geometry_container[ZONE_0][INST_0][MESH_0]->nVertex[iMarker]; diff --git a/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg b/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg index 0ebe07328e12..d90e636e19da 100644 --- a/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg +++ b/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg @@ -22,10 +22,6 @@ KIND_TURB_MODEL= SST % Mathematical problem (DIRECT, CONTINUOUS_ADJOINT, DISCRETE_ADJOINT) MATH_PROBLEM= DIRECT % -% -% Axisymmetric simulation, only compressible flows (NO, YES) -AXISYMMETRIC= NO -% % Restart solution (NO, YES) RESTART_SOL= NO % @@ -37,11 +33,11 @@ DISCARD_INFILES= NO % % System of measurements (SI, US) % International system of units (SI): ( meters, kilograms, Kelvins, -% Newtons = kg m/s^2, Pascals = N/m^2, +% Newtons = kg m/s^2, Pascals = N/m^2, % Density = kg/m^3, Speed = m/s, % Equiv. Area = m^2 ) -% United States customary units (US): ( inches, slug, Rankines, lbf = slug ft/s^2, -% psf = lbf/ft^2, Density = slug/ft^3, +% United States customary units (US): ( inches, slug, Rankines, lbf = slug ft/s^2, +% psf = lbf/ft^2, Density = slug/ft^3, % Speed = ft/s, Equiv. Area = ft^2 ) SYSTEM_MEASUREMENTS= SI @@ -179,42 +175,12 @@ RESTART_ITER= 0 TIME_ITER=9999 % ----------------------- DYNAMIC MESH DEFINITION -----------------------------% -% Type of dynamic mesh (NONE, RIGID_MOTION, DEFORMING, ROTATING_FRAME, -% MOVING_WALL, STEADY_TRANSLATION, FLUID_STRUCTURE, -% AEROELASTIC, ELASTICITY, EXTERNAL, -% AEROELASTIC_RIGID_MOTION, GUST) -SURFACE_MOVEMENT= FLUID_STRUCTURE -% -% Motion mach number (non-dimensional). Used for initializing a viscous flow -% with the Reynolds number and for computing force coeffs. with dynamic meshes. -MACH_MOTION= 0.03059 -% -% Moving wall boundary marker(s) (NONE = no marker, ignored for RIGID_MOTION) -MARKER_MOVING= ( plate ) -% -% Coordinates of the motion origin -SURFACE_MOTION_ORIGIN= -0.0028 0.0 0.0 -% -% Move Motion Origin for marker moving (1 or 0) -MOVE_MOTION_ORIGIN = 1 - -% ----------------------- BODY FORCE DEFINITION -------------------------------% % -% Apply a body force as a source term (NO, YES) -BODY_FORCE= NO +DEFORM_MESH = YES +MARKER_DEFORM_MESH = (plate) % -% Vector of body force values (BodyForce_X, BodyForce_Y, BodyForce_Z) -BODY_FORCE_VECTOR= ( 0.0, 0.0, 0.0 ) - % -------------------- BOUNDARY CONDITION DEFINITION --------------------------% % -% Euler wall boundary marker(s) (NONE = no marker) -MARKER_EULER= ( NONE ) -% -% Navier-Stokes (no-slip), constant heat flux wall marker(s) (NONE = no marker) -% Format: ( marker name, constant heat flux (J/m^2), ... ) -%MARKER_HEATFLUX= ( plate, 1000.0 ) -% % Navier-Stokes (no-slip), isothermal wall marker(s) (NONE = no marker) % Format: ( marker name, constant wall temperature (K), ... ) MARKER_ISOTHERMAL= ( plate, 293 ) @@ -230,22 +196,6 @@ MARKER_PLOTTING = ( plate ) % Marker(s) of the surface where the non-dimensional coefficients are evaluated. MARKER_MONITORING = ( plate ) % -% Viscous wall markers for which wall functions must be applied. (NONE = no marker) -% Format: ( marker name, wall function type, ... ) -MARKER_WALL_FUNCTIONS= ( plate, NO_WALL_FUNCTION ) -% -% Marker(s) of the surface where custom thermal BC's are defined. -MARKER_PYTHON_CUSTOM = (NONE) -% -% Marker(s) of the surface where obj. func. (design problem) will be evaluated -MARKER_DESIGNING = ( NONE ) -% -% Marker(s) of the surface that is going to be analyzed in detail (massflow, average pressure, distortion, etc) -MARKER_ANALYZE = ( NONE ) -% -% Method to compute the average value in MARKER_ANALYZE (AREA, MASSFLUX). -MARKER_ANALYZE_AVERAGE = MASSFLUX - % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % % Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) @@ -254,39 +204,6 @@ NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES % CFL number (initial value for the adaptive CFL number) CFL_NUMBER= 7.0 % -% Adaptive CFL number (NO, YES) -CFL_ADAPT= NO -% -% Parameters of the adaptive CFL number (factor down, factor up, CFL min value, -% CFL max value ) -CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.25, 50.0 ) -% -% Maximum Delta Time in local time stepping simulations -MAX_DELTA_TIME= 1E6 -% -% Runge-Kutta alpha coefficients -RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) -% -% Objective function in gradient evaluation (DRAG, LIFT, SIDEFORCE, MOMENT_X, -% MOMENT_Y, MOMENT_Z, EFFICIENCY, -% EQUIVALENT_AREA, NEARFIELD_PRESSURE, -% FORCE_X, FORCE_Y, FORCE_Z, THRUST, -% TORQUE, TOTAL_HEATFLUX, -% MAXIMUM_HEATFLUX, INVERSE_DESIGN_PRESSURE, -% INVERSE_DESIGN_HEATFLUX, SURFACE_TOTAL_PRESSURE, -% SURFACE_MASSFLOW, SURFACE_STATIC_PRESSURE, SURFACE_MACH) -% For a weighted sum of objectives: separate by commas, add OBJECTIVE_WEIGHT and MARKER_MONITORING in matching order. -OBJECTIVE_FUNCTION= DRAG -% -% List of weighting values when using more than one OBJECTIVE_FUNCTION. Separate by commas and match with MARKER_MONITORING. -OBJECTIVE_WEIGHT = 1.0 -% -% Reference coefficient (sensitivity) for detecting sharp edges. -REF_SHARP_EDGES= 3.0 -% -% Remove sharp edges from the sensitivity evaluation (NO, YES) -SENS_REMOVE_SHARP= NO - % ----------- SLOPE LIMITER AND DISSIPATION SENSOR DEFINITION -----------------% % % Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations. @@ -301,52 +218,6 @@ SLOPE_LIMITER_FLOW= VENKATAKRISHNAN % Required for 2nd order upwind schemes (NO, YES) MUSCL_TURB= NO % -% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, -% BARTH_JESPERSEN, VAN_ALBADA_EDGE) -SLOPE_LIMITER_TURB= VENKATAKRISHNAN -% -% Monotonic Upwind Scheme for Conservation Laws (TVD) in the adjoint flow equations. -% Required for 2nd order upwind schemes (NO, YES) -MUSCL_ADJFLOW= YES -% -% Slope limiter (NONE, VENKATAKRISHNAN, BARTH_JESPERSEN, VAN_ALBADA_EDGE, -% SHARP_EDGES, WALL_DISTANCE) -SLOPE_LIMITER_ADJFLOW= VENKATAKRISHNAN -% -% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence adjoint equations. -% Required for 2nd order upwind schemes (NO, YES) -MUSCL_ADJTURB= NO -% -% Slope limiter (NONE, VENKATAKRISHNAN, BARTH_JESPERSEN, VAN_ALBADA_EDGE) -SLOPE_LIMITER_ADJTURB= VENKATAKRISHNAN -% -% Coefficient for the Venkat's limiter (upwind scheme). A larger values decrease -% the extent of limiting, values approaching zero cause -% lower-order approximation to the solution (0.05 by default) -VENKAT_LIMITER_COEFF= 0.05 -% -% Coefficient for the adjoint sharp edges limiter (3.0 by default). -ADJ_SHARP_LIMITER_COEFF= 3.0 -% -% Freeze the value of the limiter after a number of iterations -LIMITER_ITER= 999999 -% -% 1st order artificial dissipation coefficients for -% the Lax–Friedrichs method ( 0.15 by default ) -LAX_SENSOR_COEFF= 0.15 -% -% 2nd and 4th order artificial dissipation coefficients for -% the JST method ( 0.5, 0.02 by default ) -JST_SENSOR_COEFF= ( 0.5, 0.02 ) -% -% 1st order artificial dissipation coefficients for -% the adjoint Lax–Friedrichs method ( 0.15 by default ) -ADJ_LAX_SENSOR_COEFF= 0.15 -% -% 2nd, and 4th order artificial dissipation coefficients for -% the adjoint JST method ( 0.5, 0.02 by default ) -ADJ_JST_SENSOR_COEFF= ( 0.5, 0.02 ) - % ------------------------ LINEAR SOLVER DEFINITION ---------------------------% % % Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER) @@ -407,9 +278,6 @@ CONV_NUM_METHOD_TURB= SCALAR_UPWIND % Time discretization (EULER_IMPLICIT) TIME_DISCRE_TURB= EULER_IMPLICIT % -% Reduction factor of the CFL coefficient in the turbulence problem -CFL_REDUCTION_TURB= 1.0 -% % Relaxation coefficient % ------------------------ GRID DEFORMATION PARAMETERS ------------------------% @@ -432,17 +300,10 @@ DEFORM_CONSOLE_OUTPUT= YES % Minimum residual criteria for the linear solver convergence of grid deformation DEFORM_LINEAR_SOLVER_ERROR= 1E-14 % -% Deformation coefficient (linear elasticity limits from -1.0 to 0.5, a larger -% value is also possible) -DEFORM_COEFF = 1E6 % % Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, % WALL_DISTANCE, CONSTANT_STIFFNESS) -DEFORM_STIFFNESS_TYPE= INVERSE_VOLUME -% -% Deform the grid only close to the surface. It is possible to specify how much -% of the volumetric grid is going to be deformed in meters or inches (1E6 by default) -DEFORM_LIMIT = 1E6 +DEFORM_STIFFNESS_TYPE= WALL_DISTANCE % --------------------------- CONVERGENCE PARAMETERS --------------------------% % Convergence criteria (CAUCHY, RESIDUAL) @@ -472,55 +333,22 @@ MESH_FILENAME= 2D_FlatPlate_Rounded.su2 % Mesh input file format (SU2, CGNS) MESH_FORMAT= SU2 % -% Mesh output file -MESH_OUT_FILENAME= mesh_out.su2 -% -% Restart flow input file -SOLUTION_FILENAME= restart_flow.dat -% -% Restart adjoint input file -SOLUTION_ADJ_FILENAME= solution_adj.dat % % Output file format (TECPLOT, TECPLOT_BINARY, PARAVIEW, % FIELDVIEW, FIELDVIEW_BINARY) -TABULAR_FORMAT= CSV +TABULAR_FORMAT= TECPLOT % % Output file convergence history (w/o extension) CONV_FILENAME= history % -% Output file with the forces breakdown -BREAKDOWN_FILENAME= forces_breakdown.dat -% -% Output file restart flow -RESTART_FILENAME= restart_flow.dat -% -% Output file restart adjoint -RESTART_ADJ_FILENAME= restart_adj.dat -% % Output file flow (w/o extension) variables VOLUME_FILENAME= flow % -% Output file adjoint (w/o extension) variables -VOLUME_ADJ_FILENAME= adjoint -% -% Output Objective function -VALUE_OBJFUNC_FILENAME= of_eval.dat -% -% Output objective function gradient (using continuous adjoint) -GRAD_OBJFUNC_FILENAME= of_grad.dat -% % Output file surface flow coefficient (w/o extension) SURFACE_FILENAME= surface_flow % -% Output file surface adjoint coefficient (w/o extension) -SURFACE_ADJ_FILENAME= surface_adjoint -% % Writing solution file frequency for physical time steps (dual time) OUTPUT_WRT_FREQ= 3 % -% -% Read binary restart files (YES, NO) -READ_BINARY_RESTART= YES -% % Screen output SCREEN_OUTPUT= (TIME_ITER, INNER_ITER, RMS_DENSITY, RMS_TKE, RMS_DISSIPATION, LIFT, DRAG) diff --git a/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py b/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py index 6c6a97a8f1f9..f8280f3caaed 100755 --- a/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py +++ b/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py @@ -46,7 +46,7 @@ import numpy as np # ------------------------------------------------------------------- -# Main +# Main # ------------------------------------------------------------------- def main(): @@ -54,38 +54,12 @@ def main(): # Command line options parser=OptionParser() parser.add_option("-f", "--file", dest="filename", help="Read config from FILE", metavar="FILE") - parser.add_option("--nDim", dest="nDim", default=2, help="Define the number of DIMENSIONS", - metavar="DIMENSIONS") - parser.add_option("--nZone", dest="nZone", default=1, help="Define the number of ZONES", - metavar="ZONES") parser.add_option("--parallel", action="store_true", help="Specify if we need to initialize MPI", dest="with_MPI", default=False) - parser.add_option("--fsi", dest="fsi", default="False", help="Launch the FSI driver", metavar="FSI") - - parser.add_option("--fem", dest="fem", default="False", help="Launch the FEM driver (General driver)", metavar="FEM") - - parser.add_option("--harmonic_balance", dest="harmonic_balance", default="False", - help="Launch the Harmonic Balance (HB) driver", metavar="HB") - - parser.add_option("--poisson_equation", dest="poisson_equation", default="False", - help="Launch the poisson equation driver (General driver)", metavar="POIS_EQ") - - parser.add_option("--wave_equation", dest="wave_equation", default="False", - help="Launch the wave equation driver (General driver)", metavar="WAVE_EQ") - - parser.add_option("--heat_equation", dest="heat_equation", default="False", - help="Launch the heat equation driver (General driver)", metavar="HEAT_EQ") - (options, args) = parser.parse_args() - options.nDim = int( options.nDim ) - options.nZone = int( options.nZone ) - options.fsi = options.fsi.upper() == 'TRUE' - options.fem = options.fem.upper() == 'TRUE' - options.harmonic_balance = options.harmonic_balance.upper() == 'TRUE' - options.poisson_equation = options.poisson_equation.upper() == 'TRUE' - options.wave_equation = options.wave_equation.upper() == 'TRUE' - options.heat_equation = options.heat_equation.upper() == 'TRUE' + options.nDim = int(2) + options.nZone = int(1) # Import mpi4py for parallel run if options.with_MPI == True: @@ -98,13 +72,6 @@ def main(): # Initialize the corresponding driver of SU2, this includes solver preprocessing try: - if (options.nZone == 1) and ( options.fem or options.poisson_equation or options.wave_equation or options.heat_equation ): - SU2Driver = pysu2.CSinglezoneDriver(options.filename, options.nZone, comm); - elif options.harmonic_balance: - SU2Driver = pysu2.CHBDriver(options.filename, options.nZone, comm); - elif (options.nZone == 2) and (options.fsi): - SU2Driver = pysu2.CFSIDriver(options.filename, options.nZone, comm); - else: SU2Driver = pysu2.CSinglezoneDriver(options.filename, options.nZone, comm); except TypeError as exception: print('A TypeError occured in pysu2.CDriver : ',exception) @@ -119,7 +86,7 @@ def main(): MovingMarker = 'plate' #specified by the user # Get all the tags with the moving option - MovingMarkerList = SU2Driver.GetAllMovingMarkersTag() + MovingMarkerList = SU2Driver.GetAllDeformMeshMarkersTag() # Get all the markers defined on this rank and their associated indices. allMarkerIDs = SU2Driver.GetAllBoundaryMarkers() @@ -147,9 +114,9 @@ def main(): # Extract the initial position of each node on the moving marker CoordX = np.zeros(nVertex_MovingMarker) CoordY = np.zeros(nVertex_MovingMarker) + CoordZ = np.zeros(nVertex_MovingMarker) for iVertex in range(nVertex_MovingMarker): - CoordX[iVertex] = SU2Driver.GetVertexCoordX(MovingMarkerID, iVertex) - CoordY[iVertex] = SU2Driver.GetVertexCoordY(MovingMarkerID, iVertex) + CoordX[iVertex], CoordY[iVertex], CoordZ[iVertex] = SU2Driver.GetInitialMeshCoord(MovingMarkerID, iVertex) # Time loop is defined in Python so that we have acces to SU2 functionalities at each time step if rank == 0: @@ -162,18 +129,15 @@ def main(): # Define the rigid body displacement and set the new coords of each node on the marker d_y = 0.0175*sin(2*pi*time) for iVertex in range(nVertex_MovingMarker): - newCoordX = CoordX[iVertex] - newCoordY = CoordY[iVertex] + d_y - SU2Driver.SetVertexCoordX(MovingMarkerID, iVertex, newCoordX) - SU2Driver.SetVertexCoordY(MovingMarkerID, iVertex, newCoordY) - SU2Driver.SetVertexCoordZ(MovingMarkerID, iVertex, 0.0) - SU2Driver.SetVertexVarCoord(MovingMarkerID, iVertex) + SU2Driver.SetMeshDisplacement(MovingMarkerID, int(iVertex), 0.0, d_y, 0.0) # Time iteration preprocessing SU2Driver.Preprocess(TimeIter) # Run one time iteration (e.g. dual-time) SU2Driver.Run() # Update the solver for the next time iteration SU2Driver.Update() + # Postprocess the solver + SU2Driver.Postprocess() # Monitor the solver and output solution to file if required stopCalc = SU2Driver.Monitor(TimeIter) SU2Driver.Output(TimeIter) @@ -183,9 +147,6 @@ def main(): TimeIter += 1 time += deltaT - # Postprocess the solver and exit cleanly - SU2Driver.Postprocessing() - if SU2Driver != None: del SU2Driver @@ -195,4 +156,4 @@ def main(): # this is only accessed if running from command prompt if __name__ == '__main__': - main() + main() From 52a2a69390c3ee2804b2b9ba7969ceba9b0280e7 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Sun, 14 Feb 2021 17:14:48 +0100 Subject: [PATCH 19/32] Introduced start time for forced motion --- SU2_PY/SU2_Nastran/pysu2_nastran.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/SU2_PY/SU2_Nastran/pysu2_nastran.py b/SU2_PY/SU2_Nastran/pysu2_nastran.py index cf28f42a5ebb..6f6144c4bc28 100644 --- a/SU2_PY/SU2_Nastran/pysu2_nastran.py +++ b/SU2_PY/SU2_Nastran/pysu2_nastran.py @@ -50,12 +50,14 @@ def __init__(self,time0,tipo,parameters): self.bias = parameters[0] self.amplitude = parameters[1] self.frequency = parameters[2] + self.timeStart = parameters[3] elif self.tipo == "BLENDED_STEP": self.kmax = parameters[0] self.vinf = parameters[1] self.lref = parameters[2] self.amplitude = parameters[3] + self.timeStart = parameters[4] self.tmax = 2*pi/self.kmax*self.lref/self.vinf self.omega0 = 1/2*self.kmax @@ -64,35 +66,41 @@ def __init__(self,time0,tipo,parameters): def GetDispl(self,time): - time = time - self.time0 + time = time - self.time0 - self.timeStart if self.tipo == "SINUSOIDAL": return self.bias+self.amplitude*sin(2*pi*self.frequency*time) if self.tipo == "BLENDED_STEP": - if time < self.tmax: + if time < 0: + return 0.0 + elif time < self.tmax: return self.amplitude/2.0*(1.0-cos(self.omega0*time*self.vinf/self.lref)) return self.amplitude def GetVel(self,time): - time = time - self.time0 + time = time - self.time0 - self.timeStart if self.tipo == "SINUSOIDAL": return self.amplitude*cos(2*pi*self.frequency*time)*2*pi*self.frequency if self.tipo == "BLENDED_STEP": - if time < self.tmax: + if time < 0: + return 0.0 + elif time < self.tmax: return self.amplitude/2.0*sin(self.omega0*time*self.vinf/self.lref)*(self.omega0*self.vinf/self.lref) return 0.0 def GetAcc(self,time): - time = time - self.time0 + time = time - self.time0 - self.timeStart if self.tipo == "SINUSOIDAL": return -self.amplitude*sin(2*pi*self.frequency*time)*(2*pi*self.frequency)**2 if self.tipo == "BLENDED_STEP": - if time < self.tmax: + if time < 0: + return 0.0 + elif time < self.tmax: return self.amplitude/2.0*cos(self.omega0*time*self.vinf/self.lref)*(self.omega0*self.vinf/self.lref)**2 return 0.0 From 727c129bf276dd276244d8b7798c714890efffe1 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 15 Feb 2021 09:35:15 +0100 Subject: [PATCH 20/32] Final version of new test case --- .../flatPlate_rigidMotion_Conf.cfg | 128 +----------------- .../launch_flatPlate_rigidMotion.py | 7 +- 2 files changed, 11 insertions(+), 124 deletions(-) diff --git a/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg b/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg index d90e636e19da..8d28b9561710 100644 --- a/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg +++ b/TestCases/py_wrapper/flatPlate_rigidMotion/flatPlate_rigidMotion_Conf.cfg @@ -25,22 +25,6 @@ MATH_PROBLEM= DIRECT % Restart solution (NO, YES) RESTART_SOL= NO % -% Discard the data storaged in the solution and geometry files -% e.g. AOA, dCL/dAoA, dCD/dCL, iter, etc. -% Note that AoA in the solution and geometry files is critical -% to aero design using AoA as a variable. (NO, YES) -DISCARD_INFILES= NO -% -% System of measurements (SI, US) -% International system of units (SI): ( meters, kilograms, Kelvins, -% Newtons = kg m/s^2, Pascals = N/m^2, -% Density = kg/m^3, Speed = m/s, -% Equiv. Area = m^2 ) -% United States customary units (US): ( inches, slug, Rankines, lbf = slug ft/s^2, -% psf = lbf/ft^2, Density = slug/ft^3, -% Speed = ft/s, Equiv. Area = ft^2 ) -SYSTEM_MEASUREMENTS= SI - % -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% % % Mach number (non-dimensional, based on the free-stream values) @@ -49,9 +33,6 @@ MACH_NUMBER= 0.03059 % Angle of attack (degrees, only for compressible flows) AOA= 0.0 % -% Side-slip angle (degrees, only for compressible flows) -SIDESLIP_ANGLE= 0.0 -% % Init option to choose between Reynolds (default) or thermodynamics quantities % for initializing the solution (REYNOLDS, TD_CONDITIONS) INIT_OPTION= REYNOLDS @@ -60,9 +41,6 @@ INIT_OPTION= REYNOLDS % initializing the solution (TEMPERATURE_FS, DENSITY_FS) FREESTREAM_OPTION= TEMPERATURE_FS % -% Free-stream pressure (101325.0 N/m^2, 2116.216 psf by default) -FREESTREAM_PRESSURE= 101325.0 -% % Free-stream temperature (288.15 K, 518.67 R by default) FREESTREAM_TEMPERATURE= 293.15 % @@ -71,18 +49,7 @@ REYNOLDS_NUMBER= 24407.25244 % % Reynolds length (1 m, 1 inch by default) REYNOLDS_LENGTH= 0.035 - -% -------------------- INCOMPRESSIBLE FREE-STREAM DEFINITION ------------------% -% -% Free-stream density (1.2886 Kg/m^3, 0.0025 slug/ft^3 by default) -FREESTREAM_DENSITY= 1.204 -% -% Free-stream velocity (1.0 m/s, 1.0 ft/s by default) -FREESTREAM_VELOCITY= ( 1.0, 0.00, 0.00 ) % -% Free-stream viscosity (1.853E-5 N s/m^2, 3.87E-7 lbf s/ft^2 by default) -FREESTREAM_VISCOSITY= 1.82E-5 - % ---------------------- REFERENCE VALUE DEFINITION ---------------------------% % % Reference origin for moment computation (m or in) @@ -97,60 +64,25 @@ REF_LENGTH= 0.035 % calculation) (m^2 or in^2) REF_AREA= 0.035 % -% Aircraft semi-span (0 implies automatic calculation) (m or in) -SEMI_SPAN= 0.0 -% % Flow non-dimensionalization (DIMENSIONAL, FREESTREAM_PRESS_EQ_ONE, % FREESTREAM_VEL_EQ_MACH, FREESTREAM_VEL_EQ_ONE) REF_DIMENSIONALIZATION= DIMENSIONAL - +% % ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------% % % Different gas model (STANDARD_AIR, IDEAL_GAS, VW_GAS, PR_GAS) FLUID_MODEL= STANDARD_AIR % -% Ratio of specific heats (1.4 default and the value is hardcoded -% for the model STANDARD_AIR) -GAMMA_VALUE= 1.4 -% -% Specific gas constant (287.058 J/kg*K default and this value is hardcoded -% for the model STANDARD_AIR) -GAS_CONSTANT= 287.058 -% -% Critical Temperature (131.00 K by default) -CRITICAL_TEMPERATURE= 131.00 -% -% Critical Pressure (3588550.0 N/m^2 by default) -CRITICAL_PRESSURE= 3588550.0 -% -% Acentri factor (0.035 (air)) -ACENTRIC_FACTOR= 0.035 - % --------------------------- VISCOSITY MODEL ---------------------------------% % % Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY). VISCOSITY_MODEL= SUTHERLAND % -% Molecular Viscosity that would be constant (1.716E-5 by default) -MU_CONSTANT= 1.716E-5 -% -% Sutherland Viscosity Ref (1.716E-5 default value for AIR SI) -MU_REF= 1.716E-5 -% -% Sutherland Temperature Ref (273.15 K default value for AIR SI) -MU_T_REF= 273.15 -% -% Sutherland constant (110.4 default value for AIR SI) -SUTHERLAND_CONSTANT= 110.4 - % --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% % % Conductivity model (CONSTANT_CONDUCTIVITY, CONSTANT_PRANDTL). CONDUCTIVITY_MODEL= CONSTANT_PRANDTL % -% Molecular Thermal Conductivity that would be constant (0.0257 by default) -KT_CONSTANT= 0.0257 - % ------------------------- UNSTEADY SIMULATION -------------------------------% % TIME_DOMAIN=YES @@ -163,17 +95,13 @@ TIME_STEP= 0.003 % % Total Physical Time for dual time stepping simulations (s) MAX_TIME= 1.0 -% -% Unsteady Courant-Friedrichs-Lewy number of the finest grid -UNST_CFL_NUMBER= 0.0 +TIME_ITER = 9999 % % Number of internal iterations (dual time method) INNER_ITER= 10 % % Iteration number to begin unsteady restarts RESTART_ITER= 0 - -TIME_ITER=9999 % ----------------------- DYNAMIC MESH DEFINITION -----------------------------% % DEFORM_MESH = YES @@ -208,11 +136,7 @@ CFL_NUMBER= 7.0 % % Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations. % Required for 2nd order upwind schemes (NO, YES) -MUSCL_FLOW= YES -% -% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, -% BARTH_JESPERSEN, VAN_ALBADA_EDGE) -SLOPE_LIMITER_FLOW= VENKATAKRISHNAN +MUSCL_FLOW= NO % % Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations. % Required for 2nd order upwind schemes (NO, YES) @@ -238,26 +162,8 @@ LINEAR_SOLVER_ITER= 10 % -------------------------- MULTIGRID PARAMETERS -----------------------------% % % Multi-grid levels (0 = no multi-grid) -MGLEVEL= 3 -% -% Multi-grid cycle (V_CYCLE, W_CYCLE, FULLMG_CYCLE) -MGCYCLE= W_CYCLE -% -% Multi-grid pre-smoothing level -MG_PRE_SMOOTH= ( 1, 2, 3, 3 ) -% -% Multi-grid post-smoothing level -MG_POST_SMOOTH= ( 0, 0, 0, 0 ) +MGLEVEL= 0 % -% Jacobi implicit smoothing of the correction -MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 ) -% -% Damping factor for the residual restriction -MG_DAMP_RESTRICTION= 0.75 -% -% Damping factor for the correction prolongation -MG_DAMP_PROLONGATION= 0.75 - % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % % Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC, @@ -268,8 +174,7 @@ CONV_NUM_METHOD_FLOW= JST % Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT) TIME_DISCRE_FLOW= EULER_IMPLICIT % -% Relaxation coefficient - +% % -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% % % Convective numerical method (SCALAR_UPWIND) @@ -278,8 +183,6 @@ CONV_NUM_METHOD_TURB= SCALAR_UPWIND % Time discretization (EULER_IMPLICIT) TIME_DISCRE_TURB= EULER_IMPLICIT % -% Relaxation coefficient - % ------------------------ GRID DEFORMATION PARAMETERS ------------------------% % % Linear solver or smoother for implicit formulations (FGMRES, RESTARTED_FGMRES, BCGSTAB) @@ -304,27 +207,7 @@ DEFORM_LINEAR_SOLVER_ERROR= 1E-14 % Type of element stiffness imposed for FEA mesh deformation (INVERSE_VOLUME, % WALL_DISTANCE, CONSTANT_STIFFNESS) DEFORM_STIFFNESS_TYPE= WALL_DISTANCE - -% --------------------------- CONVERGENCE PARAMETERS --------------------------% -% Convergence criteria (CAUCHY, RESIDUAL) -% -CONV_CRITERIA= CAUCHY -% -% -% Min value of the residual (log10 of the residual) -CONV_RESIDUAL_MINVAL= -10 % -% Start convergence criteria at iteration number -CONV_STARTITER= 4 -% -% Number of elements to apply the criteria -CONV_CAUCHY_ELEMS= 10 -% -% Epsilon to control the series convergence -CONV_CAUCHY_EPS= 1E-6 -% -% - % ------------------------- INPUT/OUTPUT INFORMATION --------------------------% % % Mesh input file @@ -352,3 +235,4 @@ OUTPUT_WRT_FREQ= 3 % % Screen output SCREEN_OUTPUT= (TIME_ITER, INNER_ITER, RMS_DENSITY, RMS_TKE, RMS_DISSIPATION, LIFT, DRAG) +HISTORY_OUTPUT=(ITER,RMS_RES,AERO_COEFF) diff --git a/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py b/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py index f8280f3caaed..a611e6554d79 100755 --- a/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py +++ b/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py @@ -134,10 +134,10 @@ def main(): SU2Driver.Preprocess(TimeIter) # Run one time iteration (e.g. dual-time) SU2Driver.Run() - # Update the solver for the next time iteration - SU2Driver.Update() # Postprocess the solver SU2Driver.Postprocess() + # Update the solver for the next time iteration + SU2Driver.Update() # Monitor the solver and output solution to file if required stopCalc = SU2Driver.Monitor(TimeIter) SU2Driver.Output(TimeIter) @@ -147,6 +147,9 @@ def main(): TimeIter += 1 time += deltaT + # Postprocess the solver and exit cleanly + SU2Driver.Postprocessing() + if SU2Driver != None: del SU2Driver From 905ef8d6664525b1eff21d8aa83b2239747eefc7 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 15 Feb 2021 10:50:56 +0100 Subject: [PATCH 21/32] Updated values in regression --- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index c7ec41d0a580..98d91e311d3d 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1298,7 +1298,7 @@ def main(): pywrapper_rigidMotion.cfg_dir = "py_wrapper/flatPlate_rigidMotion" pywrapper_rigidMotion.cfg_file = "flatPlate_rigidMotion_Conf.cfg" pywrapper_rigidMotion.test_iter = 5 - pywrapper_rigidMotion.test_vals = [-1.614165, 2.242641, -0.038307, 0.173866] + pywrapper_rigidMotion.test_vals = [-1.614170, 2.242953, 0.350050, 0.093137] pywrapper_rigidMotion.su2_exec = "mpirun -np 2 python launch_flatPlate_rigidMotion.py --parallel -f" pywrapper_rigidMotion.timeout = 1600 pywrapper_rigidMotion.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 55ffcd722dca..77dccd00fe71 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -86,8 +86,8 @@ def main(): viscwedge.timeout = 1600 viscwedge.new_output = True viscwedge.tol = 0.00001 - test_list.append(viscwedge) - + test_list.append(viscwedge) + ######################### ## Compressible Euler ### ######################### @@ -1877,7 +1877,7 @@ def main(): pywrapper_rigidMotion.cfg_dir = "py_wrapper/flatPlate_rigidMotion" pywrapper_rigidMotion.cfg_file = "flatPlate_rigidMotion_Conf.cfg" pywrapper_rigidMotion.test_iter = 5 - pywrapper_rigidMotion.test_vals = [-1.614167, 2.242632, -0.037871, 0.173912] + pywrapper_rigidMotion.test_vals = [-1.614170, 2.242953, 0.350050, 0.093137] pywrapper_rigidMotion.su2_exec = "python launch_flatPlate_rigidMotion.py -f" pywrapper_rigidMotion.new_output = True pywrapper_rigidMotion.timeout = 1600 From 6aca633731d46634b98e234c82e44456dfa2a994 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 15 Feb 2021 10:51:20 +0100 Subject: [PATCH 22/32] Complete reorganisation of the interface --- SU2_CFD/include/drivers/CDriver.hpp | 205 +------ SU2_CFD/src/python_wrapper_structure.cpp | 534 +++++------------- .../launch_unsteady_CHT_FlatPlate.py | 48 +- 3 files changed, 162 insertions(+), 625 deletions(-) diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index c698a85a14ee..617c23d7316c 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -97,10 +97,6 @@ class CDriver { vector > > interpolator_container; /*!< \brief Definition of the interpolation method between non-matching discretizations of the interface. */ CInterface ***interface_container; /*!< \brief Definition of the interface of information and physics. */ - su2double PyWrapVarCoord[3], /*!< \brief This is used to store the VarCoord of each vertex. */ - PyWrapNodalForce[3], /*!< \brief This is used to store the force at each vertex. */ - PyWrapNodalForceDensity[3], /*!< \brief This is used to store the force density at each vertex. */ - PyWrapNodalHeatFlux[3]; /*!< \brief This is used to store the heat flux at each vertex. */ bool dry_run; /*!< \brief Flag if SU2_CFD was started as dry-run via "SU2_CFD -d .cfg" */ public: @@ -368,7 +364,7 @@ class CDriver { /*! * \brief Process the boundary conditions and update the multigrid structure. */ - virtual void BoundaryConditionsUpdate() { } + void BoundaryConditionsUpdate(); /*! * \brief Get the total drag. @@ -468,118 +464,6 @@ class CDriver { */ vector GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex); - /*! - * \brief Get the x coordinate of a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return x coordinate of the vertex. - */ - passivedouble GetVertexCoordX(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the y coordinate of a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return y coordinate of the vertex. - */ - passivedouble GetVertexCoordY(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the z coordinate of a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return z coordinate of the vertex. - */ - passivedouble GetVertexCoordZ(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Compute the total force (pressure and shear stress) at a vertex on a specified marker (3 components). - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return True if the vertex is a halo node (non physical force). - */ - bool ComputeVertexForces(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the x component of the force at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return x component of the force at the vertex. - */ - passivedouble GetVertexForceX(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the y component of the force at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return y component of the force at the vertex. - */ - passivedouble GetVertexForceY(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the z component of the force at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return z component of the force at the vertex. - */ - passivedouble GetVertexForceZ(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the x component of the force density at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return x component of the force density at the vertex. - */ - passivedouble GetVertexForceDensityX(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the y component of the force density at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return y component of the force density at the vertex. - */ - passivedouble GetVertexForceDensityY(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the z component of the force density at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return z component of the force density at the vertex. - */ - passivedouble GetVertexForceDensityZ(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Set the x coordinate of a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \param[in] newPosX - New x coordinate of the vertex. - */ - void SetVertexCoordX(unsigned short iMarker, unsigned long iVertex, passivedouble newPosX); - - /*! - * \brief Set the y coordinate of a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \param[in] newPosY - New y coordinate of the vertex. - */ - void SetVertexCoordY(unsigned short iMarker, unsigned long iVertex, passivedouble newPosY); - - /*! - * \brief Set the z coordinate of a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \param[in] newPosZ - New z coordinate of the vertex. - */ - void SetVertexCoordZ(unsigned short iMarker, unsigned long iVertex, passivedouble newPosZ); - - /*! - * \brief Set the VarCoord of a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return Norm of the VarCoord. - */ - passivedouble SetVertexVarCoord(unsigned short iMarker, unsigned long iVertex); - /*! * \brief Get the temperature at a vertex on a specified marker. * \param[in] iMarker - Marker identifier. @@ -597,36 +481,12 @@ class CDriver { void SetVertexTemperature(unsigned short iMarker, unsigned long iVertex, passivedouble val_WallTemp); /*! - * \brief Compute the heat flux at a vertex on a specified marker (3 components). + * \brief Get the heat flux at a vertex on a specified marker (3 components). * \param[in] iMarker - Marker identifier. * \param[in] iVertex - Vertex identifier. * \return True if the vertex is a halo node. */ - bool ComputeVertexHeatFluxes(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the x component of the heat flux at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return x component of the heat flux at the vertex. - */ - passivedouble GetVertexHeatFluxX(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the y component of the heat flux at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return y component of the heat flux at the vertex. - */ - passivedouble GetVertexHeatFluxY(unsigned short iMarker, unsigned long iVertex); - - /*! - * \brief Get the z component of the heat flux at a vertex on a specified marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return z component of the heat flux at the vertex. - */ - passivedouble GetVertexHeatFluxZ(unsigned short iMarker, unsigned long iVertex); + vector GetVertexHeatFluxes(unsigned short iMarker, unsigned long iVertex); /*! * \brief Get the wall normal component of the heat flux at a vertex on a specified marker. @@ -675,24 +535,12 @@ class CDriver { */ vector GetAllBoundaryMarkersTag(); - /*! - * \brief Get all the moving boundary markers tags. - * \return List of moving boundary markers tags. - */ - vector GetAllMovingMarkersTag(); - /*! * \brief Get all the deformable boundary marker tags. * \return List of deformable boundary markers tags. */ vector GetAllDeformMeshMarkersTag(); - /*! - * \brief Get all the fluid load boundary marker tags. - * \return List of fluid load boundary markers tags. - */ - vector GetAllFluidLoadMarkersTag(); - /*! * \brief Get all the heat transfer boundary markers tags. * \return List of heat transfer boundary markers tags. @@ -816,14 +664,6 @@ class CDriver { void SetSourceTerm_DispAdjoint(unsigned short iMarker, unsigned long iVertex, passivedouble val_AdjointX, passivedouble val_AdjointY, passivedouble val_AdjointZ); - /*! - * \brief Get the undeformed mesh coordinates - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \return Undeformed Vertex Coordinates - */ - vector GetVertex_UndeformedCoord(unsigned short iMarker, unsigned long iVertex); - /*! * \brief Set the position of the heat source. * \param[in] alpha - Angle of rotation respect to Z axis. @@ -959,50 +799,11 @@ class CFluidDriver : public CDriver { */ void DynamicMeshUpdate(unsigned long TimeIter) override; - /*! - * \brief Process the boundary conditions and update the multigrid structure. - */ - void BoundaryConditionsUpdate() override; - /*! * \brief Transfer data among different zones (multiple zone). */ void Transfer_Data(unsigned short donorZone, unsigned short targetZone); - /*! - * \brief Set the total temperature of a vertex on a specified inlet marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \param[in] val_Ttotal - Value of the total (stagnation) temperature. - */ - void SetVertexTtotal(unsigned short iMarker, unsigned long iVertex, passivedouble val_Ttotal); - - /*! - * \brief Set the total pressure of a vertex on a specified inlet marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \param[in] val_Ptotal - Value of the total (stagnation) pressure. - */ - void SetVertexPtotal(unsigned short iMarker, unsigned long iVertex, passivedouble val_Ptotal); - - /*! - * \brief Set the flow direction of a vertex on a specified inlet marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \param[in] iDim - Index of the flow direction unit vector - * \param[in] val_FlowDir - Component of a unit vector representing the flow direction - */ - void SetVertexFlowDir(unsigned short iMarker, unsigned long iVertex, unsigned short iDim, passivedouble val_FlowDir); - - /*! - * \brief Set a turbulence variable on a specified inlet marker. - * \param[in] iMarker - Marker identifier. - * \param[in] iVertex - Vertex identifier. - * \param[in] iDim - Index of the turbulence variable (i.e. k is 0 in SST) - * \param[in] val_turb_var - Value of the turbulence variable to be used. - */ - void SetVertexTurbVar(unsigned short iMarker, unsigned long iVertex, unsigned short iDim, passivedouble val_tub_var); - }; diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index ff61589d1704..383ddc2eec9d 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -59,22 +59,13 @@ void CDriver::PythonInterface_Preprocessing(CConfig **config, CGeometry ****geom } } } - /*--- Initialize some variables used for external communications trough the Py wrapper. ---*/ - PyWrapVarCoord[0] = 0.0; - PyWrapVarCoord[1] = 0.0; - PyWrapVarCoord[2] = 0.0; - PyWrapNodalForce[0] = 0.0; - PyWrapNodalForce[1] = 0.0; - PyWrapNodalForce[2] = 0.0; - PyWrapNodalForceDensity[0] = 0.0; - PyWrapNodalForceDensity[1] = 0.0; - PyWrapNodalForceDensity[2] = 0.0; - PyWrapNodalHeatFlux[0] = 0.0; - PyWrapNodalHeatFlux[1] = 0.0; - PyWrapNodalHeatFlux[2] = 0.0; } +///////////////////////////////////////////////////////////////////////////// +/* Functions related to the global performance indices (Lift, Drag, ecc..) */ +///////////////////////////////////////////////////////////////////////////// + passivedouble CDriver::Get_Drag() { unsigned short val_iZone = ZONE_0; @@ -178,6 +169,10 @@ passivedouble CDriver::Get_LiftCoeff() { return SU2_TYPE::GetValue(CLift); } +///////////////////////////////////////////////////////////////////////////// +/* Functions to obtain information from the geometry/mesh */ +///////////////////////////////////////////////////////////////////////////// + unsigned long CDriver::GetNumberVertices(unsigned short iMarker){ return geometry_container[ZONE_0][INST_0][MESH_0]->nVertex[iMarker]; @@ -219,21 +214,6 @@ bool CDriver::IsAHaloNode(unsigned short iMarker, unsigned long iVertex) { } -unsigned long CDriver::GetnTimeIter() { - - return config_container[ZONE_0]->GetnTime_Iter(); -} - -unsigned long CDriver::GetTime_Iter() const{ - - return TimeIter; -} - -passivedouble CDriver::GetUnsteady_TimeStep(){ - - return SU2_TYPE::GetValue(config_container[ZONE_0]->GetTime_Step()); -} - vector CDriver::GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex) { vector coord(3,0.0); @@ -251,151 +231,50 @@ vector CDriver::GetInitialMeshCoord(unsigned short iMarker, unsig return coord_passive; } -passivedouble CDriver::GetVertexCoordX(unsigned short iMarker, unsigned long iVertex) { - - auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); - return SU2_TYPE::GetValue(Coord[0]); -} - -passivedouble CDriver::GetVertexCoordY(unsigned short iMarker, unsigned long iVertex) { - - auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); - return SU2_TYPE::GetValue(Coord[1]); -} - -passivedouble CDriver::GetVertexCoordZ(unsigned short iMarker, unsigned long iVertex) { - - if(nDim == 2) return 0.0; - - auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); - return SU2_TYPE::GetValue(Coord[2]); -} - -bool CDriver::ComputeVertexForces(unsigned short iMarker, unsigned long iVertex) { - - unsigned long iPoint; - unsigned short iDim; - su2double *Normal, Area; - - unsigned short FinestMesh = config_container[ZONE_0]->GetFinestMesh(); - - /*--- Check the kind of fluid problem ---*/ - bool compressible = (config_container[ZONE_0]->GetKind_Regime() == COMPRESSIBLE); - bool incompressible = (config_container[ZONE_0]->GetKind_Regime() == INCOMPRESSIBLE); - bool viscous_flow = config_container[ZONE_0]->GetViscous(); - - /*--- Parameters for the calculations ---*/ - // Pn: Pressure - // Pinf: Pressure_infinite - su2double Pn = 0.0; - su2double Viscosity = 0.0; - su2double Tau[3][3] = {{0.0}}; - - su2double Pinf = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetPressure_Inf(); - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); +vector CDriver::GetVertexUnitNormal(unsigned short iMarker, unsigned long iVertex){ - /*--- It is necessary to distinguish the halo nodes from the others, since they introduce non physical forces. ---*/ - if (!geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetDomain(iPoint)) - return true; + su2double *Normal; + su2double Area; + vector ret_Normal(3, 0.0); + vector ret_Normal_passive(3, 0.0); - /*--- Get the normal at the vertex: this normal goes inside the fluid domain. ---*/ Normal = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNormal(); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Get the values of pressure and viscosity ---*/ - Pn = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetPressure(iPoint); - if (viscous_flow) { - Viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); - } - - /*--- Calculate the inviscid (pressure) part of tn in the fluid nodes (force units) ---*/ - for (iDim = 0; iDim < nDim; iDim++) { - PyWrapNodalForce[iDim] = -(Pn-Pinf)*Normal[iDim]; //NB : norm(Normal) = Area - } - - /*--- Calculate the viscous (shear stress) part of tn in the fluid nodes (force units) ---*/ - if ((incompressible || compressible) && viscous_flow) { - CNumerics::ComputeStressTensor(nDim, Tau, - solver_container[ZONE_0][INST_0][FinestMesh][FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint)+1, Viscosity); - for (iDim = 0; iDim < nDim; iDim++) { - PyWrapNodalForce[iDim] += GeometryToolbox::DotProduct(nDim, Tau[iDim], Normal); - } - } - - //Divide by local are in case of force density communication. - for(iDim = 0; iDim < nDim; iDim++) - PyWrapNodalForceDensity[iDim] = PyWrapNodalForce[iDim]/Area; - return false; -} - -passivedouble CDriver::GetVertexForceX(unsigned short iMarker, unsigned long iVertex) { - return SU2_TYPE::GetValue(PyWrapNodalForce[0]); -} - -passivedouble CDriver::GetVertexForceY(unsigned short iMarker, unsigned long iVertex) { - return SU2_TYPE::GetValue(PyWrapNodalForce[1]); -} - -passivedouble CDriver::GetVertexForceZ(unsigned short iMarker, unsigned long iVertex) { - return SU2_TYPE::GetValue(PyWrapNodalForce[2]); -} + Area = GeometryToolbox::Norm(nDim, Normal); -passivedouble CDriver::GetVertexForceDensityX(unsigned short iMarker, unsigned long iVertex) { - return SU2_TYPE::GetValue(PyWrapNodalForceDensity[0]); -} + ret_Normal[0] = Normal[0]/Area; + ret_Normal[1] = Normal[1]/Area; + if(nDim>2) ret_Normal[2] = Normal[2]/Area; -passivedouble CDriver::GetVertexForceDensityY(unsigned short iMarker, unsigned long iVertex) { - return SU2_TYPE::GetValue(PyWrapNodalForceDensity[1]); -} + ret_Normal_passive[0] = SU2_TYPE::GetValue(ret_Normal[0]); + ret_Normal_passive[1] = SU2_TYPE::GetValue(ret_Normal[1]); + ret_Normal_passive[2] = SU2_TYPE::GetValue(ret_Normal[2]); -passivedouble CDriver::GetVertexForceDensityZ(unsigned short iMarker, unsigned long iVertex) { - return SU2_TYPE::GetValue(PyWrapNodalForceDensity[2]); + return ret_Normal_passive; } -void CDriver::SetVertexCoordX(unsigned short iMarker, unsigned long iVertex, passivedouble newPosX) { +////////////////////////////////////////////////////////////////////////////////// +/* Functions to obtain global parameters from SU2 (time steps, delta t, ecc...) */ +////////////////////////////////////////////////////////////////////////////////// - auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); +unsigned long CDriver::GetnTimeIter() { - PyWrapVarCoord[0] = newPosX - Coord[0]; + return config_container[ZONE_0]->GetnTime_Iter(); } -void CDriver::SetVertexCoordY(unsigned short iMarker, unsigned long iVertex, passivedouble newPosY) { - - auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); +unsigned long CDriver::GetTime_Iter() const{ - PyWrapVarCoord[1] = newPosY - Coord[1]; + return TimeIter; } -void CDriver::SetVertexCoordZ(unsigned short iMarker, unsigned long iVertex, passivedouble newPosZ) { - - auto iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - auto Coord = geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetCoord(iPoint); +passivedouble CDriver::GetUnsteady_TimeStep(){ - if(nDim > 2) { - PyWrapVarCoord[2] = newPosZ - Coord[2]; - } - else { - PyWrapVarCoord[2] = 0.0; - } + return SU2_TYPE::GetValue(config_container[ZONE_0]->GetTime_Step()); } -passivedouble CDriver::SetVertexVarCoord(unsigned short iMarker, unsigned long iVertex) { - - su2double nodalVarCoordNorm; - - geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->SetVarCoord(PyWrapVarCoord); - nodalVarCoordNorm = sqrt((PyWrapVarCoord[0])*(PyWrapVarCoord[0]) + (PyWrapVarCoord[1])*(PyWrapVarCoord[1]) + (PyWrapVarCoord[2])*(PyWrapVarCoord[2])); - - return SU2_TYPE::GetValue(nodalVarCoordNorm); - -} +/////////////////////////////////////////////////////////////////////////////// +/* Functions related to CHT solver */ +/////////////////////////////////////////////////////////////////////////////// passivedouble CDriver::GetVertexTemperature(unsigned short iMarker, unsigned long iVertex){ @@ -419,7 +298,7 @@ void CDriver::SetVertexTemperature(unsigned short iMarker, unsigned long iVertex geometry_container[ZONE_0][INST_0][MESH_0]->SetCustomBoundaryTemperature(iMarker, iVertex, val_WallTemp); } -bool CDriver::ComputeVertexHeatFluxes(unsigned short iMarker, unsigned long iVertex){ +vector CDriver::GetVertexHeatFluxes(unsigned short iMarker, unsigned long iVertex){ unsigned long iPoint; unsigned short iDim; @@ -429,45 +308,29 @@ bool CDriver::ComputeVertexHeatFluxes(unsigned short iMarker, unsigned long iVer su2double Gamma_Minus_One = Gamma - 1.0; su2double Cp = (Gamma / Gamma_Minus_One) * Gas_Constant; su2double laminar_viscosity, thermal_conductivity; - su2double GradT[3] = {0.0,0.0,0.0}; + vector GradT (3,0.0); + vector HeatFlux (3,0.0); + vector HeatFluxPassive (3,0.0); bool compressible = (config_container[ZONE_0]->GetKind_Regime() == COMPRESSIBLE); bool halo; iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - if(geometry_container[ZONE_0][INST_0][MESH_0]->nodes->GetDomain(iPoint)){ - halo = false; - } - else{ - halo = true; - } - - if(!halo && compressible){ + if(compressible){ laminar_viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); thermal_conductivity = Cp * (laminar_viscosity/Prandtl_Lam); for(iDim=0; iDim < nDim; iDim++){ GradT[iDim] = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint, 0, iDim); - PyWrapNodalHeatFlux[iDim] = -thermal_conductivity*GradT[iDim]; + HeatFlux[iDim] = -thermal_conductivity*GradT[iDim]; } } - return halo; -} + HeatFluxPassive[0] = SU2_TYPE::GetValue(HeatFlux[0]); + HeatFluxPassive[1] = SU2_TYPE::GetValue(HeatFlux[1]); + HeatFluxPassive[2] = SU2_TYPE::GetValue(HeatFlux[2]); -passivedouble CDriver::GetVertexHeatFluxX(unsigned short iMarker, unsigned long iVertex){ - - return SU2_TYPE::GetValue(PyWrapNodalHeatFlux[0]); -} - -passivedouble CDriver::GetVertexHeatFluxY(unsigned short iMarker, unsigned long iVertex){ - - return SU2_TYPE::GetValue(PyWrapNodalHeatFlux[1]); -} - -passivedouble CDriver::GetVertexHeatFluxZ(unsigned short iMarker, unsigned long iVertex){ - - return SU2_TYPE::GetValue(PyWrapNodalHeatFlux[2]); + return HeatFluxPassive; } passivedouble CDriver::GetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVertex){ @@ -536,27 +399,9 @@ passivedouble CDriver::GetThermalConductivity(unsigned short iMarker, unsigned l } -vector CDriver::GetVertexUnitNormal(unsigned short iMarker, unsigned long iVertex){ - - su2double *Normal; - su2double Area; - vector ret_Normal(3, 0.0); - vector ret_Normal_passive(3, 0.0); - - Normal = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNormal(); - - Area = GeometryToolbox::Norm(nDim, Normal); - - ret_Normal[0] = Normal[0]/Area; - ret_Normal[1] = Normal[1]/Area; - if(nDim>2) ret_Normal[2] = Normal[2]/Area; - - ret_Normal_passive[0] = SU2_TYPE::GetValue(ret_Normal[0]); - ret_Normal_passive[1] = SU2_TYPE::GetValue(ret_Normal[1]); - ret_Normal_passive[2] = SU2_TYPE::GetValue(ret_Normal[2]); - - return ret_Normal_passive; -} +//////////////////////////////////////////////////////////////////////////////// +/* Functions related to the management of markers */ +//////////////////////////////////////////////////////////////////////////////// vector CDriver::GetAllBoundaryMarkersTag(){ @@ -575,23 +420,6 @@ vector CDriver::GetAllBoundaryMarkersTag(){ return boundariesTagList; } -vector CDriver::GetAllMovingMarkersTag(){ - - vector movingBoundariesTagList; - unsigned short iMarker, nBoundariesMarker; - string Marker_Tag; - - nBoundariesMarker = config_container[ZONE_0]->GetnMarker_Moving(); - movingBoundariesTagList.resize(nBoundariesMarker); - - for(iMarker=0; iMarker < nBoundariesMarker; iMarker++){ - Marker_Tag = config_container[ZONE_0]->GetMarker_Moving_TagBound(iMarker); - movingBoundariesTagList[iMarker] = Marker_Tag; - } - - return movingBoundariesTagList; -} - vector CDriver::GetAllDeformMeshMarkersTag(){ vector interfaceBoundariesTagList; @@ -609,23 +437,6 @@ vector CDriver::GetAllDeformMeshMarkersTag(){ return interfaceBoundariesTagList; } -vector CDriver::GetAllFluidLoadMarkersTag(){ - - vector interfaceBoundariesTagList; - unsigned short iMarker, nBoundariesMarker; - string Marker_Tag; - - nBoundariesMarker = config_container[ZONE_0]->GetnMarker_Fluid_Load(); - interfaceBoundariesTagList.resize(nBoundariesMarker); - - for(iMarker=0; iMarker < nBoundariesMarker; iMarker++){ - Marker_Tag = config_container[ZONE_0]->GetMarker_Fluid_Load_TagBound(iMarker); - interfaceBoundariesTagList[iMarker] = Marker_Tag; - } - - return interfaceBoundariesTagList; -} - vector CDriver::GetAllCHTMarkersTag(){ vector CHTBoundariesTagList; @@ -633,7 +444,6 @@ vector CDriver::GetAllCHTMarkersTag(){ string Marker_Tag; nBoundariesMarker = config_container[ZONE_0]->GetnMarker_All(); - //CHTBoundariesTagList.resize(nBoundariesMarker); //The CHT markers can be identified as the markers that are customizable with a BC type HEAT_FLUX or ISOTHERMAL. for(iMarker=0; iMarker CDriver::GetAllBoundaryMarkersType(){ return allBoundariesTypeMap; } +void CDriver::SetHeatSource_Position(passivedouble alpha, passivedouble pos_x, passivedouble pos_y, passivedouble pos_z){ + + CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][RAD_SOL]; + + config_container[ZONE_0]->SetHeatSource_Rot_Z(alpha); + config_container[ZONE_0]->SetHeatSource_Center(pos_x, pos_y, pos_z); + + solver->SetVolumetricHeatSource(geometry_container[ZONE_0][INST_0][MESH_0], config_container[ZONE_0]); + +} + +void CDriver::SetInlet_Angle(unsigned short iMarker, passivedouble alpha){ + + su2double alpha_rad = alpha * PI_NUMBER/180.0; + + unsigned long iVertex; + + for (iVertex = 0; iVertex < geometry_container[ZONE_0][INST_0][MESH_0]->nVertex[iMarker]; iVertex++){ + solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, 0, cos(alpha_rad)); + solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, 1, sin(alpha_rad)); + } + +} + +///////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/* Functions related to simulation control, high level functions (reset convergence, set initial mesh, ecc...) */ +///////////////////////////////////////////////////////////////////////////////////////////////////////////////// + void CDriver::ResetConvergence() { for(iZone = 0; iZone < nZone; iZone++) { @@ -778,45 +616,7 @@ void CSinglezoneDriver::SetInitialMesh() { } } -void CFluidDriver::SetVertexTtotal(unsigned short iMarker, unsigned long iVertex, passivedouble val_Ttotal_passive){ - - su2double val_Ttotal = val_Ttotal_passive; - - solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SetInlet_Ttotal(iMarker, iVertex, val_Ttotal); - -} - -void CFluidDriver::SetVertexPtotal(unsigned short iMarker, unsigned long iVertex, passivedouble val_Ptotal_passive){ - - su2double val_Ptotal = val_Ptotal_passive; - - solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SetInlet_Ptotal(iMarker, iVertex, val_Ptotal); - -} - -void CFluidDriver::SetVertexFlowDir(unsigned short iMarker, unsigned long iVertex, unsigned short iDim, passivedouble val_FlowDir_passive){ - - su2double val_FlowDir = val_FlowDir_passive; - - solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, iDim, val_FlowDir); - -} - -void CFluidDriver::SetVertexTurbVar(unsigned short iMarker, unsigned long iVertex, unsigned short iDim, passivedouble val_turb_var_passive){ - - su2double val_turb_var = val_turb_var_passive; - - if (solver_container[ZONE_0][INST_0] == nullptr || - solver_container[ZONE_0][INST_0][MESH_0] == nullptr) { - SU2_MPI::Error("Could not find an appropriate solver.", CURRENT_FUNCTION); - } else if (solver_container[ZONE_0][INST_0][MESH_0][TURB_SOL] == nullptr) { - SU2_MPI::Error("Tried to set turbulence variables without a turbulence solver.", CURRENT_FUNCTION); - } - solver_container[ZONE_0][INST_0][MESH_0][TURB_SOL]->SetInlet_TurbVar(iMarker, iVertex, iDim, val_turb_var); - -} - -void CFluidDriver::BoundaryConditionsUpdate(){ +void CDriver::BoundaryConditionsUpdate(){ int rank = MASTER_NODE; unsigned short iZone; @@ -829,63 +629,21 @@ void CFluidDriver::BoundaryConditionsUpdate(){ } } -void CDriver::SetMeshDisplacement(unsigned short iMarker, unsigned long iVertex, passivedouble DispX, passivedouble DispY, passivedouble DispZ) { - - unsigned long iPoint; - PyWrapVarCoord[0] = DispX; - PyWrapVarCoord[1] = DispY; - PyWrapVarCoord[2] = DispZ; - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - - solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->SetBound_Disp(iPoint,PyWrapVarCoord); - -} - -void CDriver::CommunicateMeshDisplacement(void) { - - solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->InitiateComms(geometry_container[ZONE_0][INST_0][MESH_0], - config_container[ZONE_0], MESH_DISPLACEMENTS); - solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->CompleteComms(geometry_container[ZONE_0][INST_0][MESH_0], - config_container[ZONE_0], MESH_DISPLACEMENTS); - -} - -vector CDriver::GetMeshDisp_Sensitivity(unsigned short iMarker, unsigned long iVertex) { - - unsigned long iPoint; - vector Disp_Sens(3, 0.0); - vector Disp_Sens_passive(3, 0.0); - - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][ADJMESH_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; - - Disp_Sens[0] = solver->GetNodes()->GetBoundDisp_Sens(iPoint, 0); - Disp_Sens[1] = solver->GetNodes()->GetBoundDisp_Sens(iPoint, 1); - if (geometry->GetnDim() == 3) - Disp_Sens[2] = solver->GetNodes()->GetBoundDisp_Sens(iPoint, 2); - else - Disp_Sens[2] = 0.0; - - Disp_Sens_passive[0] = SU2_TYPE::GetValue(Disp_Sens[0]); - Disp_Sens_passive[1] = SU2_TYPE::GetValue(Disp_Sens[1]); - Disp_Sens_passive[2] = SU2_TYPE::GetValue(Disp_Sens[2]); - - return Disp_Sens_passive; - -} +//////////////////////////////////////////////////////////////////////////////// +/* Functions related to finite elements */ +//////////////////////////////////////////////////////////////////////////////// void CDriver::SetFEA_Loads(unsigned short iMarker, unsigned long iVertex, passivedouble LoadX, passivedouble LoadY, passivedouble LoadZ) { unsigned long iPoint; - PyWrapNodalForce[0] = LoadX; - PyWrapNodalForce[1] = LoadY; - PyWrapNodalForce[2] = LoadZ; + vector NodalForce (3,0.0); + NodalForce[0] = LoadX; + NodalForce[1] = LoadY; + NodalForce[2] = LoadZ; iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]->GetNodes()->Set_FlowTraction(iPoint,PyWrapNodalForce); + solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]->GetNodes()->Set_FlowTraction(iPoint,NodalForce); } @@ -967,6 +725,37 @@ vector CDriver::GetFEA_Velocity_n(unsigned short iMarker, unsigne } +//////////////////////////////////////////////////////////////////////////////// +/* Functions related to adjoint simulations */ +//////////////////////////////////////////////////////////////////////////////// + +vector CDriver::GetMeshDisp_Sensitivity(unsigned short iMarker, unsigned long iVertex) { + + unsigned long iPoint; + vector Disp_Sens(3, 0.0); + vector Disp_Sens_passive(3, 0.0); + + iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); + CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][ADJMESH_SOL]; + CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + + Disp_Sens[0] = solver->GetNodes()->GetBoundDisp_Sens(iPoint, 0); + Disp_Sens[1] = solver->GetNodes()->GetBoundDisp_Sens(iPoint, 1); + if (geometry->GetnDim() == 3) + Disp_Sens[2] = solver->GetNodes()->GetBoundDisp_Sens(iPoint, 2); + else + Disp_Sens[2] = 0.0; + + Disp_Sens_passive[0] = SU2_TYPE::GetValue(Disp_Sens[0]); + Disp_Sens_passive[1] = SU2_TYPE::GetValue(Disp_Sens[1]); + Disp_Sens_passive[2] = SU2_TYPE::GetValue(Disp_Sens[2]); + + return Disp_Sens_passive; + +} + + + vector CDriver::GetFlowLoad_Sensitivity(unsigned short iMarker, unsigned long iVertex) { unsigned long iPoint; @@ -992,31 +781,6 @@ vector CDriver::GetFlowLoad_Sensitivity(unsigned short iMarker, u } -vector CDriver::GetFlowLoad(unsigned short iMarker, unsigned long iVertex) { - - vector FlowLoad(3, 0.0); - vector FlowLoad_passive(3, 0.0); - - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; - - if (config_container[ZONE_0]->GetSolid_Wall(iMarker)) { - FlowLoad[0] = solver->GetVertexTractions(iMarker, iVertex, 0); - FlowLoad[1] = solver->GetVertexTractions(iMarker, iVertex, 1); - if (geometry->GetnDim() == 3) - FlowLoad[2] = solver->GetVertexTractions(iMarker, iVertex, 2); - else - FlowLoad[2] = 0.0; - } - - FlowLoad_passive[0] = SU2_TYPE::GetValue(FlowLoad[0]); - FlowLoad_passive[1] = SU2_TYPE::GetValue(FlowLoad[1]); - FlowLoad_passive[2] = SU2_TYPE::GetValue(FlowLoad[2]); - - return FlowLoad_passive; - -} - void CDriver::SetFlowLoad_Adjoint(unsigned short iMarker, unsigned long iVertex, passivedouble val_AdjointX, passivedouble val_AdjointY, passivedouble val_AdjointZ) { @@ -1046,53 +810,59 @@ void CDriver::SetSourceTerm_DispAdjoint(unsigned short iMarker, unsigned long iV } -vector CDriver::GetVertex_UndeformedCoord(unsigned short iMarker, unsigned long iVertex) { +//////////////////////////////////////////////////////////////////////////////// +/* Functions related to mesh deformation */ +//////////////////////////////////////////////////////////////////////////////// - unsigned long iPoint; - vector MeshCoord(3, 0.0); - vector MeshCoord_passive(3, 0.0); +void CDriver::SetMeshDisplacement(unsigned short iMarker, unsigned long iVertex, passivedouble DispX, passivedouble DispY, passivedouble DispZ) { - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]; - CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; - iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); + unsigned long iPoint; + vector MeshDispl (3,0.0); - if (solver != nullptr) { - MeshCoord[0] = solver->GetNodes()->GetMesh_Coord(iPoint,0); - MeshCoord[1] = solver->GetNodes()->GetMesh_Coord(iPoint,1); - if (geometry->GetnDim() == 3) - MeshCoord[2] = solver->GetNodes()->GetMesh_Coord(iPoint,2); - else - MeshCoord[2] = 0.0; - } + MeshDispl[0] = DispX; + MeshDispl[1] = DispY; + MeshDispl[2] = DispZ; - MeshCoord_passive[0] = SU2_TYPE::GetValue(MeshCoord[0]); - MeshCoord_passive[1] = SU2_TYPE::GetValue(MeshCoord[1]); - MeshCoord_passive[2] = SU2_TYPE::GetValue(MeshCoord[2]); + iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - return MeshCoord_passive; + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->SetBound_Disp(iPoint,MeshDispl); } -void CDriver::SetHeatSource_Position(passivedouble alpha, passivedouble pos_x, passivedouble pos_y, passivedouble pos_z){ - - CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][RAD_SOL]; - - config_container[ZONE_0]->SetHeatSource_Rot_Z(alpha); - config_container[ZONE_0]->SetHeatSource_Center(pos_x, pos_y, pos_z); +void CDriver::CommunicateMeshDisplacement(void) { - solver->SetVolumetricHeatSource(geometry_container[ZONE_0][INST_0][MESH_0], config_container[ZONE_0]); + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->InitiateComms(geometry_container[ZONE_0][INST_0][MESH_0], + config_container[ZONE_0], MESH_DISPLACEMENTS); + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->CompleteComms(geometry_container[ZONE_0][INST_0][MESH_0], + config_container[ZONE_0], MESH_DISPLACEMENTS); } -void CDriver::SetInlet_Angle(unsigned short iMarker, passivedouble alpha){ +//////////////////////////////////////////////////////////////////////////////// +/* Functions related to flow loads */ +//////////////////////////////////////////////////////////////////////////////// - su2double alpha_rad = alpha * PI_NUMBER/180.0; +vector CDriver::GetFlowLoad(unsigned short iMarker, unsigned long iVertex) { - unsigned long iVertex; + vector FlowLoad(3, 0.0); + vector FlowLoad_passive(3, 0.0); - for (iVertex = 0; iVertex < geometry_container[ZONE_0][INST_0][MESH_0]->nVertex[iMarker]; iVertex++){ - solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, 0, cos(alpha_rad)); - solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, 1, sin(alpha_rad)); + CSolver *solver = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]; + CGeometry *geometry = geometry_container[ZONE_0][INST_0][MESH_0]; + + if (config_container[ZONE_0]->GetSolid_Wall(iMarker)) { + FlowLoad[0] = solver->GetVertexTractions(iMarker, iVertex, 0); + FlowLoad[1] = solver->GetVertexTractions(iMarker, iVertex, 1); + if (geometry->GetnDim() == 3) + FlowLoad[2] = solver->GetVertexTractions(iMarker, iVertex, 2); + else + FlowLoad[2] = 0.0; } + FlowLoad_passive[0] = SU2_TYPE::GetValue(FlowLoad[0]); + FlowLoad_passive[1] = SU2_TYPE::GetValue(FlowLoad[1]); + FlowLoad_passive[2] = SU2_TYPE::GetValue(FlowLoad[2]); + + return FlowLoad_passive; + } diff --git a/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py b/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py index 1ed717dd6974..06a3edd2b7eb 100755 --- a/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py +++ b/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py @@ -45,7 +45,7 @@ from math import * # ------------------------------------------------------------------- -# Main +# Main # ------------------------------------------------------------------- def main(): @@ -53,38 +53,12 @@ def main(): # Command line options parser=OptionParser() parser.add_option("-f", "--file", dest="filename", help="Read config from FILE", metavar="FILE") - parser.add_option("--nDim", dest="nDim", default=2, help="Define the number of DIMENSIONS", - metavar="DIMENSIONS") - parser.add_option("--nZone", dest="nZone", default=1, help="Define the number of ZONES", - metavar="ZONES") parser.add_option("--parallel", action="store_true", help="Specify if we need to initialize MPI", dest="with_MPI", default=False) - parser.add_option("--fsi", dest="fsi", default="False", help="Launch the FSI driver", metavar="FSI") - - parser.add_option("--fem", dest="fem", default="False", help="Launch the FEM driver (General driver)", metavar="FEM") - - parser.add_option("--harmonic_balance", dest="harmonic_balance", default="False", - help="Launch the Harmonic Balance (HB) driver", metavar="HB") - - parser.add_option("--poisson_equation", dest="poisson_equation", default="False", - help="Launch the poisson equation driver (General driver)", metavar="POIS_EQ") - - parser.add_option("--wave_equation", dest="wave_equation", default="False", - help="Launch the wave equation driver (General driver)", metavar="WAVE_EQ") - - parser.add_option("--heat_equation", dest="heat_equation", default="False", - help="Launch the heat equation driver (General driver)", metavar="HEAT_EQ") - (options, args) = parser.parse_args() - options.nDim = int( options.nDim ) - options.nZone = int( options.nZone ) - options.fsi = options.fsi.upper() == 'TRUE' - options.fem = options.fem.upper() == 'TRUE' - options.harmonic_balance = options.harmonic_balance.upper() == 'TRUE' - options.poisson_equation = options.poisson_equation.upper() == 'TRUE' - options.wave_equation = options.wave_equation.upper() == 'TRUE' - options.heat_equation = options.heat_equation.upper() == 'TRUE' + options.nDim = int(2) + options.nZone = int(1) # Import mpi4py for parallel run if options.with_MPI == True: @@ -92,18 +66,11 @@ def main(): comm = MPI.COMM_WORLD rank = comm.Get_rank() else: - comm = 0 + comm = 0 rank = 0 # Initialize the corresponding driver of SU2, this includes solver preprocessing try: - if (options.nZone == 1) and ( options.fem or options.poisson_equation or options.wave_equation or options.heat_equation ): - SU2Driver = pysu2.CGeneralDriver(options.filename, options.nZone, comm); - elif options.harmonic_balance: - SU2Driver = pysu2.CHBDriver(options.filename, options.nZone, comm); - elif (options.nZone == 2) and (options.fsi): - SU2Driver = pysu2.CFSIDriver(options.filename, options.nZone, comm); - else: SU2Driver = pysu2.CSinglezoneDriver(options.filename, options.nZone, comm); except TypeError as exception: print('A TypeError occured in pysu2.CDriver : ',exception) @@ -162,6 +129,8 @@ def main(): SU2Driver.BoundaryConditionsUpdate() # Run one time iteration (e.g. dual-time) SU2Driver.Run() + # Postprocess the solver and exit cleanly + SU2Driver.Postprocess() # Update the solver for the next time iteration SU2Driver.Update() # Monitor the solver and output solution to file if required @@ -173,9 +142,6 @@ def main(): TimeIter += 1 time += deltaT - # Postprocess the solver and exit cleanly - SU2Driver.Postprocessing() - if SU2Driver != None: del SU2Driver @@ -185,4 +151,4 @@ def main(): # this is only accessed if running from command prompt if __name__ == '__main__': - main() + main() From 96e23ea2a9a29046344c924ec10b469fe7ef404f Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 15 Feb 2021 11:06:13 +0100 Subject: [PATCH 23/32] Fixes to vector data --- SU2_CFD/src/python_wrapper_structure.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 383ddc2eec9d..14d9b5bc362e 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -643,7 +643,7 @@ void CDriver::SetFEA_Loads(unsigned short iMarker, unsigned long iVertex, passiv NodalForce[2] = LoadZ; iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]->GetNodes()->Set_FlowTraction(iPoint,NodalForce); + solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]->GetNodes()->Set_FlowTraction(iPoint,NodalForce.data()); } @@ -825,7 +825,7 @@ void CDriver::SetMeshDisplacement(unsigned short iMarker, unsigned long iVertex, iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->SetBound_Disp(iPoint,MeshDispl); + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->SetBound_Disp(iPoint,MeshDispl.data()); } From fbcca5870e04284026bddb1ad2010e27755738b9 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 15 Feb 2021 11:29:22 +0100 Subject: [PATCH 24/32] Better initialisation of variables --- SU2_CFD/src/python_wrapper_structure.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 14d9b5bc362e..80b4d9aa98e2 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -637,13 +637,13 @@ void CDriver::SetFEA_Loads(unsigned short iMarker, unsigned long iVertex, passiv passivedouble LoadY, passivedouble LoadZ) { unsigned long iPoint; - vector NodalForce (3,0.0); + su2double NodalForce[3] = {0.0,0.0,0.0}; NodalForce[0] = LoadX; NodalForce[1] = LoadY; NodalForce[2] = LoadZ; iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]->GetNodes()->Set_FlowTraction(iPoint,NodalForce.data()); + solver_container[ZONE_0][INST_0][MESH_0][FEA_SOL]->GetNodes()->Set_FlowTraction(iPoint,NodalForce); } @@ -817,7 +817,7 @@ void CDriver::SetSourceTerm_DispAdjoint(unsigned short iMarker, unsigned long iV void CDriver::SetMeshDisplacement(unsigned short iMarker, unsigned long iVertex, passivedouble DispX, passivedouble DispY, passivedouble DispZ) { unsigned long iPoint; - vector MeshDispl (3,0.0); + su2double MeshDispl[3] = {0.0,0.0,0.0}; MeshDispl[0] = DispX; MeshDispl[1] = DispY; @@ -825,7 +825,7 @@ void CDriver::SetMeshDisplacement(unsigned short iMarker, unsigned long iVertex, iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); - solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->SetBound_Disp(iPoint,MeshDispl.data()); + solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->SetBound_Disp(iPoint,MeshDispl); } From 8a14d4df0a386944d06ec498494179b0e52501a7 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 15 Feb 2021 11:37:39 +0100 Subject: [PATCH 25/32] Small typo in parallel regression --- TestCases/parallel_regression.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 98d91e311d3d..ab0f720073f2 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -59,7 +59,7 @@ def main(): thermalbath_frozen.cfg_dir = "nonequilibrium/thermalbath/frozen" thermalbath_frozen.cfg_file = "thermalbath_frozen.cfg" thermalbath_frozen.test_iter = 10 - thermalbath_frozen.test_vals = [ -32.000000, -32.000000, -11.92359, -11.962329, -32.000000, 10.813864] + thermalbath_frozen.test_vals = [ -32.000000, -32.000000, -11.92359, -11.962329, -32.000000, 10.813864] thermalbath_frozen.su2_exec = "mpirun -n 2 SU2_CFD" thermalbath_frozen.timeout = 1600 thermalbath_frozen.new_output = True @@ -1298,7 +1298,7 @@ def main(): pywrapper_rigidMotion.cfg_dir = "py_wrapper/flatPlate_rigidMotion" pywrapper_rigidMotion.cfg_file = "flatPlate_rigidMotion_Conf.cfg" pywrapper_rigidMotion.test_iter = 5 - pywrapper_rigidMotion.test_vals = [-1.614170, 2.242953, 0.350050, 0.093137] + pywrapper_rigidMotion.test_vals = [-1.614170, 2.242953, 0.350036, 0.093137] pywrapper_rigidMotion.su2_exec = "mpirun -np 2 python launch_flatPlate_rigidMotion.py --parallel -f" pywrapper_rigidMotion.timeout = 1600 pywrapper_rigidMotion.tol = 0.00001 From 3aaf175504ab56489dcb0a10d60c9c2937128fb4 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Mon, 15 Feb 2021 17:06:28 +0100 Subject: [PATCH 26/32] Fixed regression values after PR#1199 --- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index ab0f720073f2..b5c198277f5d 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1298,7 +1298,7 @@ def main(): pywrapper_rigidMotion.cfg_dir = "py_wrapper/flatPlate_rigidMotion" pywrapper_rigidMotion.cfg_file = "flatPlate_rigidMotion_Conf.cfg" pywrapper_rigidMotion.test_iter = 5 - pywrapper_rigidMotion.test_vals = [-1.614170, 2.242953, 0.350036, 0.093137] + pywrapper_rigidMotion.test_vals = [-1.551335, 2.295594, 0.350036, 0.093081] pywrapper_rigidMotion.su2_exec = "mpirun -np 2 python launch_flatPlate_rigidMotion.py --parallel -f" pywrapper_rigidMotion.timeout = 1600 pywrapper_rigidMotion.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 77dccd00fe71..daeb8e95cb51 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1877,7 +1877,7 @@ def main(): pywrapper_rigidMotion.cfg_dir = "py_wrapper/flatPlate_rigidMotion" pywrapper_rigidMotion.cfg_file = "flatPlate_rigidMotion_Conf.cfg" pywrapper_rigidMotion.test_iter = 5 - pywrapper_rigidMotion.test_vals = [-1.614170, 2.242953, 0.350050, 0.093137] + pywrapper_rigidMotion.test_vals = [-1.551335, 2.295594, 0.350050, 0.093081] pywrapper_rigidMotion.su2_exec = "python launch_flatPlate_rigidMotion.py -f" pywrapper_rigidMotion.new_output = True pywrapper_rigidMotion.timeout = 1600 From 52ba9900bad912fe4590f30748feafd1c5779d1e Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Tue, 16 Feb 2021 09:35:35 +0100 Subject: [PATCH 27/32] array deleted twice --- SU2_PY/FSI_tools/FSIInterface.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/SU2_PY/FSI_tools/FSIInterface.py b/SU2_PY/FSI_tools/FSIInterface.py index 702ad4d78b0a..9fe034d1e1fa 100644 --- a/SU2_PY/FSI_tools/FSIInterface.py +++ b/SU2_PY/FSI_tools/FSIInterface.py @@ -1392,9 +1392,7 @@ def interpolateFluidLoadsOnSolidMesh(self, FSI_config): del sendBuff_X del sendBuff_Y del sendBuff_Z - del self.solidLoads_array_X_recon - del self.solidLoads_array_Y_recon - del self.solidLoads_array_Z_recon + self.comm.barrier() else: self.localSolidLoads_array_X = self.solidLoads_array_X.getArray().copy() self.localSolidLoads_array_Y = self.solidLoads_array_Y.getArray().copy() From 2f1e61dfbd45e320082b4d4c8969e92a6ef76e4b Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Tue, 16 Feb 2021 14:10:01 +0100 Subject: [PATCH 28/32] Multiple superposed forced motions --- SU2_PY/SU2_Nastran/pysu2_nastran.py | 74 +++++++++++++++++------------ 1 file changed, 43 insertions(+), 31 deletions(-) diff --git a/SU2_PY/SU2_Nastran/pysu2_nastran.py b/SU2_PY/SU2_Nastran/pysu2_nastran.py index 6f6144c4bc28..437453fb136a 100644 --- a/SU2_PY/SU2_Nastran/pysu2_nastran.py +++ b/SU2_PY/SU2_Nastran/pysu2_nastran.py @@ -43,21 +43,24 @@ class ImposedMotionFunction: - def __init__(self,time0,tipo,parameters): + def __init__(self,time0,typeOfMotion,parameters,mode): + self.time0 = time0 - self.tipo = tipo - if self.tipo == "SINUSOIDAL": - self.bias = parameters[0] - self.amplitude = parameters[1] - self.frequency = parameters[2] - self.timeStart = parameters[3] - - elif self.tipo == "BLENDED_STEP": - self.kmax = parameters[0] - self.vinf = parameters[1] - self.lref = parameters[2] - self.amplitude = parameters[3] - self.timeStart = parameters[4] + self.typeOfMotion = typeOfMotion + self.mode = mode + + if self.typeOfMotion == "SINUSOIDAL": + self.bias = parameters["BIAS"] + self.amplitude = parameters["AMPLITUDE"] + self.frequency = parameters["FREQUENCY"] + self.timeStart = parameters["TIME_0"] + + elif self.typeOfMotion == "BLENDED_STEP": + self.kmax = parameters["K_MAX"] + self.vinf = parameters["V_INF"] + self.lref = parameters["L_REF"] + self.amplitude = parameters["AMPLITUDE"] + self.timeStart = parameters["TIME_0"] self.tmax = 2*pi/self.kmax*self.lref/self.vinf self.omega0 = 1/2*self.kmax @@ -67,10 +70,10 @@ def __init__(self,time0,tipo,parameters): def GetDispl(self,time): time = time - self.time0 - self.timeStart - if self.tipo == "SINUSOIDAL": + if self.typeOfMotion == "SINUSOIDAL": return self.bias+self.amplitude*sin(2*pi*self.frequency*time) - if self.tipo == "BLENDED_STEP": + if self.typeOfMotion == "BLENDED_STEP": if time < 0: return 0.0 elif time < self.tmax: @@ -81,10 +84,10 @@ def GetDispl(self,time): def GetVel(self,time): time = time - self.time0 - self.timeStart - if self.tipo == "SINUSOIDAL": + if self.typeOfMotion == "SINUSOIDAL": return self.amplitude*cos(2*pi*self.frequency*time)*2*pi*self.frequency - if self.tipo == "BLENDED_STEP": + if self.typeOfMotion == "BLENDED_STEP": if time < 0: return 0.0 elif time < self.tmax: @@ -94,10 +97,10 @@ def GetVel(self,time): def GetAcc(self,time): time = time - self.time0 - self.timeStart - if self.tipo == "SINUSOIDAL": + if self.typeOfMotion == "SINUSOIDAL": return -self.amplitude*sin(2*pi*self.frequency*time)*(2*pi*self.frequency)**2 - if self.tipo == "BLENDED_STEP": + if self.typeOfMotion == "BLENDED_STEP": if time < 0: return 0.0 elif time < self.tmax: @@ -289,7 +292,7 @@ def __init__(self, config_fileName, ImposedMotion): self.markers = {} self.refsystems = [] self.ImposedMotionToSet = True - self.ImposedMotionFunction = {} + self.ImposedMotionFunction = [] print("\n") print(" Reading the mesh ".center(80,"-")) @@ -733,14 +736,17 @@ def __temporalIteration(self,time): This method integrates in time the solution. """ + self.__reset(self.q) + self.__reset(self.qdot) + self.__reset(self.qddot) + self.__reset(self.a) + if not self.ImposedMotion: eps = 1e-6 self.__SetLoads() # Prediction step - self.__reset(self.qddot) - self.__reset(self.a) self.a += (self.alpha_f)/(1-self.alpha_m)*self.qddot_n self.a -= (self.alpha_m)/(1-self.alpha_m)*self.a_n @@ -768,14 +774,20 @@ def __temporalIteration(self,time): self.a += (1-self.alpha_f)/(1-self.alpha_m)*self.qddot else: if self.ImposedMotionToSet: - for imode in self.Config["IMPOSED_MODES"].keys(): - self.ImposedMotionFunction[imode] = ImposedMotionFunction(time,self.Config["IMPOSED_MODES"][imode],self.Config["IMPOSED_PARAMETERS"][imode]) - self.ImposedMotionToSet = False - for imode in self.Config["IMPOSED_MODES"].keys(): - self.q[imode] = self.ImposedMotionFunction[imode].GetDispl(time) - self.qdot[imode] = self.ImposedMotionFunction[imode].GetVel(time) - self.qddot[imode] = self.ImposedMotionFunction[imode].GetAcc(time) - self.a = np.copy(self.qddot) + iImposedFunc = 0 + for imode in self.Config["IMPOSED_MODES"].keys(): + for isuperposed in range(len(self.Config["IMPOSED_MODES"][imode])): + typeOfMotion = self.Config["IMPOSED_MODES"][imode][isuperposed] + parameters = self.Config["IMPOSED_PARAMETERS"][imode][isuperposed] + self.ImposedMotionFunction[iImposedFunc] = ImposedMotionFunction(time, typeOfMotion, parameters, imode) + iImposedFunc += 1 + self.ImposedMotionToSet = False + for iImposedFunc in range(len(self.ImposedMotionFunction)): + imode = self.ImposedMotionFunction[iImposedFunc].mode + self.q[imode] += self.ImposedMotionFunction[iImposedFunc].GetDispl(time) + self.qdot[imode] += self.ImposedMotionFunction[iImposedFunc].GetVel(time) + self.qddot[imode] += self.ImposedMotionFunction[iImposedFunc].GetAcc(time) + self.a = np.copy(self.qddot) def __SetLoads(self): From ae338089ba31a2118c6535f629eadaa52fd79607 Mon Sep 17 00:00:00 2001 From: Nicola-Fonzi Date: Tue, 16 Feb 2021 14:45:19 +0100 Subject: [PATCH 29/32] Error with array dimension --- SU2_PY/SU2_Nastran/pysu2_nastran.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_PY/SU2_Nastran/pysu2_nastran.py b/SU2_PY/SU2_Nastran/pysu2_nastran.py index 437453fb136a..aa83f44c5dd2 100644 --- a/SU2_PY/SU2_Nastran/pysu2_nastran.py +++ b/SU2_PY/SU2_Nastran/pysu2_nastran.py @@ -41,7 +41,7 @@ # Config class # ---------------------------------------------------------------------- -class ImposedMotionFunction: +class ImposedMotionClass: def __init__(self,time0,typeOfMotion,parameters,mode): @@ -779,7 +779,7 @@ def __temporalIteration(self,time): for isuperposed in range(len(self.Config["IMPOSED_MODES"][imode])): typeOfMotion = self.Config["IMPOSED_MODES"][imode][isuperposed] parameters = self.Config["IMPOSED_PARAMETERS"][imode][isuperposed] - self.ImposedMotionFunction[iImposedFunc] = ImposedMotionFunction(time, typeOfMotion, parameters, imode) + self.ImposedMotionFunction.append(ImposedMotionClass(time, typeOfMotion, parameters, imode)) iImposedFunc += 1 self.ImposedMotionToSet = False for iImposedFunc in range(len(self.ImposedMotionFunction)): From 9d78c06d1eff53ef53b5c6ec26ab549a63aef6ce Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 19 Feb 2021 14:28:14 +0000 Subject: [PATCH 30/32] constification, remove legacy python FSI --- Common/include/CConfig.hpp | 4 +- Common/include/option_structure.hpp | 4 -- Common/src/CConfig.cpp | 15 +----- SU2_CFD/include/drivers/CDriver.hpp | 47 ++++++++-------- SU2_CFD/src/iteration/CIteration.cpp | 31 ----------- SU2_CFD/src/iteration/CIterationFactory.cpp | 2 +- SU2_CFD/src/python_wrapper_structure.cpp | 54 +++++++++---------- .../cont_adj_rans/oneram6/turb_ONERAM6.cfg | 7 +-- TestCases/gust/inv_gust_NACA0012.cfg | 2 +- TestCases/moving_wall/cavity/lam_cavity.cfg | 3 +- .../spinning_cylinder/spinning_cylinder.cfg | 2 +- .../pitching_NACA64A010.cfg | 2 +- .../pitching_oneram6/pitching_ONERAM6.cfg | 2 +- .../rotating_naca0012/rotating_NACA0012.cfg | 2 +- 14 files changed, 58 insertions(+), 119 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 3aa76c0faa40..91684f20d32d 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -657,7 +657,7 @@ class CConfig { nMarker_ZoneInterface, /*!< \brief Number of markers in the zone interface. */ nMarker_Plotting, /*!< \brief Number of markers to plot. */ nMarker_Analyze, /*!< \brief Number of markers to analyze. */ - nMarker_Moving, /*!< \brief Number of markers in motion (DEFORMING, MOVING_WALL, or FLUID_STRUCTURE). */ + nMarker_Moving, /*!< \brief Number of markers in motion (DEFORMING, MOVING_WALL). */ nMarker_PyCustom, /*!< \brief Number of markers that are customizable in Python. */ nMarker_DV, /*!< \brief Number of markers affected by the design variables. */ nMarker_WallFunctions; /*!< \brief Number of markers for which wall functions must be applied. */ @@ -667,7 +667,7 @@ class CConfig { *Marker_Plotting, /*!< \brief Markers to plot. */ *Marker_Analyze, /*!< \brief Markers to analyze. */ *Marker_ZoneInterface, /*!< \brief Markers in the FSI interface. */ - *Marker_Moving, /*!< \brief Markers in motion (DEFORMING, MOVING_WALL, or FLUID_STRUCTURE). */ + *Marker_Moving, /*!< \brief Markers in motion (DEFORMING, MOVING_WALL). */ *Marker_PyCustom, /*!< \brief Markers that are customizable in Python. */ *Marker_DV, /*!< \brief Markers affected by the design variables. */ *Marker_WallFunctions; /*!< \brief Markers for which wall functions must be applied. */ diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 5a3c8213c379..c7044f8b75ea 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -176,7 +176,6 @@ enum ENUM_MAIN_SOLVER { INC_NAVIER_STOKES =5, /*!< \brief Definition of the incompressible Navier-Stokes' solver. */ INC_RANS = 6, /*!< \brief Definition of the incompressible Reynolds-averaged Navier-Stokes' (RANS) solver. */ HEAT_EQUATION = 7, /*!< \brief Definition of the finite volume heat solver. */ - FLUID_STRUCTURE_INTERACTION = 8, /*!< \brief Definition of a FSI solver. */ FEM_ELASTICITY = 9, /*!< \brief Definition of a FEM solver. */ ADJ_EULER = 10, /*!< \brief Definition of the continuous adjoint Euler's solver. */ ADJ_NAVIER_STOKES = 11, /*!< \brief Definition of the continuous adjoint Navier-Stokes' solver. */ @@ -231,7 +230,6 @@ static const MapType Solver_Map = { MakePair("DISC_ADJ_FEM_RANS", DISC_ADJ_FEM_RANS) MakePair("DISC_ADJ_FEM_NS", DISC_ADJ_FEM_NS) MakePair("DISC_ADJ_FEM", DISC_ADJ_FEM) - MakePair("FLUID_STRUCTURE_INTERACTION", FLUID_STRUCTURE_INTERACTION) MakePair("TEMPLATE_SOLVER", TEMPLATE_SOLVER) MakePair("MULTIPHYSICS", MULTIPHYSICS) }; @@ -699,7 +697,6 @@ enum ENUM_SURFACEMOVEMENT { MOVING_WALL = 2, /*!< \brief Simulation with moving wall. */ AEROELASTIC = 3, /*!< \brief Simulation with aeroelastic motion. */ AEROELASTIC_RIGID_MOTION = 4, /*!< \brief Simulation with rotation and aeroelastic motion. */ - FLUID_STRUCTURE = 5, /*!< \brief Fluid structure deformation. */ EXTERNAL = 6, /*!< \brief Simulation with external motion. */ EXTERNAL_ROTATION = 7, /*!< \brief Simulation with external rotation motion. */ }; @@ -708,7 +705,6 @@ static const MapType SurfaceMovement_Map = { MakePair("MOVING_WALL", MOVING_WALL) MakePair("AEROELASTIC_RIGID_MOTION", AEROELASTIC_RIGID_MOTION) MakePair("AEROELASTIC", AEROELASTIC) - MakePair("FLUID_STRUCTURE", FLUID_STRUCTURE) MakePair("EXTERNAL", EXTERNAL) MakePair("EXTERNAL_ROTATION", EXTERNAL_ROTATION) }; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 01137a94aae3..b434a2f6edd2 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -3114,13 +3114,6 @@ void CConfig::SetnZone(){ } - /*--- Temporary fix until Multizone Disc. Adj. solver is ready ---- */ - - if (Kind_Solver == FLUID_STRUCTURE_INTERACTION){ - - nZone = GetnZone(Mesh_FileName, Mesh_FileFormat); - - } } void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_izone, unsigned short val_nDim) { @@ -3563,11 +3556,7 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ } } - if ((nKind_SurfaceMovement > 1) && GetSurface_Movement(FLUID_STRUCTURE)) { - SU2_MPI::Error("FSI in combination with moving surfaces is currently not supported.", CURRENT_FUNCTION); - } - - if ((nKind_SurfaceMovement != nMarker_Moving) && !GetSurface_Movement(FLUID_STRUCTURE)) { + if (nKind_SurfaceMovement != nMarker_Moving) { SU2_MPI::Error("Number of KIND_SURFACE_MOVEMENT must match number of MARKER_MOVING", CURRENT_FUNCTION); } @@ -5590,7 +5579,6 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) { case RIGID_MOTION: cout << "rigid mesh motion." << endl; break; case MOVING_HTP: cout << "HTP moving." << endl; break; case ROTATING_FRAME: cout << "rotating reference frame." << endl; break; - case FLUID_STRUCTURE: cout << "fluid-structure motion." << endl; break; case EXTERNAL: cout << "externally prescribed motion." << endl; break; } } @@ -8302,7 +8290,6 @@ bool CConfig::GetVolumetric_Movement() const { if (GetSurface_Movement(AEROELASTIC) || GetSurface_Movement(AEROELASTIC_RIGID_MOTION)|| - GetSurface_Movement(FLUID_STRUCTURE) || GetSurface_Movement(EXTERNAL) || GetSurface_Movement(EXTERNAL_ROTATION)){ volumetric_movement = true; diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index 617c23d7316c..b4c1bc0c868d 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -428,13 +428,13 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return True if the specified vertex is a halo node. */ - bool IsAHaloNode(unsigned short iMarker, unsigned long iVertex); + bool IsAHaloNode(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Get the number of external iterations. * \return Number of external iterations. */ - unsigned long GetnTimeIter(); + unsigned long GetnTimeIter() const; /*! * \brief Get the current external iteration. @@ -446,7 +446,7 @@ class CDriver { * \brief Get the unsteady time step. * \return Unsteady time step. */ - passivedouble GetUnsteady_TimeStep(); + passivedouble GetUnsteady_TimeStep() const; /*! * \brief Get the global index of a vertex on a specified marker. @@ -454,7 +454,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return Vertex global index. */ - unsigned long GetVertexGlobalIndex(unsigned short iMarker, unsigned long iVertex); + unsigned long GetVertexGlobalIndex(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Get undeformed coordinates from the mesh solver. @@ -462,7 +462,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return x,y,z coordinates of the vertex. */ - vector GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex); + vector GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Get the temperature at a vertex on a specified marker. @@ -470,7 +470,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return Temperature of the vertex. */ - passivedouble GetVertexTemperature(unsigned short iMarker, unsigned long iVertex); + passivedouble GetVertexTemperature(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Set the temperature of a vertex on a specified marker. @@ -486,7 +486,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return True if the vertex is a halo node. */ - vector GetVertexHeatFluxes(unsigned short iMarker, unsigned long iVertex); + vector GetVertexHeatFluxes(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Get the wall normal component of the heat flux at a vertex on a specified marker. @@ -494,7 +494,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return Wall normal component of the heat flux at the vertex. */ - passivedouble GetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVertex); + passivedouble GetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Set the wall normal component of the heat flux at a vertex on a specified marker. @@ -510,7 +510,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return Thermal conductivity at the vertex. */ - passivedouble GetThermalConductivity(unsigned short iMarker, unsigned long iVertex); + passivedouble GetThermalConductivity(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Preprocess the inlets via file input for all solvers. @@ -518,8 +518,7 @@ class CDriver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void Inlet_Preprocessing(CSolver ***solver, CGeometry **geometry, - CConfig *config) const; + void Inlet_Preprocessing(CSolver ***solver, CGeometry **geometry, CConfig *config) const; /*! * \brief Get the unit normal (vector) at a vertex on a specified marker. @@ -527,43 +526,43 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return Unit normal (vector) at the vertex. */ - vector GetVertexUnitNormal(unsigned short iMarker, unsigned long iVertex); + vector GetVertexUnitNormal(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Get all the boundary markers tags. * \return List of boundary markers tags. */ - vector GetAllBoundaryMarkersTag(); + vector GetAllBoundaryMarkersTag() const; /*! * \brief Get all the deformable boundary marker tags. * \return List of deformable boundary markers tags. */ - vector GetAllDeformMeshMarkersTag(); + vector GetAllDeformMeshMarkersTag() const; /*! * \brief Get all the heat transfer boundary markers tags. * \return List of heat transfer boundary markers tags. */ - vector GetAllCHTMarkersTag(); + vector GetAllCHTMarkersTag() const; /*! * \brief Get all the (subsonic) inlet boundary markers tags. * \return List of inlet boundary markers tags. */ - vector GetAllInletMarkersTag(); + vector GetAllInletMarkersTag() const; /*! * \brief Get all the boundary markers tags with their associated indices. * \return List of boundary markers tags with their indices. */ - map GetAllBoundaryMarkers(); + map GetAllBoundaryMarkers() const; /*! * \brief Get all the boundary markers tags with their associated types. * \return List of boundary markers tags with their types. */ - map GetAllBoundaryMarkersType(); + map GetAllBoundaryMarkersType() const; /*! * \brief Set the mesh displacement for the elasticity mesh solver. @@ -586,7 +585,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return Vector of sensitivities. */ - vector GetMeshDisp_Sensitivity(unsigned short iMarker, unsigned long iVertex); + vector GetMeshDisp_Sensitivity(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Set the load in X direction for the structural solver. @@ -605,7 +604,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return Vector of displacements. */ - vector GetFEA_Displacements(unsigned short iMarker, unsigned long iVertex); + vector GetFEA_Displacements(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Return the velocities from the FEA Solver. @@ -613,7 +612,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return Vector of velocities. */ - vector GetFEA_Velocity(unsigned short iMarker, unsigned long iVertex); + vector GetFEA_Velocity(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Return the velocities from the FEA Solver. @@ -621,7 +620,7 @@ class CDriver { * \param[in] iVertex - Vertex identifier. * \return Vector of velocities at time n. */ - vector GetFEA_Velocity_n(unsigned short iMarker, unsigned long iVertex); + vector GetFEA_Velocity_n(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Get the sensitivity of the flow loads for the structural solver. @@ -631,7 +630,7 @@ class CDriver { * \param[in] LoadX - Value of the load in the direction Y. * \param[in] LoadX - Value of the load in the direction Z. */ - vector GetFlowLoad_Sensitivity(unsigned short iMarker, unsigned long iVertex); + vector GetFlowLoad_Sensitivity(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Get the flow load (from the extra step - the repeated methods should be unified once the postprocessing @@ -639,7 +638,7 @@ class CDriver { * \param[in] iMarker - Marker identifier. * \param[in] iVertex - Vertex identifier. */ - vector GetFlowLoad(unsigned short iMarker, unsigned long iVertex); + vector GetFlowLoad(unsigned short iMarker, unsigned long iVertex) const; /*! * \brief Set the adjoint of the flow tractions (from the extra step - diff --git a/SU2_CFD/src/iteration/CIteration.cpp b/SU2_CFD/src/iteration/CIteration.cpp index b33b39131f92..d76b49ef187e 100644 --- a/SU2_CFD/src/iteration/CIteration.cpp +++ b/SU2_CFD/src/iteration/CIteration.cpp @@ -34,13 +34,8 @@ void CIteration::SetGrid_Movement(CGeometry** geometry, CSurfaceMovement* surfac CVolumetricMovement* grid_movement, CSolver*** solver, CConfig* config, unsigned long IntIter, unsigned long TimeIter) { unsigned short Kind_Grid_Movement = config->GetKind_GridMovement(); - unsigned long nIterMesh; - bool stat_mesh = true; bool adjoint = config->GetContinuous_Adjoint(); - /*--- Only write to screen if this option is enabled ---*/ - bool Screen_Output = config->GetDeform_Output(); - unsigned short val_iZone = config->GetiZone(); /*--- Perform mesh movement depending on specified type ---*/ @@ -120,32 +115,6 @@ void CIteration::SetGrid_Movement(CGeometry** geometry, CSurfaceMovement* surfac } } - if (config->GetSurface_Movement(FLUID_STRUCTURE)) { - if (rank == MASTER_NODE && Screen_Output) - cout << endl << "Deforming the grid for Fluid-Structure Interaction applications." << endl; - - /*--- Deform the volume grid around the new boundary locations ---*/ - - if (rank == MASTER_NODE && Screen_Output) cout << "Deforming the volume grid." << endl; - grid_movement->SetVolume_Deformation(geometry[MESH_0], config, true, false); - - nIterMesh = grid_movement->Get_nIterMesh(); - stat_mesh = (nIterMesh == 0); - - if (!adjoint && !stat_mesh) { - if (rank == MASTER_NODE && Screen_Output) cout << "Computing grid velocities by finite differencing." << endl; - geometry[MESH_0]->SetGridVelocity(config, TimeIter); - } else if (stat_mesh) { - if (rank == MASTER_NODE && Screen_Output) - cout << "The mesh is up-to-date. Using previously stored grid velocities." << endl; - } - - /*--- Update the multigrid structure after moving the finest grid, - including computing the grid velocities on the coarser levels. ---*/ - - grid_movement->UpdateMultiGrid(geometry, config); - } - if (config->GetSurface_Movement(EXTERNAL) || config->GetSurface_Movement(EXTERNAL_ROTATION)) { /*--- Apply rigid rotation to entire grid first, if necessary ---*/ diff --git a/SU2_CFD/src/iteration/CIterationFactory.cpp b/SU2_CFD/src/iteration/CIterationFactory.cpp index d55188cf2abc..e330c9434e60 100644 --- a/SU2_CFD/src/iteration/CIterationFactory.cpp +++ b/SU2_CFD/src/iteration/CIterationFactory.cpp @@ -112,7 +112,7 @@ CIteration* CIterationFactory::CreateIteration(ENUM_MAIN_SOLVER kindSolver, cons iteration = new CDiscAdjHeatIteration(config); break; - case NO_SOLVER: case FLUID_STRUCTURE_INTERACTION: case TEMPLATE_SOLVER: case MULTIPHYSICS: + case NO_SOLVER: case TEMPLATE_SOLVER: case MULTIPHYSICS: SU2_MPI::Error("No iteration found for specified solver.", CURRENT_FUNCTION); break; } diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 80b4d9aa98e2..3299b75484c8 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -33,10 +33,7 @@ void CDriver::PythonInterface_Preprocessing(CConfig **config, CGeometry ****geometry, CSolver *****solver){ int rank = MASTER_NODE; - -#ifdef HAVE_MPI - MPI_Comm_rank(SU2_MPI::GetComm(), &rank); -#endif + SU2_MPI::Comm_rank(SU2_MPI::GetComm(), &rank); /* --- Initialize boundary conditions customization, this is achieve through the Python wrapper --- */ for(iZone=0; iZone < nZone; iZone++){ @@ -193,7 +190,7 @@ unsigned long CDriver::GetNumberHaloVertices(unsigned short iMarker){ } -unsigned long CDriver::GetVertexGlobalIndex(unsigned short iMarker, unsigned long iVertex) { +unsigned long CDriver::GetVertexGlobalIndex(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint, GlobalIndex; @@ -204,7 +201,7 @@ unsigned long CDriver::GetVertexGlobalIndex(unsigned short iMarker, unsigned lon } -bool CDriver::IsAHaloNode(unsigned short iMarker, unsigned long iVertex) { +bool CDriver::IsAHaloNode(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint; @@ -214,7 +211,7 @@ bool CDriver::IsAHaloNode(unsigned short iMarker, unsigned long iVertex) { } -vector CDriver::GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex) { +vector CDriver::GetInitialMeshCoord(unsigned short iMarker, unsigned long iVertex) const { vector coord(3,0.0); vector coord_passive(3, 0.0); @@ -231,7 +228,7 @@ vector CDriver::GetInitialMeshCoord(unsigned short iMarker, unsig return coord_passive; } -vector CDriver::GetVertexUnitNormal(unsigned short iMarker, unsigned long iVertex){ +vector CDriver::GetVertexUnitNormal(unsigned short iMarker, unsigned long iVertex) const { su2double *Normal; su2double Area; @@ -257,9 +254,9 @@ vector CDriver::GetVertexUnitNormal(unsigned short iMarker, unsig /* Functions to obtain global parameters from SU2 (time steps, delta t, ecc...) */ ////////////////////////////////////////////////////////////////////////////////// -unsigned long CDriver::GetnTimeIter() { +unsigned long CDriver::GetnTimeIter() const { - return config_container[ZONE_0]->GetnTime_Iter(); + return config_container[ZONE_0]->GetnTime_Iter(); } unsigned long CDriver::GetTime_Iter() const{ @@ -267,7 +264,7 @@ unsigned long CDriver::GetTime_Iter() const{ return TimeIter; } -passivedouble CDriver::GetUnsteady_TimeStep(){ +passivedouble CDriver::GetUnsteady_TimeStep() const { return SU2_TYPE::GetValue(config_container[ZONE_0]->GetTime_Step()); } @@ -276,7 +273,7 @@ passivedouble CDriver::GetUnsteady_TimeStep(){ /* Functions related to CHT solver */ /////////////////////////////////////////////////////////////////////////////// -passivedouble CDriver::GetVertexTemperature(unsigned short iMarker, unsigned long iVertex){ +passivedouble CDriver::GetVertexTemperature(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint; su2double vertexWallTemp(0.0); @@ -298,7 +295,7 @@ void CDriver::SetVertexTemperature(unsigned short iMarker, unsigned long iVertex geometry_container[ZONE_0][INST_0][MESH_0]->SetCustomBoundaryTemperature(iMarker, iVertex, val_WallTemp); } -vector CDriver::GetVertexHeatFluxes(unsigned short iMarker, unsigned long iVertex){ +vector CDriver::GetVertexHeatFluxes(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint; unsigned short iDim; @@ -313,7 +310,6 @@ vector CDriver::GetVertexHeatFluxes(unsigned short iMarker, unsig vector HeatFluxPassive (3,0.0); bool compressible = (config_container[ZONE_0]->GetKind_Regime() == COMPRESSIBLE); - bool halo; iPoint = geometry_container[ZONE_0][INST_0][MESH_0]->vertex[iMarker][iVertex]->GetNode(); @@ -333,7 +329,7 @@ vector CDriver::GetVertexHeatFluxes(unsigned short iMarker, unsig return HeatFluxPassive; } -passivedouble CDriver::GetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVertex){ +passivedouble CDriver::GetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVertex) const{ unsigned long iPoint; unsigned short iDim; @@ -381,7 +377,7 @@ void CDriver::SetVertexNormalHeatFlux(unsigned short iMarker, unsigned long iVer geometry_container[ZONE_0][INST_0][MESH_0]->SetCustomBoundaryHeatFlux(iMarker, iVertex, val_WallHeatFlux); } -passivedouble CDriver::GetThermalConductivity(unsigned short iMarker, unsigned long iVertex){ +passivedouble CDriver::GetThermalConductivity(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint; su2double Prandtl_Lam = config_container[ZONE_0]->GetPrandtl_Lam(); @@ -403,7 +399,7 @@ passivedouble CDriver::GetThermalConductivity(unsigned short iMarker, unsigned l /* Functions related to the management of markers */ //////////////////////////////////////////////////////////////////////////////// -vector CDriver::GetAllBoundaryMarkersTag(){ +vector CDriver::GetAllBoundaryMarkersTag() const { vector boundariesTagList; unsigned short iMarker,nBoundariesMarkers; @@ -420,7 +416,7 @@ vector CDriver::GetAllBoundaryMarkersTag(){ return boundariesTagList; } -vector CDriver::GetAllDeformMeshMarkersTag(){ +vector CDriver::GetAllDeformMeshMarkersTag() const { vector interfaceBoundariesTagList; unsigned short iMarker, nBoundariesMarker; @@ -437,7 +433,7 @@ vector CDriver::GetAllDeformMeshMarkersTag(){ return interfaceBoundariesTagList; } -vector CDriver::GetAllCHTMarkersTag(){ +vector CDriver::GetAllCHTMarkersTag() const { vector CHTBoundariesTagList; unsigned short iMarker, nBoundariesMarker; @@ -456,7 +452,7 @@ vector CDriver::GetAllCHTMarkersTag(){ return CHTBoundariesTagList; } -vector CDriver::GetAllInletMarkersTag(){ +vector CDriver::GetAllInletMarkersTag() const { vector BoundariesTagList; unsigned short iMarker, nBoundariesMarker; @@ -476,7 +472,7 @@ vector CDriver::GetAllInletMarkersTag(){ return BoundariesTagList; } -map CDriver::GetAllBoundaryMarkers(){ +map CDriver::GetAllBoundaryMarkers() const { map allBoundariesMap; unsigned short iMarker, nBoundaryMarkers; @@ -492,7 +488,7 @@ map CDriver::GetAllBoundaryMarkers(){ return allBoundariesMap; } -map CDriver::GetAllBoundaryMarkersType(){ +map CDriver::GetAllBoundaryMarkersType() const { map allBoundariesTypeMap; unsigned short iMarker, KindBC; @@ -647,7 +643,7 @@ void CDriver::SetFEA_Loads(unsigned short iMarker, unsigned long iVertex, passiv } -vector CDriver::GetFEA_Displacements(unsigned short iMarker, unsigned long iVertex) { +vector CDriver::GetFEA_Displacements(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint; vector Displacements(3, 0.0); @@ -672,7 +668,7 @@ vector CDriver::GetFEA_Displacements(unsigned short iMarker, unsi } -vector CDriver::GetFEA_Velocity(unsigned short iMarker, unsigned long iVertex) { +vector CDriver::GetFEA_Velocity(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint; vector Velocity(3, 0.0); @@ -698,7 +694,7 @@ vector CDriver::GetFEA_Velocity(unsigned short iMarker, unsigned return Velocity_passive; } -vector CDriver::GetFEA_Velocity_n(unsigned short iMarker, unsigned long iVertex) { +vector CDriver::GetFEA_Velocity_n(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint; vector Velocity_n(3, 0.0); @@ -729,7 +725,7 @@ vector CDriver::GetFEA_Velocity_n(unsigned short iMarker, unsigne /* Functions related to adjoint simulations */ //////////////////////////////////////////////////////////////////////////////// -vector CDriver::GetMeshDisp_Sensitivity(unsigned short iMarker, unsigned long iVertex) { +vector CDriver::GetMeshDisp_Sensitivity(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint; vector Disp_Sens(3, 0.0); @@ -754,9 +750,7 @@ vector CDriver::GetMeshDisp_Sensitivity(unsigned short iMarker, u } - - -vector CDriver::GetFlowLoad_Sensitivity(unsigned short iMarker, unsigned long iVertex) { +vector CDriver::GetFlowLoad_Sensitivity(unsigned short iMarker, unsigned long iVertex) const { unsigned long iPoint; vector FlowLoad_Sens(3, 0.0); @@ -842,7 +836,7 @@ void CDriver::CommunicateMeshDisplacement(void) { /* Functions related to flow loads */ //////////////////////////////////////////////////////////////////////////////// -vector CDriver::GetFlowLoad(unsigned short iMarker, unsigned long iVertex) { +vector CDriver::GetFlowLoad(unsigned short iMarker, unsigned long iVertex) const { vector FlowLoad(3, 0.0); vector FlowLoad_passive(3, 0.0); diff --git a/TestCases/cont_adj_rans/oneram6/turb_ONERAM6.cfg b/TestCases/cont_adj_rans/oneram6/turb_ONERAM6.cfg index 5c7f5c8dcdd8..ebfc62d342cd 100644 --- a/TestCases/cont_adj_rans/oneram6/turb_ONERAM6.cfg +++ b/TestCases/cont_adj_rans/oneram6/turb_ONERAM6.cfg @@ -11,12 +11,7 @@ % ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% % -% Physical governing equations (EULER, NAVIER_STOKES, -% PLASMA_EULER, PLASMA_NAVIER_STOKES, -% FREE_SURFACE_EULER, FREE_SURFACE_NAVIER_STOKES, -% FLUID_STRUCTURE_EULER, FLUID_STRUCTURE_NAVIER_STOKES, -% AEROACOUSTIC_EULER, AEROACOUSTIC_NAVIER_STOKES, -% WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY) +% Physical governing equations (EULER, NAVIER_STOKES, etc.) SOLVER= NAVIER_STOKES % % Specify turbulence model (NONE, SA, SA_NEG, SST) diff --git a/TestCases/gust/inv_gust_NACA0012.cfg b/TestCases/gust/inv_gust_NACA0012.cfg index 61950875cf14..c3d44d246925 100644 --- a/TestCases/gust/inv_gust_NACA0012.cfg +++ b/TestCases/gust/inv_gust_NACA0012.cfg @@ -72,7 +72,7 @@ INNER_ITER= 100 % ----------------------- DYNAMIC MESH DEFINITION -----------------------------% % % Type of dynamic mesh (NONE, RIGID_MOTION, DEFORMING, ROTATING_FRAME, -% MOVING_WALL, FLUID_STRUCTURE, AEROELASTIC, ELASTICITY, +% MOVING_WALL, AEROELASTIC, ELASTICITY, % EXTERNAL, AEROELASTIC_RIGID_MOTION) GRID_MOVEMENT= GUST % diff --git a/TestCases/moving_wall/cavity/lam_cavity.cfg b/TestCases/moving_wall/cavity/lam_cavity.cfg index 545e44b297b5..f6f7378821e0 100644 --- a/TestCases/moving_wall/cavity/lam_cavity.cfg +++ b/TestCases/moving_wall/cavity/lam_cavity.cfg @@ -48,8 +48,7 @@ REYNOLDS_LENGTH= 1.0 % ----------------------- DYNAMIC MESH DEFINITION -----------------------------% % % Type of dynamic mesh (NONE, RIGID_MOTION, DEFORMING, ROTATING_FRAME, -% MOVING_WALL, FLUID_STRUCTURE, AEROELASTIC, ELASTICITY, -% EXTERNAL) +% MOVING_WALL, AEROELASTIC, ELASTICITY, EXTERNAL) SURFACE_MOVEMENT= MOVING_WALL % % Motion mach number (non-dimensional). Used for initializing a viscous flow diff --git a/TestCases/moving_wall/spinning_cylinder/spinning_cylinder.cfg b/TestCases/moving_wall/spinning_cylinder/spinning_cylinder.cfg index 4dc1a80d0a7b..534193074e4e 100644 --- a/TestCases/moving_wall/spinning_cylinder/spinning_cylinder.cfg +++ b/TestCases/moving_wall/spinning_cylinder/spinning_cylinder.cfg @@ -48,7 +48,7 @@ REYNOLDS_LENGTH= 1.0 % ----------------------- DYNAMIC MESH DEFINITION -----------------------------% % % Type of dynamic mesh (NONE, RIGID_MOTION, DEFORMING, ROTATING_FRAME, -% MOVING_WALL, FLUID_STRUCTURE, AEROELASTIC, EXTERNAL) +% MOVING_WALL, AEROELASTIC, EXTERNAL) SURFACE_MOVEMENT= MOVING_WALL % % Motion mach number (non-dimensional). Used for intitializing a viscous flow diff --git a/TestCases/optimization_euler/pitching_naca64a010/pitching_NACA64A010.cfg b/TestCases/optimization_euler/pitching_naca64a010/pitching_NACA64A010.cfg index 8d22b6507446..f5fa022f1b4e 100644 --- a/TestCases/optimization_euler/pitching_naca64a010/pitching_NACA64A010.cfg +++ b/TestCases/optimization_euler/pitching_naca64a010/pitching_NACA64A010.cfg @@ -51,7 +51,7 @@ UNST_ADJOINT_ITER= 251 % Dynamic mesh simulation (NO, YES) GRID_MOVEMENT= YES % -% Type of mesh motion (NONE, FLUTTER, RIGID_MOTION, FLUID_STRUCTURE) +% Type of mesh motion (NONE, FLUTTER, RIGID_MOTION) GRID_MOVEMENT_KIND= RIGID_MOTION % % Motion mach number (non-dimensional). Used for initializing a viscous flow diff --git a/TestCases/optimization_euler/pitching_oneram6/pitching_ONERAM6.cfg b/TestCases/optimization_euler/pitching_oneram6/pitching_ONERAM6.cfg index 8a1ab1851c2b..d9ca2a72dd7a 100644 --- a/TestCases/optimization_euler/pitching_oneram6/pitching_ONERAM6.cfg +++ b/TestCases/optimization_euler/pitching_oneram6/pitching_ONERAM6.cfg @@ -65,7 +65,7 @@ UNST_ADJOINT_ITER= 251 GRID_MOVEMENT= YES % % Type of dynamic mesh (NONE, RIGID_MOTION, DEFORMING, ROTATING_FRAME, -% MOVING_WALL, FLUID_STRUCTURE, AEROELASTIC, EXTERNAL) +% MOVING_WALL, AEROELASTIC, EXTERNAL) GRID_MOVEMENT_KIND= RIGID_MOTION % % Motion mach number (non-dimensional). Used for intitializing a viscous flow diff --git a/TestCases/optimization_euler/rotating_naca0012/rotating_NACA0012.cfg b/TestCases/optimization_euler/rotating_naca0012/rotating_NACA0012.cfg index ff3a00c420f0..83801d29d8c5 100644 --- a/TestCases/optimization_euler/rotating_naca0012/rotating_NACA0012.cfg +++ b/TestCases/optimization_euler/rotating_naca0012/rotating_NACA0012.cfg @@ -58,7 +58,7 @@ REF_AREA= 1.0 GRID_MOVEMENT= YES % % Type of dynamic mesh (NONE, RIGID_MOTION, DEFORMING, ROTATING_FRAME, -% MOVING_WALL, FLUID_STRUCTURE, AEROELASTIC, EXTERNAL) +% MOVING_WALL, AEROELASTIC, EXTERNAL) GRID_MOVEMENT_KIND= ROTATING_FRAME % % Motion mach number (non-dimensional). Used for intitializing a viscous flow From f956724c55365ae53b9df7577a48db3ea9ff965f Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 19 Feb 2021 14:35:24 +0000 Subject: [PATCH 31/32] fix #1202 --- SU2_CFD/src/output/CFlowCompOutput.cpp | 7 ------- SU2_CFD/src/output/CFlowOutput.cpp | 4 ++++ SU2_CFD/src/output/CNEMOCompOutput.cpp | 7 ------- 3 files changed, 4 insertions(+), 14 deletions(-) diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index 3b6f328b533f..14240f531057 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -279,9 +279,6 @@ void CFlowCompOutput::SetHistoryOutputFields(CConfig *config){ Add_CpInverseDesignOutput(config); - /*--- Add combo obj value --- */ - - AddHistoryOutput("COMBO", "ComboObj", ScreenOutputFormat::SCIENTIFIC, "COMBO", "Combined obj. function value.", HistoryFieldType::COEFFICIENT); } void CFlowCompOutput::SetVolumeOutputFields(CConfig *config){ @@ -718,10 +715,6 @@ void CFlowCompOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSol Set_CpInverseDesign(flow_solver, geometry, config); - /*--- Set combo obj value --- */ - - SetHistoryOutputValue("COMBO", flow_solver->GetTotal_ComboObj()); - } bool CFlowCompOutput::SetInit_Residuals(CConfig *config){ diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index e0b28799cf5e..4c33d0354a26 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -638,6 +638,8 @@ void CFlowOutput::AddAerodynamicCoefficients(CConfig *config){ /// DESCRIPTION: Angle of attack AddHistoryOutput("AOA", "AoA", ScreenOutputFormat::FIXED, "AOA", "Angle of attack"); + + AddHistoryOutput("COMBO", "ComboObj", ScreenOutputFormat::SCIENTIFIC, "COMBO", "Combined obj. function value.", HistoryFieldType::COEFFICIENT); } void CFlowOutput::SetAerodynamicCoefficients(CConfig *config, CSolver *flow_solver){ @@ -680,6 +682,8 @@ void CFlowOutput::SetAerodynamicCoefficients(CConfig *config, CSolver *flow_solv } SetHistoryOutputValue("AOA", config->GetAoA()); + + SetHistoryOutputValue("COMBO", flow_solver->GetTotal_ComboObj()); } void CFlowOutput::SetRotatingFrameCoefficients(CConfig *config, CSolver *flow_solver) { diff --git a/SU2_CFD/src/output/CNEMOCompOutput.cpp b/SU2_CFD/src/output/CNEMOCompOutput.cpp index 917823fa4fd4..43686451a968 100644 --- a/SU2_CFD/src/output/CNEMOCompOutput.cpp +++ b/SU2_CFD/src/output/CNEMOCompOutput.cpp @@ -279,9 +279,6 @@ void CNEMOCompOutput::SetHistoryOutputFields(CConfig *config){ Add_CpInverseDesignOutput(config); - /*--- Add combo obj value --- */ - - AddHistoryOutput("COMBO", "ComboObj", ScreenOutputFormat::SCIENTIFIC, "COMBO", "Combined obj. function value.", HistoryFieldType::COEFFICIENT); } void CNEMOCompOutput::SetVolumeOutputFields(CConfig *config){ @@ -686,10 +683,6 @@ void CNEMOCompOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSol Set_CpInverseDesign(NEMO_solver, geometry, config); - /*--- Set combo obj value --- */ - - SetHistoryOutputValue("COMBO", NEMO_solver->GetTotal_ComboObj()); - } bool CNEMOCompOutput::SetInit_Residuals(CConfig *config){ From 392e31f89dcfede060dad906421f4d42e96dfcc9 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 19 Feb 2021 14:48:32 +0000 Subject: [PATCH 32/32] more const --- SU2_CFD/include/drivers/CDriver.hpp | 18 +++++++++--------- SU2_CFD/src/python_wrapper_structure.cpp | 18 +++++++++--------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index b4c1bc0c868d..954f943036fd 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -370,57 +370,57 @@ class CDriver { * \brief Get the total drag. * \return Total drag. */ - passivedouble Get_Drag(); + passivedouble Get_Drag() const; /*! * \brief Get the total lift. * \return Total lift. */ - passivedouble Get_Lift(); + passivedouble Get_Lift() const; /*! * \brief Get the total x moment. * \return Total x moment. */ - passivedouble Get_Mx(); + passivedouble Get_Mx() const; /*! * \brief Get the total y moment. * \return Total y moment. */ - passivedouble Get_My(); + passivedouble Get_My() const; /*! * \brief Get the total z moment. * \return Total z moment. */ - passivedouble Get_Mz(); + passivedouble Get_Mz() const; /*! * \brief Get the total drag coefficient. * \return Total drag coefficient. */ - passivedouble Get_DragCoeff(); + passivedouble Get_DragCoeff() const; /*! * \brief Get the total lift coefficient. * \return Total lift coefficient. */ - passivedouble Get_LiftCoeff(); + passivedouble Get_LiftCoeff() const; /*! * \brief Get the number of vertices (halo nodes included) from a specified marker. * \param[in] iMarker - Marker identifier. * \return Number of vertices. */ - unsigned long GetNumberVertices(unsigned short iMarker); + unsigned long GetNumberVertices(unsigned short iMarker) const; /*! * \brief Get the number of halo vertices from a specified marker. * \param[in] iMarker - Marker identifier. * \return Number of vertices. */ - unsigned long GetNumberHaloVertices(unsigned short iMarker); + unsigned long GetNumberHaloVertices(unsigned short iMarker) const; /*! * \brief Check if a vertex is physical or not (halo node) on a specified marker. diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 3299b75484c8..ebf9a926d5d2 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -63,7 +63,7 @@ void CDriver::PythonInterface_Preprocessing(CConfig **config, CGeometry ****geom /* Functions related to the global performance indices (Lift, Drag, ecc..) */ ///////////////////////////////////////////////////////////////////////////// -passivedouble CDriver::Get_Drag() { +passivedouble CDriver::Get_Drag() const { unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); @@ -78,7 +78,7 @@ passivedouble CDriver::Get_Drag() { return SU2_TYPE::GetValue(val_Drag); } -passivedouble CDriver::Get_Lift() { +passivedouble CDriver::Get_Lift() const { unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); @@ -93,7 +93,7 @@ passivedouble CDriver::Get_Lift() { return SU2_TYPE::GetValue(val_Lift); } -passivedouble CDriver::Get_Mx(){ +passivedouble CDriver::Get_Mx() const { unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); @@ -110,7 +110,7 @@ passivedouble CDriver::Get_Mx(){ return SU2_TYPE::GetValue(val_Mx); } -passivedouble CDriver::Get_My(){ +passivedouble CDriver::Get_My() const { unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); @@ -127,7 +127,7 @@ passivedouble CDriver::Get_My(){ return SU2_TYPE::GetValue(val_My); } -passivedouble CDriver::Get_Mz() { +passivedouble CDriver::Get_Mz() const { unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); @@ -144,7 +144,7 @@ passivedouble CDriver::Get_Mz() { return SU2_TYPE::GetValue(val_Mz); } -passivedouble CDriver::Get_DragCoeff() { +passivedouble CDriver::Get_DragCoeff() const { unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); @@ -155,7 +155,7 @@ passivedouble CDriver::Get_DragCoeff() { return SU2_TYPE::GetValue(CDrag); } -passivedouble CDriver::Get_LiftCoeff() { +passivedouble CDriver::Get_LiftCoeff() const { unsigned short val_iZone = ZONE_0; unsigned short FinestMesh = config_container[val_iZone]->GetFinestMesh(); @@ -170,13 +170,13 @@ passivedouble CDriver::Get_LiftCoeff() { /* Functions to obtain information from the geometry/mesh */ ///////////////////////////////////////////////////////////////////////////// -unsigned long CDriver::GetNumberVertices(unsigned short iMarker){ +unsigned long CDriver::GetNumberVertices(unsigned short iMarker) const { return geometry_container[ZONE_0][INST_0][MESH_0]->nVertex[iMarker]; } -unsigned long CDriver::GetNumberHaloVertices(unsigned short iMarker){ +unsigned long CDriver::GetNumberHaloVertices(unsigned short iMarker) const { unsigned long nHaloVertices, iVertex, iPoint;