diff --git a/modules/aerodyn/python-lib/aerodyn_inflow_library.py b/modules/aerodyn/python-lib/aerodyn_inflow_library.py index 508d42c47a..9539ddff68 100644 --- a/modules/aerodyn/python-lib/aerodyn_inflow_library.py +++ b/modules/aerodyn/python-lib/aerodyn_inflow_library.py @@ -38,7 +38,6 @@ c_char_p, c_wchar, c_wchar_p, - c_bool ) import numpy as np import datetime @@ -100,8 +99,9 @@ def __init__(self, library_path): self.MSL2SWL = 0.0 # Offset between still-water level and mean sea level (m) [positive upward] # flags - self.storeHHVel = False - self.transposeDCM= False + self.storeHHVel = 1 # 0=false, 1=true + self.transposeDCM= 1 # 0=false, 1=true + self.debuglevel = 0 # 0-4 levels # VTK self.WrVTK = 0 # default of no vtk output @@ -126,12 +126,10 @@ def __init__(self, library_path): # number of output channels self.numChannels = 0 # Number of channels returned - # Aero calculation method -- AeroProjMod - # APM_BEM_NoSweepPitchTwist - 1 - "Original AeroDyn model where momentum balance is done in the WithoutSweepPitchTwist system" - # APM_BEM_Polar - 2 - "Use staggered polar grid for momentum balance in each annulus" - # APM_LiftingLine - 3 - "Use the blade lifting line (i.e. the structural) orientation (currently for OLAF with VAWT)" - self.AeroProjMod = 1 + # Number of turbines + self.numTurbines = 1 +#FIXME: some assumptions about a single turbine here. # Initial position of hub and blades # used for setup of AD, not used after init. self.initHubPos = np.zeros(shape=(3),dtype=c_float) @@ -159,11 +157,41 @@ def __init__(self, library_path): # _initialize_routines() ------------------------------------------------------------------------------------------------------------ def _initialize_routines(self): - self.AeroDyn_Inflow_C_Init.argtypes = [ - POINTER(c_bool), # AD input file passed as string + # initialize data storage in library for multiple turbines + self.ADI_C_PreInit.argtypes = [ + POINTER(c_int), # numTurbines + POINTER(c_int), # transposeDCM + POINTER(c_int), # debuglevel + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C + ] + self.ADI_C_PreInit.restype = c_int + + # setup one rotor + self.ADI_C_SetupRotor.argtypes = [ + POINTER(c_int), # iturb + POINTER(c_int), # isHAWT + POINTER(c_float), # Turb_RefPos + POINTER(c_float), # initHubPos + POINTER(c_double), # initHubOrient_flat + POINTER(c_float), # initNacellePos + POINTER(c_double), # initNacelleOrient_flat + POINTER(c_int), # numBlades + POINTER(c_float), # initRootPos_flat + POINTER(c_double), # initRootOrient_flat + POINTER(c_int), # numMeshPts + POINTER(c_float), # initMeshPos_flat + POINTER(c_double), # initMeshOrient_flat + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C + ] + + # initialize ADI with data set by PreInit and SetupRotor + self.ADI_C_Init.argtypes = [ + POINTER(c_int), # AD input file passed as string POINTER(c_char_p), # AD input file as string POINTER(c_int), # AD input file string length - POINTER(c_bool), # IfW input file passed as string + POINTER(c_int), # IfW input file passed as string POINTER(c_char_p), # IfW input file as string POINTER(c_int), # IfW input file string length POINTER(c_char), # OutRootName @@ -175,47 +203,35 @@ def _initialize_routines(self): POINTER(c_float), # defPvap POINTER(c_float), # WtrDpth POINTER(c_float), # MSL2SWL - POINTER(c_int), # AeroProjMod POINTER(c_int), # InterpOrder POINTER(c_double), # dt POINTER(c_double), # tmax - POINTER(c_bool), # storeHHVel - POINTER(c_bool), # transposeDCM + POINTER(c_int), # storeHHVel POINTER(c_int), # WrVTK POINTER(c_int), # WrVTK_Type POINTER(c_float), # VTKNacDim POINTER(c_float), # VTKHubRad POINTER(c_int), # wrOuts -- file format for writing outputs POINTER(c_double), # DT_Outs -- timestep for outputs to file - POINTER(c_float), # initHubPos - POINTER(c_double), # initHubOrient_flat - POINTER(c_float), # initNacellePos - POINTER(c_double), # initNacelleOrient_flat - POINTER(c_int), # numBlades - POINTER(c_float), # initRootPos_flat - POINTER(c_double), # initRootOrient_flat - POINTER(c_int), # numMeshPts - POINTER(c_float), # initMeshPos_flat - POINTER(c_double), # initMeshOrient_flat POINTER(c_int), # number of channels POINTER(c_char), # output channel names POINTER(c_char), # output channel units POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.AeroDyn_Inflow_C_Init.restype = c_int + self.ADI_C_Init.restype = c_int - #self.AeroDyn_Inflow_C_ReInit.argtypes = [ + #self.ADI_C_ReInit.argtypes = [ # POINTER(c_double), # t_initial # POINTER(c_double), # dt # POINTER(c_double), # tmax # POINTER(c_int), # ErrStat_C # POINTER(c_char) # ErrMsg_C #] - #self.AeroDyn_Inflow_C_ReInit.restype = c_int + #self.ADI_C_ReInit.restype = c_int - self.AeroDyn_Inflow_C_CalcOutput.argtypes = [ - POINTER(c_double), # Time_C + self.ADI_C_SetRotorMotion.argtypes = [ + POINTER(c_int), # iturb POINTER(c_float), # HubPos POINTER(c_double), # HubOrient_flat POINTER(c_float), # HubVel @@ -233,50 +249,104 @@ def _initialize_routines(self): POINTER(c_double), # MeshOrient_flat POINTER(c_float), # MeshVel POINTER(c_float), # MeshAcc + ] + + + self.ADI_C_GetRotorLoads.argtypes = [ + POINTER(c_int), # iturb + POINTER(c_int), # numMeshPts POINTER(c_float), # meshFrc -- mesh forces/moments in flat array of 6*numMeshPts + POINTER(c_int), # ErrStat_C + POINTER(c_char) # ErrMsg_C + ] + self.ADI_C_GetRotorLoads.restype = c_int + + + self.ADI_C_CalcOutput.argtypes = [ + POINTER(c_double), # Time_C POINTER(c_float), # Output Channel Values POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.AeroDyn_Inflow_C_CalcOutput.restype = c_int + self.ADI_C_CalcOutput.restype = c_int + - self.AeroDyn_Inflow_C_UpdateStates.argtypes = [ + self.ADI_C_UpdateStates.argtypes = [ POINTER(c_double), # Time_C POINTER(c_double), # TimeNext_C - POINTER(c_float), # HubPos - POINTER(c_double), # HubOrient_flat - POINTER(c_float), # HubVel - POINTER(c_float), # HubAcc - POINTER(c_float), # NacPos - POINTER(c_double), # NacOrient_flat - POINTER(c_float), # NacVel - POINTER(c_float), # NacAcc - POINTER(c_float), # RootPos - POINTER(c_double), # RootOrient_flat - POINTER(c_float), # RootVel - POINTER(c_float), # RootAcc - POINTER(c_int), # numMeshPts - POINTER(c_float), # MeshPos - POINTER(c_double), # MeshOrient_flat - POINTER(c_float), # MeshVel - POINTER(c_float), # MeshAcc POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.AeroDyn_Inflow_C_UpdateStates.restype = c_int + self.ADI_C_UpdateStates.restype = c_int - self.AeroDyn_Inflow_C_End.argtypes = [ + self.ADI_C_End.argtypes = [ POINTER(c_int), # ErrStat_C POINTER(c_char) # ErrMsg_C ] - self.AeroDyn_Inflow_C_End.restype = c_int + self.ADI_C_End.restype = c_int - # aerodyn_inflow_init ------------------------------------------------------------------------------------------------------------ - def aerodyn_inflow_init(self, AD_input_string_array, IfW_input_string_array): - # some bookkeeping initialization - self._numChannels_c = c_int(0) + # adi_init ------------------------------------------------------------------------------------------------------------ + def adi_preinit(self): + # Pass number of turbines over to setup arrays. + + # call ADI_C_PreInit + self.ADI_C_PreInit( + byref(c_int(self.numTurbines)), # IN: numTurbines + byref(c_int(self.transposeDCM)), # IN: transposeDCM + byref(c_int(self.debuglevel)), # IN: debuglevel + byref(self.error_status_c), # OUT: ErrStat_C + self.error_message_c # OUT: ErrMsg_C + ) + + self.check_error() + + def adi_setuprotor(self,iturb,isHAWT,turbRefPos): + # setup one rotor with initial root/mesh info self._initNumMeshPts = self.numMeshPts self._initNumBlades = self.numBlades + _turbRefPos = (c_float * len(turbRefPos))(*turbRefPos) + + # check hub and root points for initialization + self.check_init_hubroot() + + # Check initial mesh positions + self.check_init_mesh() + + # Flatten arrays to pass + # [x2,y1,z1, x2,y2,z2 ...] + initHubPos_c = (c_float * len(self.initHubPos ))(*self.initHubPos ) + initHubOrient_c = (c_double * len(self.initHubOrient ))(*self.initHubOrient ) + initNacellePos_c = (c_float * len(self.initNacellePos ))(*self.initNacellePos ) + initNacelleOrient_c = (c_double * len(self.initNacelleOrient))(*self.initNacelleOrient) + initRootPos_flat_c = self.flatPosArr( self._initNumBlades, self.numBlades,self.initRootPos, 'RootPos') + initRootOrient_flat_c = self.flatOrientArr(self._initNumBlades, self.numBlades,self.initRootOrient, 'RootOrient') + initMeshPos_flat_c = self.flatPosArr( self._initNumMeshPts,self.numMeshPts,self.initMeshPos, 'MeshPos') + initMeshOrient_flat_c = self.flatOrientArr(self._initNumMeshPts,self.numMeshPts,self.initMeshOrient,'MeshOrient') + + self.ADI_C_SetupRotor( + c_int(iturb), # IN: iturb -- current turbine number + c_int(isHAWT), # IN: 1: is HAWT, 0: VAWT or cross-flow + _turbRefPos, # IN: turbine reference position + initHubPos_c, # IN: initHubPos -- initial hub position + initHubOrient_c, # IN: initHubOrient -- initial hub orientation DCM in flat array of 9 elements + initNacellePos_c, # IN: initNacellePos -- initial hub position + initNacelleOrient_c, # IN: initNacelleOrient -- initial hub orientation DCM in flat array of 9 elements + byref(c_int(self.numBlades)), # IN: number of blades (matches number of blade root positions) + initRootPos_flat_c, # IN: initBladeRootPos -- initial node positions in flat array of 3*numBlades + initRootOrient_flat_c, # IN: initBladeRootOrient -- initial blade root orientation DCMs in flat array of 9*numBlades + byref(c_int(self.numMeshPts)), # IN: number of mesh points expected + initMeshPos_flat_c, # IN: initMeshPos -- initial node positions in flat array of 3*numMeshPts + initMeshOrient_flat_c, # IN: initMeshOrient -- initial node orientation DCMs in flat array of 9*numMeshPts + byref(self.error_status_c), # OUT: ErrStat_C + self.error_message_c # OUT: ErrMsg_C + ) + + self.check_error() + + + def adi_init(self, AD_input_string_array, IfW_input_string_array): + # some bookkeeping initialization + self._numChannels_c = c_int(0) # Primary input file will be passed as a single string joined by # C_NULL_CHAR. @@ -293,31 +363,16 @@ def aerodyn_inflow_init(self, AD_input_string_array, IfW_input_string_array): # Rootname for ADI output files (echo etc). _outRootName_c = create_string_buffer((self.outRootName.ljust(self.default_str_c_len)).encode('utf-8')) - # check hub and root points for initialization - self.check_init_hubroot() - - # Check initial mesh positions - self.check_init_mesh() - # Flatten arrays to pass # [x2,y1,z1, x2,y2,z2 ...] VTKNacDim_c = (c_float * len(self.VTKNacDim ))(*self.VTKNacDim ) - initHubPos_c = (c_float * len(self.initHubPos ))(*self.initHubPos ) - initHubOrient_c = (c_double * len(self.initHubOrient ))(*self.initHubOrient ) - initNacellePos_c = (c_float * len(self.initNacellePos ))(*self.initNacellePos ) - initNacelleOrient_c = (c_double * len(self.initNacelleOrient))(*self.initNacelleOrient) - initRootPos_flat_c = self.flatPosArr( self._initNumBlades, self.numBlades,self.initRootPos, 'Init','RootPos') - initRootOrient_flat_c = self.flatOrientArr(self._initNumBlades, self.numBlades,self.initRootOrient, 'Init','RootOrient') - initMeshPos_flat_c = self.flatPosArr( self._initNumMeshPts,self.numMeshPts,self.initMeshPos, 'Init','MeshPos') - initMeshOrient_flat_c = self.flatOrientArr(self._initNumMeshPts,self.numMeshPts,self.initMeshOrient,'Init','MeshOrient') - - # call AeroDyn_Inflow_C_Init - self.AeroDyn_Inflow_C_Init( - byref(c_bool(self.ADinputPass)), # IN: AD input file is passed + # call ADI_C_Init + self.ADI_C_Init( + byref(c_int(self.ADinputPass)), # IN: AD input file is passed c_char_p(AD_input_string), # IN: AD input file as string (or filename if ADinputPass is false) byref(c_int(AD_input_string_length)), # IN: AD input file string length - byref(c_bool(self.IfWinputPass)), # IN: IfW input file is passed + byref(c_int(self.IfWinputPass)), # IN: IfW input file is passed c_char_p(IfW_input_string), # IN: IfW input file as string (or filename if IfWinputPass is false) byref(c_int(IfW_input_string_length)), # IN: IfW input file string length _outRootName_c, # IN: rootname for ADI file writing @@ -329,28 +384,16 @@ def aerodyn_inflow_init(self, AD_input_string_array, IfW_input_string_array): byref(c_float(self.defPvap)), # IN: defPvap byref(c_float(self.WtrDpth)), # IN: WtrDpth byref(c_float(self.MSL2SWL)), # IN: MSL2SWL - byref(c_int(self.AeroProjMod)), # IN: AeroProjMod byref(c_int(self.InterpOrder)), # IN: InterpOrder (1: linear, 2: quadratic) byref(c_double(self.dt)), # IN: time step (dt) byref(c_double(self.tmax)), # IN: tmax - byref(c_bool(self.storeHHVel)), # IN: storeHHVel - byref(c_bool(self.transposeDCM)), # IN: transposeDCM + byref(c_int(self.storeHHVel)), # IN: storeHHVel byref(c_int(self.WrVTK)), # IN: WrVTK byref(c_int(self.WrVTK_Type)), # IN: WrVTK_Type VTKNacDim_c, # IN: VTKNacDim byref(c_float(self.VTKHubRad)), # IN: VTKHubRad byref(c_int(self.wrOuts)), # IN: wrOuts -- file format for writing outputs byref(c_double(self.DT_Outs)), # IN: DT_Outs -- timestep for outputs to file - initHubPos_c, # IN: initHubPos -- initial hub position - initHubOrient_c, # IN: initHubOrient -- initial hub orientation DCM in flat array of 9 elements - initNacellePos_c, # IN: initNacellePos -- initial hub position - initNacelleOrient_c, # IN: initNacelleOrient -- initial hub orientation DCM in flat array of 9 elements - byref(c_int(self.numBlades)), # IN: number of blades (matches number of blade root positions) - initRootPos_flat_c, # IN: initBladeRootPos -- initial node positions in flat array of 3*numBlades - initRootOrient_flat_c, # IN: initBladeRootOrient -- initial blade root orientation DCMs in flat array of 9*numBlades - byref(c_int(self.numMeshPts)), # IN: number of attachment points expected (where motions are transferred into HD) - initMeshPos_flat_c, # IN: initMeshPos -- initial node positions in flat array of 3*numMeshPts - initMeshOrient_flat_c, # IN: initMeshOrient -- initial node orientation DCMs in flat array of 9*numMeshPts byref(self._numChannels_c), # OUT: number of channels self._channel_names_c, # OUT: output channel names self._channel_units_c, # OUT: output channel units @@ -364,12 +407,12 @@ def aerodyn_inflow_init(self, AD_input_string_array, IfW_input_string_array): self.numChannels = self._numChannels_c.value - ## aerodyn_inflow_reinit ------------------------------------------------------------------------------------------------------------ - #def aerodyn_inflow_reinit(self): - # #FIXME: need to pass something in here I think. Not sure what. + ## adi_reinit ------------------------------------------------------------------------------------------------------------ + #FIXME: this routine is not setup + #def adi_reinit(self): # - # # call AeroDyn_Inflow_C_ReInit - # self.AeroDyn_Inflow_C_ReInit( + # # call ADI_C_ReInit + # self.ADI_C_ReInit( # byref(c_double(self.dt)), # IN: time step (dt) # byref(c_double(self.tmax)), # IN: tmax # byref(self.error_status_c), # OUT: ErrStat_C @@ -377,15 +420,14 @@ def aerodyn_inflow_init(self, AD_input_string_array, IfW_input_string_array): # ) # # self.check_error() - # #FIXME: anything coming out that needs handling/passing? - # aerodyn_inflow_calcOutput ------------------------------------------------------------------------------------------------------------ - def aerodyn_inflow_calcOutput(self, time, hubPos, hubOrient, hubVel, hubAcc, \ + # adi_setrotormotion ------------------------------------------------------------------------------------------------------------ + def adi_setrotormotion(self, iturb, \ + hubPos, hubOrient, hubVel, hubAcc, \ nacPos, nacOrient, nacVel, nacAcc, \ rootPos, rootOrient, rootVel, rootAcc, \ - meshPos, meshOrient, meshVel, meshAcc, \ - meshFrcMom, outputChannelValues): + meshPos, meshOrient, meshVel, meshAcc): # Check input motion info self.check_input_motions_hubNac(hubPos,hubOrient,hubVel,hubAcc,'hub') @@ -403,26 +445,19 @@ def aerodyn_inflow_calcOutput(self, time, hubPos, hubOrient, hubVel, hubAcc, \ _nacAcc_c = (c_float * len(np.squeeze(nacAcc) ))(*np.squeeze(nacAcc) ) # Make a flat 1D arrays of motion info: # [x2,y1,z1, x2,y2,z2 ...] - _rootPos_flat_c = self.flatPosArr( self._initNumBlades,self.numBlades,rootPos, time,'MeshPos') - _rootOrient_flat_c = self.flatOrientArr(self._initNumBlades,self.numBlades,rootOrient,time,'MeshOrient') - _rootVel_flat_c = self.flatVelAccArr(self._initNumBlades,self.numBlades,rootVel, time,'MeshVel') - _rootAcc_flat_c = self.flatVelAccArr(self._initNumBlades,self.numBlades,rootAcc, time,'MeshAcc') + _rootPos_flat_c = self.flatPosArr( self._initNumBlades,self.numBlades,rootPos, 'MeshPos') + _rootOrient_flat_c = self.flatOrientArr(self._initNumBlades,self.numBlades,rootOrient,'MeshOrient') + _rootVel_flat_c = self.flatVelAccArr(self._initNumBlades,self.numBlades,rootVel, 'MeshVel') + _rootAcc_flat_c = self.flatVelAccArr(self._initNumBlades,self.numBlades,rootAcc, 'MeshAcc') # Make a flat 1D arrays of motion info: # [x2,y1,z1, x2,y2,z2 ...] - _meshPos_flat_c = self.flatPosArr( self._initNumMeshPts,self.numMeshPts,meshPos, time,'MeshPos') - _meshOrient_flat_c = self.flatOrientArr(self._initNumMeshPts,self.numMeshPts,meshOrient,time,'MeshOrient') - _meshVel_flat_c = self.flatVelAccArr(self._initNumMeshPts,self.numMeshPts,meshVel, time,'MeshVel') - _meshAcc_flat_c = self.flatVelAccArr(self._initNumMeshPts,self.numMeshPts,meshAcc, time,'MeshAcc') - - # Resulting Forces/moments -- [Fx1,Fy1,Fz1,Mx1,My1,Mz1, Fx2,Fy2,Fz2,Mx2,My2,Mz2 ...] - _meshFrc_flat_c = (c_float * (6 * self.numMeshPts))(0.0,) - - # Set up output channels - outputChannelValues_c = (c_float * self.numChannels)(0.0,) + _meshPos_flat_c = self.flatPosArr( self._initNumMeshPts,self.numMeshPts,meshPos, 'MeshPos') + _meshOrient_flat_c = self.flatOrientArr(self._initNumMeshPts,self.numMeshPts,meshOrient,'MeshOrient') + _meshVel_flat_c = self.flatVelAccArr(self._initNumMeshPts,self.numMeshPts,meshVel, 'MeshVel') + _meshAcc_flat_c = self.flatVelAccArr(self._initNumMeshPts,self.numMeshPts,meshAcc, 'MeshAcc') - # Run AeroDyn_Inflow_C_CalcOutput - self.AeroDyn_Inflow_C_CalcOutput( - byref(c_double(time)), # IN: time at which to calculate output forces + self.ADI_C_SetRotorMotion( + c_int(iturb), # IN: iturb -- current turbine number _hubPos_c, # IN: hub positions _hubOrient_c, # IN: hub orientations _hubVel_c, # IN: hub velocity [TVx,TVy,TVz,RVx,RVy,RVz] @@ -440,14 +475,29 @@ def aerodyn_inflow_calcOutput(self, time, hubPos, hubOrient, hubVel, hubAcc, \ _meshOrient_flat_c, # IN: Orientations (DCM) _meshVel_flat_c, # IN: velocities at desired positions _meshAcc_flat_c, # IN: accelerations at desired positions - _meshFrc_flat_c, # OUT: resulting forces/moments array - outputChannelValues_c, # OUT: output channel values as described in input file byref(self.error_status_c), # OUT: ErrStat_C self.error_message_c # OUT: ErrMsg_C ) self.check_error() + + # adi_calcOutput ------------------------------------------------------------------------------------------------------------ + def adi_getrotorloads(self, iturb, meshFrcMom): + # Resulting Forces/moments -- [Fx1,Fy1,Fz1,Mx1,My1,Mz1, Fx2,Fy2,Fz2,Mx2,My2,Mz2 ...] + _meshFrc_flat_c = (c_float * (6 * self.numMeshPts))(0.0,) + + # Run ADI_C_GetRotorLoads + self.ADI_C_GetRotorLoads( + c_int(iturb), # IN: iturb -- current turbine number + byref(c_int(self.numMeshPts)), # IN: number of attachment points expected (where motions are transferred into HD) + _meshFrc_flat_c, # OUT: resulting forces/moments array + byref(self.error_status_c), # OUT: ErrStat_C + self.error_message_c # OUT: ErrMsg_C + ) + + self.check_error() + ## Reshape Force/Moment into [N,6] count = 0 for j in range(0,self.numMeshPts): @@ -458,81 +508,47 @@ def aerodyn_inflow_calcOutput(self, time, hubPos, hubOrient, hubVel, hubAcc, \ meshFrcMom[j,4] = _meshFrc_flat_c[count+4] meshFrcMom[j,5] = _meshFrc_flat_c[count+5] count = count + 6 - - # Convert output channel values back into python - for k in range(0,self.numChannels): - outputChannelValues[k] = float(outputChannelValues_c[k]) - # aerodyn_inflow_updateStates ------------------------------------------------------------------------------------------------------------ - def aerodyn_inflow_updateStates(self, time, timeNext, \ - hubPos, hubOrient, hubVel, hubAcc, \ - nacPos, nacOrient, nacVel, nacAcc, \ - rootPos, rootOrient, rootVel, rootAcc, \ - meshPos, meshOrient, meshVel, meshAcc): - # Check input motion info - self.check_input_motions_hubNac(hubPos,hubOrient,hubVel,hubAcc,'hub') - self.check_input_motions_hubNac(nacPos,nacOrient,nacVel,nacAcc,'nacelle') - self.check_input_motions_root(rootPos,rootOrient,rootVel,rootAcc) - self.check_input_motions_mesh(meshPos,meshOrient,meshVel,meshAcc) + # adi_calcOutput ------------------------------------------------------------------------------------------------------------ + def adi_calcOutput(self, time, outputChannelValues): - _hubPos_c = (c_float * len(np.squeeze(hubPos) ))(*np.squeeze(hubPos) ) - _hubOrient_c = (c_double * len(np.squeeze(hubOrient)))(*np.squeeze(hubOrient)) - _hubVel_c = (c_float * len(np.squeeze(hubVel) ))(*np.squeeze(hubVel) ) - _hubAcc_c = (c_float * len(np.squeeze(hubAcc) ))(*np.squeeze(hubAcc) ) - _nacPos_c = (c_float * len(np.squeeze(nacPos) ))(*np.squeeze(nacPos) ) - _nacOrient_c = (c_double * len(np.squeeze(nacOrient)))(*np.squeeze(nacOrient)) - _nacVel_c = (c_float * len(np.squeeze(nacVel) ))(*np.squeeze(nacVel) ) - _nacAcc_c = (c_float * len(np.squeeze(nacAcc) ))(*np.squeeze(nacAcc) ) - # Make a flat 1D arrays of motion info: - # [x2,y1,z1, x2,y2,z2 ...] - _rootPos_flat_c = self.flatPosArr( self._initNumBlades,self.numBlades,rootPos, time,'MeshPos') - _rootOrient_flat_c = self.flatOrientArr(self._initNumBlades,self.numBlades,rootOrient,time,'MeshOrient') - _rootVel_flat_c = self.flatVelAccArr(self._initNumBlades,self.numBlades,rootVel, time,'MeshVel') - _rootAcc_flat_c = self.flatVelAccArr(self._initNumBlades,self.numBlades,rootAcc, time,'MeshAcc') - # Make a flat 1D arrays of motion info: - # [x2,y1,z1, x2,y2,z2 ...] - _meshPos_flat_c = self.flatPosArr( self._initNumMeshPts,self.numMeshPts,meshPos, time,'MeshPos') - _meshOrient_flat_c = self.flatOrientArr(self._initNumMeshPts,self.numMeshPts,meshOrient,time,'MeshOrient') - _meshVel_flat_c = self.flatVelAccArr(self._initNumMeshPts,self.numMeshPts,meshVel, time,'MeshVel') - _meshAcc_flat_c = self.flatVelAccArr(self._initNumMeshPts,self.numMeshPts,meshAcc, time,'MeshAcc') + # Set up output channels + outputChannelValues_c = (c_float * self.numChannels)(0.0,) - # Resulting Forces/moments -- [Fx1,Fy1,Fz1,Mx1,My1,Mz1, Fx2,Fy2,Fz2,Mx2,My2,Mz2 ...] - _meshFrc_flat_c = (c_float * (6 * self.numMeshPts))(0.0,) + # Run ADI_C_CalcOutput + self.ADI_C_CalcOutput( + byref(c_double(time)), # IN: time at which to calculate output forces + outputChannelValues_c, # OUT: output channel values as described in input file + byref(self.error_status_c), # OUT: ErrStat_C + self.error_message_c # OUT: ErrMsg_C + ) + + self.check_error() + + # Convert output channel values back into python + for k in range(0,self.numChannels): + outputChannelValues[k] = float(outputChannelValues_c[k]) + + # adi_updateStates ------------------------------------------------------------------------------------------------------------ + def adi_updateStates(self, time, timeNext): # Run AeroDyn_Inflow_UpdateStates_c - self.AeroDyn_Inflow_C_UpdateStates( + self.ADI_C_UpdateStates( byref(c_double(time)), # IN: time at which to calculate output forces byref(c_double(timeNext)), # IN: time T+dt we are stepping to - _hubPos_c, # IN: hub positions - _hubOrient_c, # IN: hub orientations - _hubVel_c, # IN: hub velocity [TVx,TVy,TVz,RVx,RVy,RVz] - _hubAcc_c, # IN: hub acclerations [TAx,TAy,TAz,RAx,RAy,RAz] - _nacPos_c, # IN: nac positions - _nacOrient_c, # IN: nac orientations - _nacVel_c, # IN: nac velocity [TVx,TVy,TVz,RVx,RVy,RVz] - _nacAcc_c, # IN: nac acclerations [TAx,TAy,TAz,RAx,RAy,RAz] - _rootPos_flat_c, # IN: positions - _rootOrient_flat_c, # IN: Orientations (DCM) - _rootVel_flat_c, # IN: velocities at desired positions - _rootAcc_flat_c, # IN: accelerations at desired positions - byref(c_int(self.numMeshPts)), # IN: number of attachment points expected (where motions are transferred into HD) - _meshPos_flat_c, # IN: positions - _meshOrient_flat_c, # IN: Orientations (DCM) - _meshVel_flat_c, # IN: velocities at desired positions - _meshAcc_flat_c, # IN: accelerations at desired positions byref(self.error_status_c), # OUT: ErrStat_C self.error_message_c # OUT: ErrMsg_C ) self.check_error() - # aerodyn_inflow_end ------------------------------------------------------------------------------------------------------------ - def aerodyn_inflow_end(self): + # adi_end ------------------------------------------------------------------------------------------------------------ + def adi_end(self): if not self.ended: self.ended = True - # Run AeroDyn_Inflow_C_End - self.AeroDyn_Inflow_C_End( + # Run ADI_C_End + self.ADI_C_End( byref(self.error_status_c), self.error_message_c ) @@ -547,14 +563,14 @@ def check_error(self): print(f"AeroDyn/InflowWind error status: {self.error_levels[self.error_status_c.value]}: {self.error_message_c.value.decode('ascii')}") else: print(f"AeroDyn/InflowWind error status: {self.error_levels[self.error_status_c.value]}: {self.error_message_c.value.decode('ascii')}") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/InflowWind terminated prematurely.") - def flatPosArr(self,initNumMeshPts,numPts,MeshPosArr,time,name): + def flatPosArr(self,initNumMeshPts,numPts,MeshPosArr,name): if initNumMeshPts != numPts: - print(f"At time {time}, the number of {name} points changed from initial value of {initNumMeshPts}. This is not permitted during the simulation.") - self.aerodyn_inflow_end() + print(f"The number of {name} points changed from initial value of {initNumMeshPts}. This is not permitted during the simulation.") + self.adi_end() raise Exception("\nError in calling AeroDyn/InflowWind library.") meshPos_flat = [pp for p in MeshPosArr for pp in p] meshPos_flat_c = (c_float * (3 * numPts))(0.0,) @@ -563,10 +579,10 @@ def flatPosArr(self,initNumMeshPts,numPts,MeshPosArr,time,name): return meshPos_flat_c - def flatOrientArr(self,initNumMeshPts,numPts,MeshOrientArr,time,name): + def flatOrientArr(self,initNumMeshPts,numPts,MeshOrientArr,name): if initNumMeshPts != numPts: - print(f"At time {time}, the number of {name} points changed from initial value of {initNumMeshPts}. This is not permitted during the simulation.") - self.aerodyn_inflow_end() + print(f"The number of {name} points changed from initial value of {initNumMeshPts}. This is not permitted during the simulation.") + self.adi_end() raise Exception("\nError in calling AeroDyn/InflowWind library.") meshOrient_flat = [pp for p in MeshOrientArr for pp in p] meshOrient_flat_c = (c_double * (9 * numPts))(0.0,) @@ -575,10 +591,10 @@ def flatOrientArr(self,initNumMeshPts,numPts,MeshOrientArr,time,name): return meshOrient_flat_c - def flatVelAccArr(self,initNumMeshPts,numPts,MeshArr,time,name): + def flatVelAccArr(self,initNumMeshPts,numPts,MeshArr,name): if initNumMeshPts != numPts: - print(f"At time {time}, the number of {name} points changed from initial value of {initNumMeshPts}. This is not permitted during the simulation.") - self.aerodyn_inflow_end() + print(f"The number of {name} points changed from initial value of {initNumMeshPts}. This is not permitted during the simulation.") + self.adi_end() raise Exception("\nError in calling AeroDyn/InflowWind library.") # Velocity -- [Vx2,Vy1,Vz1,RVx1,RVy1,RVz1, Vx2,Vy2,Vz2,RVx2,RVy2,RVz2 ...] meshVel_flat = [pp for p in MeshArr for pp in p] @@ -612,39 +628,39 @@ def check_init_hubroot(self): #print(" size 0 ", self.initNacelleOrient.shape[0]) if self.numBlades < 1: print("No blades. Set numBlades to number of AD blades in the model") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if self.initRootPos.shape[1] != 3: print("Expecting a Nx3 array of blade root positions (initRootPos) with second index for [x,y,z]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if self.initRootPos.shape[0] != self.numBlades: print("Expecting a Nx3 array of blade root positions (initRootPos) with first index for blade number") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if self.initRootOrient.shape[1] != 9: print("Expecting a Nx9 array of blade root orientations as DCMs (initRootOrient) with second index for [r11,r12,r13,r21,r22,r23,r31,r32,r33]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if self.initRootOrient.shape[0] != self.numBlades: print("Expecting a Nx3 array of blade root orientations (initRootOrient) with first index for blade number") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if np.squeeze(self.initHubPos.ndim) > 1 or self.initHubPos.shape[0] != 3: print("Expecting a 3 element array for initHubPos [x,y,z]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if np.squeeze(self.initHubOrient.ndim) > 1 or self.initHubOrient.shape[0] != 9: print("Expecting a 9 element array for initHubOrient DCM [r11,r12,r13,r21,r22,r23,r31,r32,r33]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if np.squeeze(self.initNacellePos.ndim) > 1 or self.initNacellePos.shape[0] != 3: print("Expecting a 3 element array for initNacellePos [x,y,z]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if np.squeeze(self.initNacelleOrient.ndim) > 1 or self.initNacelleOrient.shape[0] != 9: print("Expecting a 9 element array for initNacelleOrient DCM [r11,r12,r13,r21,r22,r23,r31,r32,r33]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") @@ -660,23 +676,23 @@ def check_init_mesh(self): # Verify that the shape of initMeshPos is correct if self.initMeshPos.shape[0] != self.initMeshOrient.shape[0]: print("Different number of meshs in inital position and orientation arrays") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if self.initMeshPos.shape[1] != 3: print("Expecting a Nx3 array of initial mesh positions (initMeshPos) with second index for [x,y,z]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if self.initMeshPos.shape[0] != self.numMeshPts: print("Expecting a Nx3 array of initial mesh positions (initMeshPos) with first index for mesh number.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if self.initMeshOrient.shape[1] != 9: print("Expecting a Nx9 array of initial mesh orientations as DCMs (initMeshOrient) with second index for [r11,r12,r13,r21,r22,r23,r31,r32,r33]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") if self.initMeshOrient.shape[0] != self.numMeshPts: print("Expecting a Nx3 array of initial mesh orientations (initMeshOrient) with first index for mesh number.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn terminated prematurely.") @@ -684,72 +700,72 @@ def check_input_motions_hubNac(self,nodePos,nodeOrient,nodeVel,nodeAcc,_name): # Verify that the shape of positions array is correct if nodePos.size != 3: print("Expecting a Nx3 array of "+_name+" positions with second index for [x,y,z]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") # Verify that the shape of orientations array is correct if nodeOrient.size != 9: print("Expecting a Nx9 array of "+_name+" orientations with second index for [r11,r12,r13,r21,r22,r23,r31,r32,r33]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") # Verify that the shape of velocities array is correct if nodeVel.size != 6: print("Expecting a Nx6 array of "+_name+" velocities with second index for [x,y,z,Rx,Ry,Rz]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") # Verify that the shape of accelerations array is correct if nodeAcc.size != 6: print("Expecting a Nx6 array of "+_name+" accelerations with second index for [x,y,z,Rx,Ry,Rz]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") def check_input_motions_root(self,rootPos,rootOrient,rootVel,rootAcc): # make sure number of roots didn't change for some reason if self._initNumBlades != self.numBlades: print(f"At time {time}, the number of root points changed from initial value of {self._initNumBlades}. This is not permitted during the simulation.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nError in calling AeroDyn/AeroDyn library.") # Verify that the shape of positions array is correct if rootPos.shape[1] != 3: print("Expecting a Nx3 array of root positions (rootOrient) with second index for [x,y,z]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") if rootPos.shape[0] != self.numBlades: print("Expecting a Nx3 array of root positions (rootOrient) with first index for root number.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") # Verify that the shape of orientations array is correct if rootOrient.shape[1] != 9: print("Expecting a Nx9 array of root orientations (rootPos) with second index for [r11,r12,r13,r21,r22,r23,r31,r32,r33]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") if rootOrient.shape[0] != self.numBlades: print("Expecting a Nx9 array of root orientations (rootPos) with first index for root number.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") # Verify that the shape of velocities array is correct if rootVel.shape[1] != 6: print("Expecting a Nx6 array of root velocities (rootVel) with second index for [x,y,z,Rx,Ry,Rz]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") if rootVel.shape[0] != self.numBlades: print("Expecting a Nx6 array of root velocities (rootVel) with first index for root number.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") # Verify that the shape of accelerations array is correct if rootAcc.shape[1] != 6: print("Expecting a Nx6 array of root accelerations (rootAcc) with second index for [x,y,z,Rx,Ry,Rz]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") if rootAcc.shape[0] != self.numBlades: print("Expecting a Nx6 array of root accelerations (rootAcc) with first index for root number.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") @@ -757,47 +773,47 @@ def check_input_motions_mesh(self,meshPos,meshOrient,meshVel,meshAcc): # make sure number of meshs didn't change for some reason if self._initNumMeshPts != self.numMeshPts: print(f"At time {time}, the number of mesh points changed from initial value of {self._initNumMeshPts}. This is not permitted during the simulation.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nError in calling AeroDyn/AeroDyn library.") # Verify that the shape of positions array is correct if meshPos.shape[1] != 3: print("Expecting a Nx3 array of mesh positions (meshOrient) with second index for [x,y,z]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") if meshPos.shape[0] != self.numMeshPts: print("Expecting a Nx3 array of mesh positions (meshOrient) with first index for mesh number.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") # Verify that the shape of orientations array is correct if meshOrient.shape[1] != 9: print("Expecting a Nx9 array of mesh orientations (meshPos) with second index for [r11,r12,r13,r21,r22,r23,r31,r32,r33]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") if meshOrient.shape[0] != self.numMeshPts: print("Expecting a Nx9 array of mesh orientations (meshPos) with first index for mesh number.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") # Verify that the shape of velocities array is correct if meshVel.shape[1] != 6: print("Expecting a Nx6 array of mesh velocities (meshVel) with second index for [x,y,z,Rx,Ry,Rz]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") if meshVel.shape[0] != self.numMeshPts: print("Expecting a Nx6 array of mesh velocities (meshVel) with first index for mesh number.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") # Verify that the shape of accelerations array is correct if meshAcc.shape[1] != 6: print("Expecting a Nx6 array of mesh accelerations (meshAcc) with second index for [x,y,z,Rx,Ry,Rz]") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") if meshAcc.shape[0] != self.numMeshPts: print("Expecting a Nx6 array of mesh accelerations (meshAcc) with first index for mesh number.") - self.aerodyn_inflow_end() + self.adi_end() raise Exception("\nAeroDyn/AeroDyn terminated prematurely.") @@ -826,7 +842,8 @@ def output_channel_units(self): # correctly, this will be identical to the corresponding values in the # AeroDyn/InflowWind output channels. -#FIXME: this is incorrect +#FIXME: this may not output everything in the interface (updates have been made +# since writing this, but this routine was not updated accordingly class DriverDbg(): """ This is only for debugging purposes only. The input motions and resulting @@ -841,8 +858,8 @@ def __init__(self,filename,numMeshPts): # write file header t_string=datetime.datetime.now() dt_string=datetime.date.today() - self.DbgFile.write(f"## This file was generated by aerodyn_inflow_c_lib on {dt_string.strftime('%b-%d-%Y')} at {t_string.strftime('%H:%M:%S')}\n") - self.DbgFile.write(f"## This file contains the resulting forces/moments at each of {self.numMeshPts} mesh points passed into the aerodyn_inflow_c_lib\n") + self.DbgFile.write(f"## This file was generated by adi_c_lib on {dt_string.strftime('%b-%d-%Y')} at {t_string.strftime('%H:%M:%S')}\n") + self.DbgFile.write(f"## This file contains the resulting forces/moments at each of {self.numMeshPts} mesh points passed into the adi_c_lib\n") self.DbgFile.write("#\n") self.DbgFile.write("#\n") self.DbgFile.write("#\n") @@ -927,7 +944,7 @@ def end(self): # Helper class for writing channels to file. # for the regression testing to mirror the output from the AD15 and InfowWind # from an OpenFAST simulation. This may also have value for debugging -# interfacing to the AeroDyn_Inflow_C_Binding library +# interfacing to the ADI_C_Binding library class WriteOutChans(): """ diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index c622f0d12e..748e6ab26d 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -1120,7 +1120,7 @@ subroutine Init_u( u, p, p_AD, InputFileData, MHK, WtrDpth, InitInp, errStat, er if (errStat >= AbortErrLev) return ! set node initial position/orientation - position = 0.0_ReKi + position = InitInp%originInit do j=1,p%NumTwrNds IF ( MHK == MHK_FixedBottom ) THEN position(3) = InputFileData%TwrElev(j) - WtrDpth diff --git a/modules/aerodyn/src/AeroDyn_Driver_Registry.txt b/modules/aerodyn/src/AeroDyn_Driver_Registry.txt index b216fc47f3..8cd3d8dd7d 100644 --- a/modules/aerodyn/src/AeroDyn_Driver_Registry.txt +++ b/modules/aerodyn/src/AeroDyn_Driver_Registry.txt @@ -138,7 +138,7 @@ typedef ^ ^ ReKi Patm typedef ^ ^ ReKi Pvap - - - "Vapour pressure of working fluid" "Pa" typedef ^ ^ ReKi WtrDpth - - - "Water depth" "m" typedef ^ ^ ReKi MSL2SWL - - - "Offset between still-water level and mean sea level" "m" -typedef ^ ^ IntKi numTurbines - - - "number of blades on turbine" "-" +typedef ^ ^ IntKi numTurbines - -9999 - "number of turbine rotors" "-" typedef ^ ^ WTData WT : - - "Wind turbine data for driver" "-" typedef ^ ^ DbKi dT - - - "time increment" "s" typedef ^ ^ DbKi tMax - - - "time increment" "s" diff --git a/modules/aerodyn/src/AeroDyn_Driver_Types.f90 b/modules/aerodyn/src/AeroDyn_Driver_Types.f90 index 5cfa46275a..71835dac10 100644 --- a/modules/aerodyn/src/AeroDyn_Driver_Types.f90 +++ b/modules/aerodyn/src/AeroDyn_Driver_Types.f90 @@ -172,7 +172,7 @@ MODULE AeroDyn_Driver_Types REAL(ReKi) :: Pvap = 0.0_ReKi !< Vapour pressure of working fluid [Pa] REAL(ReKi) :: WtrDpth = 0.0_ReKi !< Water depth [m] REAL(ReKi) :: MSL2SWL = 0.0_ReKi !< Offset between still-water level and mean sea level [m] - INTEGER(IntKi) :: numTurbines = 0_IntKi !< number of blades on turbine [-] + INTEGER(IntKi) :: numTurbines = -9999 !< number of turbine rotors [-] TYPE(WTData) , DIMENSION(:), ALLOCATABLE :: WT !< Wind turbine data for driver [-] REAL(DbKi) :: dT = 0.0_R8Ki !< time increment [s] REAL(DbKi) :: tMax = 0.0_R8Ki !< time increment [s] diff --git a/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90 b/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90 index 65e970e32f..0b81ec342c 100644 --- a/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90 +++ b/modules/aerodyn/src/AeroDyn_Inflow_C_Binding.f90 @@ -30,11 +30,15 @@ MODULE AeroDyn_Inflow_C_BINDING IMPLICIT NONE - PUBLIC :: AeroDyn_Inflow_C_Init - !PUBLIC :: AeroDyn_Inflow_C_ReInit - PUBLIC :: AeroDyn_Inflow_C_CalcOutput - PUBLIC :: AeroDyn_Inflow_C_UpdateStates - PUBLIC :: AeroDyn_Inflow_C_End + PUBLIC :: ADI_C_Init + !PUBLIC :: ADI_C_ReInit + PUBLIC :: ADI_C_CalcOutput + PUBLIC :: ADI_C_UpdateStates + PUBLIC :: ADI_C_End + PUBLIC :: ADI_C_PreInit ! Initial call to setup number of turbines + PUBLIC :: ADI_C_SetupRotor ! Initial node positions etc for a rotor + PUBLIC :: ADI_C_SetRotorMotion ! Set motions for a given rotor + PUBLIC :: ADI_C_GetRotorLoads ! Retrieve loads for a given rotor !------------------------------------------------------------------------------------ ! Version info for display @@ -42,13 +46,13 @@ MODULE AeroDyn_Inflow_C_BINDING !------------------------------------------------------------------------------------ - ! Debugging: debugverbose + ! Debugging: debugverbose -- passed at PreInit ! 0 - none ! 1 - some summary info ! 2 - above + all position/orientation info ! 3 - above + input files (if direct passed) ! 4 - above + meshes - integer(IntKi), parameter :: debugverbose = 0 + integer(IntKi) :: debugverbose = 0 !------------------------------------------------------------------------------------ ! Error handling @@ -90,6 +94,7 @@ MODULE AeroDyn_Inflow_C_BINDING type(ADI_Data) :: ADI type(ADI_InitInputType) :: InitInp !< Initialization data type(ADI_InitOutputType) :: InitOutData !< Initial output data -- Names, units, and version info. + type(ADI_InputType) :: ADI_u !< ADI inputs -- set by AD_SetInputMotion. Copied as needed (necessary for correction steps) !------------------------------ ! Simulation data type(Dvr_SimData) :: Sim !< data about the simulation @@ -142,22 +147,22 @@ MODULE AeroDyn_Inflow_C_BINDING ! - 1 point -- rigid floating body assumption ! - N points -- flexible structure (either floating or fixed bottom) logical :: TransposeDCM !< Transpose DCMs as passed in -- test the vtk outputs to see if needed - integer(IntKi) :: NumMeshPts ! Number of mesh points we are interfacing motions/loads to/from AD - type(MeshType) :: BldPtMotionMesh ! mesh for motions of external nodes - type(MeshType) :: BldPtLoadMesh ! mesh for loads for external nodes - type(MeshType) :: BldPtLoadMesh_tmp ! mesh for loads for external nodes -- temporary -! type(MeshType) :: NacMotionMesh ! mesh for motion of nacelle -- TODO: add this mesh for nacelle load transfers -! type(MeshType) :: NacLoadMesh ! mesh for loads for nacelle loads -- TODO: add this mesh for nacelle load transfers + integer(IntKi), allocatable :: NumMeshPts(:) ! Number of mesh points we are interfacing motions/loads to/from AD for each rotor + type(MeshType), allocatable :: BldPtMotionMesh(:) ! mesh for motions of external nodes + type(MeshType), allocatable :: BldPtLoadMesh(:) ! mesh for loads for external nodes + type(MeshType), allocatable :: BldPtLoadMesh_tmp(:) ! mesh for loads for external nodes -- temporary +! type(MeshType), allocatable :: NacMotionMesh(:) ! mesh for motion of nacelle -- TODO: add this mesh for nacelle load transfers +! type(MeshType), allocatable :: NacLoadMesh(:) ! mesh for loads for nacelle loads -- TODO: add this mesh for nacelle load transfers !------------------------------ ! Mesh mapping: motions ! The mapping of motions from the nodes passed in to the corresponding AD meshes - type(MeshMapType), allocatable :: Map_BldPtMotion_2_AD_Blade(:) ! Mesh mapping between input motion mesh for blade -! type(MeshMapType) :: Map_AD_Nac_2_NacPtLoad ! Mesh mapping between input motion mesh for nacelle + type(MeshMapType), allocatable :: Map_BldPtMotion_2_AD_Blade(:,:) ! Mesh mapping between input motion mesh for blade + type(MeshMapType), allocatable :: Map_AD_Nac_2_NacPtLoad(:) ! Mesh mapping between input motion mesh for nacelle !------------------------------ ! Mesh mapping: loads ! The mapping of loads from the AD meshes to the corresponding external nodes - type(MeshMapType), allocatable :: Map_AD_BldLoad_P_2_BldPtLoad(:) ! Mesh mapping between AD output blade line2 load to BldPtLoad for return -! type(MeshMapType) :: Map_NacPtMotion_2_AD_Nac ! Mesh mapping between AD output nacelle pt load to NacLoad for return + type(MeshMapType), allocatable :: Map_AD_BldLoad_P_2_BldPtLoad(:,:) ! Mesh mapping between AD output blade line2 load to BldPtLoad for return +! type(MeshMapType) :: Map_NacPtMotion_2_AD_Nac(:) ! Mesh mapping between AD output nacelle pt load to NacLoad for return ! Motions input (so we don't have to reallocate all the time real(ReKi), allocatable :: tmpBldPtMeshPos(:,:) ! temp array. Probably don't need this, but makes conversion from C clearer. real(ReKi), allocatable :: tmpBldPtMeshOri(:,:,:) ! temp array. Probably don't need this, but makes conversion from C clearer. @@ -167,6 +172,12 @@ MODULE AeroDyn_Inflow_C_BINDING !------------------------------------------------------------------------------------ + ! NOTE on turbine origin + ! The turbine origin is set by TurbOrigin_C during the ADI_C_SetupRotor routine. This is the tower base location. All + ! blade, tower, nacelle, and hub coordinates are relative to this location. Since AD15 and IfW use absolute positioning, + ! the reference positions for the blades, tower, nacelle, and hub are set by the values passed into ADI_C_SetupRotor + + ! TurbOrigin_C (stored as Sim%WT(iWT)%OriginInit). When the mesh and other points are passed in, they are relative to + ! their respective rotor origin. CONTAINS @@ -187,37 +198,173 @@ subroutine SetErr(ErrStat, ErrMsg, ErrStat_C, ErrMsg_C) else ErrMsg_C = TRANSFER( trim(ErrMsg)//C_NULL_CHAR, ErrMsg_C ) endif + if (ErrStat /= ErrID_None) call WrScr(NewLine//'ADI_C_Binding: '//trim(ErrMsg)//NewLine) end subroutine SetErr +!=============================================================================================================== +!--------------------------------------------- AeroDyn PreInit ------------------------------------------------- +!=============================================================================================================== +!> Allocate all the arrays for data storage for all turbine rotors +subroutine ADI_C_PreInit(NumTurbines_C,TransposeDCM_in,debuglevel,ErrStat_C,ErrMsg_C) BIND (C, NAME='ADI_C_PreInit') + implicit none +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_PreInit +!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_PreInit +#endif + integer(c_int), intent(in ) :: NumTurbines_C + integer(c_int), intent(in ) :: TransposeDCM_in !< Transpose DCMs as they are passed i + integer(c_int), intent(in ) :: debuglevel + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + + ! Local variables + integer(IntKi) :: iWT !< current turbine + integer :: ErrStat !< aggregated error status + character(ErrMsgLen) :: ErrMsg !< aggregated error message + integer :: ErrStat2 !< temporary error status from a call + character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call + character(*), parameter :: RoutineName = 'ADI_C_PreInit' !< for error handling + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + CALL NWTC_Init( ProgNameIn=version%Name ) + CALL DispCopyrightLicense( version%Name ) + CALL DispCompileRuntimeInfo( version%Name ) + + ! interface debugging + debugverbose = int(debuglevel,IntKi) + ! if non-zero, show all passed data here. Then check valid values + if (debugverbose /= 0_IntKi) then + call WrScr(" Interface debugging level "//trim(Num2Lstr(debugverbose))//" requested.") + call ShowPassedData() + endif + + ! check valid debug level + if (debugverbose < 0_IntKi .or. debugverbose > 4_IntKi) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Interface debug level must be between 0 and 4"//NewLine// & + " 0 - none"//NewLine// & + " 1 - some summary info and variables passed through interface"//NewLine// & + " 2 - above + all position/orientation info"//NewLine// & + " 3 - above + input files (if direct passed)"//NewLine// & + " 4 - above + meshes" + if (Failed()) return; + endif + + + ! Set number of turbines + Sim%NumTurbines = int(NumTurbines_C,IntKi) + + if (Sim%NumTurbines < 1_IntKi .or. Sim%NumTurbines > 9_IntKi) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = 'AeroDyn_Inflow simulates between 1 and 9 turbines, but '//trim(Num2LStr(Sim%NumTurbines))//' was specified' + if (Failed()) return; + endif + + ! Flag to transpose DCMs as they are passed in + TransposeDCM = TransposeDCM_in==1_c_int + + ! Allocate arrays and meshes for the number of turbines + if (allocated(InitInp%AD%rotors)) deallocate(InitInp%AD%rotors) + allocate(InitInp%AD%rotors(Sim%NumTurbines),stat=errStat2); if (Failed0('rotors')) return + + ! Allocate data storage for turbine info + if (allocated(Sim%WT)) deallocate(Sim%WT) + allocate(Sim%WT(Sim%NumTurbines),stat=errStat2); if (Failed0('wind turbines')) return + do iWT=1,Sim%NumTurbines + Sim%WT(iWT)%NumBlades = -999 + enddo + + ! storage for numnber of meshpoints + if (allocated(NumMeshPts)) deallocate(NumMeshPts) + allocate(NumMeshPts(Sim%NumTurbines),stat=errStat2); if (Failed0('NumMeshPts')) return + NumMeshPts = -999 + + ! allocate meshes and meshmappings + if (allocated(BldPtMotionMesh )) deallocate(BldPtMotionMesh ) + if (allocated(BldPtLoadMesh )) deallocate(BldPtLoadMesh ) + if (allocated(BldPtLoadMesh_tmp)) deallocate(BldPtLoadMesh_tmp) + !if (allocated(NacMotionMesh )) deallocate(NacMotionMesh ) + !if (allocated(NacLoadMesh )) deallocate(NacLoadMesh ) + if (allocated(Map_BldPtMotion_2_AD_Blade )) deallocate(Map_BldPtMotion_2_AD_Blade ) + if (allocated(Map_AD_BldLoad_P_2_BldPtLoad)) deallocate(Map_AD_BldLoad_P_2_BldPtLoad) + !if (allocated(Map_NacPtMotion_2_AD_Nac )) deallocate(Map_NacPtMotion_2_AD_Nac ) + allocate(BldPtMotionMesh( Sim%NumTurbines), STAT=ErrStat2); if (Failed0('BldPtMotionMesh' )) return + allocate(BldPtLoadMesh( Sim%NumTurbines), STAT=ErrStat2); if (Failed0('BldPtLoadMesh' )) return + allocate(BldPtLoadMesh_tmp(Sim%NumTurbines), STAT=ErrStat2); if (Failed0('BldPtLoadMesh_tmp')) return + !allocate(NacMotionMesh( Sim%NumTurbines), STAT=ErrStat2); if (Failed0('NacMotionMesh' )) return + !allocate(NacLoadMesh( Sim%NumTurbines), STAT=ErrStat2); if (Failed0('NacLoadMesh' )) return + !allocate(Map_NacPtMotion_2_AD_Nac(Sim%NumTurbines),STAT=ErrStat2); if (Failed0('Map_AD_BldLoad_P_2_BldPtLoad')) return + + + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + +contains + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) then + call ClearTmpStorage() + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + endif + end function Failed + + ! check for failed where /= 0 is fatal + logical function Failed0(txt) + character(*), intent(in) :: txt + if (errStat /= 0) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Could not allocate "//trim(txt) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + endif + Failed0 = ErrStat >= AbortErrLev + if(Failed0) call ClearTmpStorage() + end function Failed0 + + !> This subroutine prints out all the variables that are passed in. Use this only + !! for debugging the interface on the Fortran side. + subroutine ShowPassedData() + character(1) :: TmpFlag + call WrSCr("") + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: Variables passed in through interface") + call WrScr(" ADI_C_PreInit") + call WrScr(" --------------------------------------") + call WrScr(" NumTurbines_C "//trim(Num2LStr( NumTurbines_C )) ) + TmpFlag="F"; if (TransposeDCM_in==1_c_int) TmpFlag="T" + call WrScr(" TransposeDCM_in "//TmpFlag ) + call WrScr(" debuglevel "//trim(Num2LStr( debuglevel )) ) + call WrScr("-----------------------------------------------------------") + end subroutine ShowPassedData + +end subroutine ADI_C_PreInit + !=============================================================================================================== !--------------------------------------------- AeroDyn Init---------------------------------------------------- !=============================================================================================================== -SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinputFileStringLength_C, & +SUBROUTINE ADI_C_Init( ADinputFilePassed, ADinputFileString_C, ADinputFileStringLength_C, & IfWinputFilePassed, IfWinputFileString_C, IfWinputFileStringLength_C, OutRootName_C, & gravity_C, defFldDens_C, defKinVisc_C, defSpdSound_C, & defPatm_C, defPvap_C, WtrDpth_C, MSL2SWL_C, & - AeroProjMod_C, & InterpOrder_C, DT_C, TMax_C, & - storeHHVel, TransposeDCM_in, & + storeHHVel, & WrVTK_in, WrVTK_inType, VTKNacDim_in, VTKHubRad_in, & wrOuts_C, DT_Outs_C, & - HubPos_C, HubOri_C, & - NacPos_C, NacOri_C, & - NumBlades_C, BldRootPos_C, BldRootOri_C, & - NumMeshPts_C, InitMeshPos_C, InitMeshOri_C, & NumChannels_C, OutputChannelNames_C, OutputChannelUnits_C, & - ErrStat_C, ErrMsg_C) BIND (C, NAME='AeroDyn_Inflow_C_Init') + ErrStat_C, ErrMsg_C) BIND (C, NAME='ADI_C_Init') implicit none #ifndef IMPLICIT_DLLEXPORT -!DEC$ ATTRIBUTES DLLEXPORT :: AeroDyn_Inflow_C_Init -!GCC$ ATTRIBUTES DLLEXPORT :: AeroDyn_Inflow_C_Init +!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_Init +!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_Init #endif ! Input file info - logical(c_bool), intent(in ) :: ADinputFilePassed !< Write VTK outputs [0: none, 1: init only, 2: animation] + integer(c_int), intent(in ) :: ADinputFilePassed !< Write VTK outputs [0: none, 1: init only, 2: animation] type(c_ptr), intent(in ) :: ADinputFileString_C !< Input file as a single string with lines deliniated by C_NULL_CHAR integer(c_int), intent(in ) :: ADinputFileStringLength_C !< lenght of the input file string - logical(c_bool), intent(in ) :: IfWinputFilePassed !< Write VTK outputs [0: none, 1: init only, 2: animation] + integer(c_int), intent(in ) :: IfWinputFilePassed !< Write VTK outputs [0: none, 1: init only, 2: animation] type(c_ptr), intent(in ) :: IfWinputFileString_C !< Input file as a single string with lines deliniated by C_NULL_CHAR integer(c_int), intent(in ) :: IfWinputFileStringLength_C !< lenght of the input file string character(kind=c_char), intent(in ) :: OutRootName_C(IntfStrLen) !< Root name to use for echo files and other @@ -230,31 +377,13 @@ SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinpu real(c_float), intent(in ) :: defPvap_C !< Vapour pressure of working fluid (Pa) [used only for an MHK turbine cavitation check] real(c_float), intent(in ) :: WtrDpth_C !< Water depth (m) real(c_float), intent(in ) :: MSL2SWL_C !< Offset between still-water level and mean sea level (m) [positive upward] - ! Aero calculation method -- AeroProjMod - ! APM_BEM_NoSweepPitchTwist - 1 - "Original AeroDyn model where momentum balance is done in the WithoutSweepPitchTwist system" - ! APM_BEM_Polar - 2 - "Use staggered polar grid for momentum balance in each annulus" - ! APM_LiftingLine - 3 - "Use the blade lifting line (i.e. the structural) orientation (currently for OLAF with VAWT)" - integer(c_int), intent(in ) :: AeroProjMod_C !< Type of aerodynamic projection - ! Initial hub and blade root positions/orientations - real(c_float), intent(in ) :: HubPos_C( 3 ) !< Hub position - real(c_double), intent(in ) :: HubOri_C( 9 ) !< Hub orientation - real(c_float), intent(in ) :: NacPos_C( 3 ) !< Nacelle position - real(c_double), intent(in ) :: NacOri_C( 9 ) !< Nacelle orientation - integer(c_int), intent(in ) :: NumBlades_C !< Number of blades - real(c_float), intent(in ) :: BldRootPos_C( 3*NumBlades_C ) !< Blade root positions - real(c_double), intent(in ) :: BldRootOri_C( 9*NumBlades_C ) !< Blade root orientations - ! Initial nodes - integer(c_int), intent(in ) :: NumMeshPts_C !< Number of mesh points we are transfering motions to and output loads to - real(c_float), intent(in ) :: InitMeshPos_C( 3*NumMeshPts_C ) !< A 3xNumMeshPts_C array [x,y,z] - real(c_double), intent(in ) :: InitMeshOri_C( 9*NumMeshPts_C ) !< A 9xNumMeshPts_C array [r11,r12,r13,r21,r22,r23,r31,r32,r33] ! Interpolation integer(c_int), intent(in ) :: InterpOrder_C !< Interpolation order to use (must be 1 or 2) ! Time real(c_double), intent(in ) :: DT_C !< Timestep used with AD for stepping forward from t to t+dt. Must be constant. real(c_double), intent(in ) :: TMax_C !< Maximum time for simulation ! Flags - logical(c_bool), intent(in ) :: storeHHVel !< Store hub height time series from IfW - logical(c_bool), intent(in ) :: TransposeDCM_in !< Transpose DCMs as they are passed in + integer(c_int), intent(in ) :: storeHHVel !< Store hub height time series from IfW ! VTK integer(c_int), intent(in ) :: WrVTK_in !< Write VTK outputs [0: none, 1: init only, 2: animation] integer(c_int), intent(in ) :: WrVTK_inType !< Write VTK outputs as [1: surface, 2: lines, 3: both] @@ -280,18 +409,45 @@ SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinpu integer(IntKi) :: ErrStat2 !< temporary error status from a call character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call integer(IntKi) :: i,j,k !< generic counters - character(*), parameter :: RoutineName = 'AeroDyn_Inflow_C_Init' !< for error handling + integer(IntKi) :: iWT !< current turbine number (iterate through during setup for ADI_Init call) + integer(IntKi) :: AeroProjMod !< for checking that all turbines use the same AeroProjMod + character(*), parameter :: RoutineName = 'ADI_C_Init' !< for error handling ! Initialize error handling ErrStat = ErrID_None ErrMsg = "" + ErrStat2 = ErrID_None + ErrMsg2 = "" NumChannels_C = 0_c_int OutputChannelNames_C(:) = '' OutputChannelUnits_C(:) = '' - CALL NWTC_Init( ProgNameIn=version%Name ) - CALL DispCopyrightLicense( version%Name ) - CALL DispCompileRuntimeInfo( version%Name ) + + ! check if Pre-Init was called + if (Sim%NumTurbines < 0_IntKi) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Call ADI_C_PreInit and ADI_C_SetupRotor prior to calling ADI_C_Init" + if (Failed()) return + endif + + do iWT=1,Sim%NumTurbines + if (Sim%WT(iWT)%NumBlades < 0) call SetErrStat(ErrID_Fatal,"Rotor "//trim(Num2LStr(iWT))//" not initialized. Call ADI_C_SetupRotor prior to calling ADI_C_Init",ErrStat,ErrMsg,RoutineName) + enddo + if (Failed()) return + + + ! Check that all turbines are using the same AeroProjMod (mixing projection modes is not currently supported) + AeroProjMod = InitInp%AD%rotors(1)%AeroProjMod + do iWT = 2,Sim%NumTurbines + if(AeroProjMod /= InitInp%AD%rotors(iWT)%AeroProjMod) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Different AeroProjMod values for each turbine (set from TurbineIsHAWT flag). Check that all turbines are of the same type (HAWT or not)." + if (Failed()) return + endif + enddo + + ! Setup temporary storage arrays for simpler transfers + call SetTempStorage(ErrStat2,ErrMsg2); if (Failed()) return !-------------------------- @@ -315,7 +471,7 @@ SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinpu ! Format AD input file contents InitInp%AD%RootName = OutRootName - if (ADinputFilePassed) then + if (ADinputFilePassed==1_c_int) then InitInp%AD%UsePrimaryInputFile = .FALSE. ! Don't try to read an input -- use passed data instead (blades and AF tables not passed) InitInp%AD%InputFile = "passed_ad_file" ! not actually used call InitFileInfo(ADinputFileString, InitInp%AD%PassedPrimaryInputData, ErrStat2, ErrMsg2); if (Failed()) return @@ -332,7 +488,7 @@ SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinpu ! Format IfW input file contents ! RootName is set in ADI_Init using InitInp%RootName InitInp%RootName = OutRootName - if (IfWinputFilePassed) then + if (IfWinputFilePassed==1_c_int) then InitInp%IW_InitInp%UseInputFile = .FALSE. ! Don't try to read an input -- use passed data instead (blades and AF tables not passed) InitInp%IW_InitInp%InputFile = "passed_ifw_file" ! not actually used call InitFileInfo(IfWinputFileString, InitInp%IW_InitInp%PassedFileData, ErrStat2, ErrMsg2); if (Failed()) return @@ -351,24 +507,20 @@ SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinpu ! of the InFileInfo data structure. ! CU is the screen -- system dependent. if (debugverbose >= 3) then - if (ADinputFilePassed) call Print_FileInfo_Struct( CU, InitInp%AD%PassedPrimaryInputData ) - if (IfWinputFilePassed) call Print_FileInfo_Struct( CU, InitInp%IW_InitInp%PassedFileData ) + if (ADinputFilePassed==1_c_int) call Print_FileInfo_Struct( CU, InitInp%AD%PassedPrimaryInputData ) + if (IfWinputFilePassed==1_c_int) call Print_FileInfo_Struct( CU, InitInp%IW_InitInp%PassedFileData ) endif ! Store data about the simulation (NOTE: we are not fully populating the Sim data structure) - allocate (Sim%WT(1),stat=errStat2); if (Failed0('wind turbines')) return Sim%dT = REAL(DT_C, DbKi) Sim%TMax = REAL(TMax_C, DbKi) Sim%numSteps = ceiling(Sim%tMax/Sim%dt) - Sim%NumTurbines = 1_IntKi ! only one turbine for now - Sim%WT(1)%NumBlades = int(NumBlades_C, IntKi) Sim%root = trim(OutRootName) - Sim%WT(1)%OriginInit = (/ 0.0_SiKi, 0.0_SiKi, 0.0_SiKi /) !TODO: should this be an input? ! Timekeeping n_Global = 0_IntKi ! Assume we are on timestep 0 at start n_VTK = -1_IntKi ! Set VTK output to T=0 at first call ! Interpolation order - InterpOrder = int(InterpOrder_C, IntKi) + InterpOrder = int(InterpOrder_C, IntKi) ! VTK outputs WrOutputsData%WrVTK = int(WrVTK_in, IntKi) @@ -386,9 +538,6 @@ SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinpu ! Validate and set some inputs (moved to subroutine to make cleaner to read call ValidateSetInputs(ErrStat2,ErrMsg2); if(Failed()) return - ! Flag to transpose DCMs as they are passed in - TransposeDCM = TransposeDCM_in - ! Linearization ! for now, set linearization to false. Pass this in later when interface supports it InitInp%AD%Linearize = .FALSE. @@ -406,56 +555,16 @@ SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinpu InitInp%AD%defPvap = REAL(defPvap_C, ReKi) InitInp%AD%WtrDpth = REAL(WtrDpth_C, ReKi) InitInp%AD%MSL2SWL = REAL(MSL2SWL_C, ReKi) - InitInp%storeHHVel = storeHHVel + InitInp%storeHHVel = storeHHVel==1_c_int InitInp%WrVTK = WrOutputsData%WrVTK InitInp%WrVTK_Type = WrOutputsData%WrVTK_Type InitInp%IW_InitInp%CompInflow = 1 ! Use InflowWind ! setup rotors for AD -- interface only supports one rotor at present - allocate (InitInp%AD%rotors(1),stat=errStat2); if (Failed0('rotors')) return - InitInp%AD%rotors(1)%AeroProjMod = int(AeroProjMod_C, IntKi) - InitInp%AD%rotors(1)%numBlades = Sim%WT(1)%NumBlades - call AllocAry(InitInp%AD%rotors(1)%BladeRootPosition, 3, Sim%WT(1)%NumBlades, 'BldRootPos', errStat2, errMsg2 ); if (Failed()) return - call AllocAry(InitInp%AD%rotors(1)%BladeRootOrientation, 3, 3, Sim%WT(1)%NumBlades, 'BldRootOri', errStat2, errMsg2 ); if (Failed()) return - InitInp%AD%rotors(1)%HubPosition = real(HubPos_C(1:3),ReKi) - InitInp%AD%rotors(1)%HubOrientation = reshape( real(HubOri_C(1:9),R8Ki), (/3,3/) ) - InitInp%AD%rotors(1)%NacellePosition = real(NacPos_C(1:3),ReKi) - InitInp%AD%rotors(1)%NacelleOrientation = reshape( real(NacOri_C(1:9),R8Ki), (/3,3/) ) - InitInp%AD%rotors(1)%BladeRootPosition = reshape( real(BldRootPos_C(1:3*Sim%WT(1)%NumBlades),ReKi), (/ 3,Sim%WT(1)%NumBlades/) ) - InitInp%AD%rotors(1)%BladeRootOrientation = reshape( real(BldRootOri_C(1:9*Sim%WT(1)%NumBlades),R8Ki), (/3,3,Sim%WT(1)%NumBlades/) ) - if (TransposeDCM) then - InitInp%AD%rotors(1)%HubOrientation = transpose(InitInp%AD%rotors(1)%HubOrientation) - InitInp%AD%rotors(1)%NacelleOrientation = transpose(InitInp%AD%rotors(1)%NacelleOrientation) - do i=1,Sim%WT(1)%NumBlades - InitInp%AD%rotors(1)%BladeRootOrientation(1:3,1:3,i) = transpose(InitInp%AD%rotors(1)%BladeRootOrientation(1:3,1:3,i)) - enddo - endif - - ! Remap the orientation DCM just in case there is some issue with passed - call OrientRemap(InitInp%AD%rotors(1)%HubOrientation) - call OrientRemap(InitInp%AD%rotors(1)%NacelleOrientation) - do i=1,Sim%WT(1)%NumBlades - call OrientRemap(InitInp%AD%rotors(1)%BladeRootOrientation(1:3,1:3,i)) + do iWT=1,Sim%NumTurbines + InitInp%AD%rotors(iWT)%numBlades = Sim%WT(iWT)%NumBlades enddo - ! Number of blades and initial positions - ! - NumMeshPts is the number of interface Mesh points we are expecting on the python - ! side. Will validate this against what AD reads from the initialization info. - NumMeshPts = int(NumMeshPts_C, IntKi) - if (NumMeshPts < 1) then - ErrStat2 = ErrID_Fatal - ErrMsg2 = "At least one node point must be specified" - if (Failed()) return - endif - ! Allocate temporary arrays to simplify data conversions - call AllocAry( tmpBldPtMeshPos, 3, NumMeshPts, "tmpBldPtMeshPos", ErrStat2, ErrMsg2 ); if (Failed()) return - call AllocAry( tmpBldPtMeshOri, 3, 3, NumMeshPts, "tmpBldPtMeshOri", ErrStat2, ErrMsg2 ); if (Failed()) return - call AllocAry( tmpBldPtMeshVel, 6, NumMeshPts, "tmpBldPtMeshVel", ErrStat2, ErrMsg2 ); if (Failed()) return - call AllocAry( tmpBldPtMeshAcc, 6, NumMeshPts, "tmpBldPtMeshAcc", ErrStat2, ErrMsg2 ); if (Failed()) return - call AllocAry( tmpBldPtMeshFrc, 6, NumMeshPts, "tmpBldPtMeshFrc", ErrStat2, ErrMsg2 ); if (Failed()) return - tmpBldPtMeshPos( 1:3,1:NumMeshPts) = reshape( real(InitMeshPos_C(1:3*NumMeshPts),ReKi), (/ 3,NumMeshPts/) ) - tmpBldPtMeshOri(1:3,1:3,1:NumMeshPts) = reshape( real(InitMeshOri_C(1:9*NumMeshPts),ReKi), (/3,3,NumMeshPts/) ) - !---------------------------------------------------- ! Allocate input array u and corresponding InputTimes @@ -485,17 +594,25 @@ SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinpu !------------------------------------------------------------- ! Sanity checks !------------------------------------------------------------- - call CheckNodes(ErrStat2,ErrMsg2); if (Failed()) return + do iWT=1,Sim%NumTurbines + call CheckNodes(iWT); if (Failed()) return + enddo !------------------------------------------------------------- ! Set the interface meshes for motion inputs and loads output !------------------------------------------------------------- - call SetMotionLoadsInterfaceMeshes(ErrStat2,ErrMsg2); if (Failed()) return + call SetupMotionLoadsInterfaceMeshes(); if (Failed()) return + ! setup meshes if (WrOutputsData%WrVTK > 0_IntKi) then - call setVTKParameters(WrOutputsData, Sim, ADI, errStat, errMsg, 'vtk-ADI') + call setVTKParameters(WrOutputsData, Sim, ADI, ErrStat2, ErrMsg2, 'vtk-ADI') if (Failed()) return - call WrVTK_refMeshes(ADI%u(1)%AD%rotors(:),WrOutputsData%VTKRefPoint,ErrStat2,ErrMsg2) + endif + ! write meshes for this rotor + if (WrOutputsData%WrVTK > 0_IntKi) then + do iWT=1,Sim%NumTurbines + call WrVTK_refMeshes(ADI%u(1)%AD%rotors(:),WrOutputsData%VTKRefPoint,ErrStat2,ErrMsg2) + enddo if (Failed()) return endif @@ -515,6 +632,13 @@ SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinpu enddo InputTimePrev = ADI%InputTimes(1) - Sim%dT ! Initialize for UpdateStates + !------------------------------------------------------------- + ! copy of ADI inputs. AD_SetInputMotion will set this mesh. When CalcOutput is called, + ! this data is used. When UpdateStates is called, this data is copied over to the ADI%u + !------------------------------------------------------------- + call ADI_CopyInput (ADI%u(1), ADI_u, MESH_NEWCOPY, Errstat2, ErrMsg2) + if (Failed()) return + !------------------------------------------------------------- ! Initial setup of other pieces of x,xd,z,OtherState @@ -562,6 +686,11 @@ SUBROUTINE AeroDyn_Inflow_C_Init( ADinputFilePassed, ADinputFileString_C, ADinpu call SetupFileOutputs() endif + + ! destroy the InitInp and InitOutput + call ADI_DestroyInitInput( InitInp, Errstat2, ErrMsg2); if (Failed()) return + call ADI_DestroyInitOutput(InitOutData, Errstat2, ErrMsg2); if (Failed()) return + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) @@ -570,7 +699,7 @@ logical function Failed() CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev if (Failed) then - call FailCleanup() + call ClearTmpStorage() call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) endif end function Failed @@ -578,31 +707,30 @@ end function Failed ! check for failed where /= 0 is fatal logical function Failed0(txt) character(*), intent(in) :: txt - if (errStat /= 0) then + if (errStat2 /= 0) then ErrStat2 = ErrID_Fatal ErrMsg2 = "Could not allocate "//trim(txt) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) endif Failed0 = ErrStat >= AbortErrLev - if(Failed0) call FailCleanup() + if(Failed0) then + call ClearTmpStorage() + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + endif end function Failed0 - subroutine FailCleanup() - if (allocated(tmpBldPtMeshPos)) deallocate(tmpBldPtMeshPos) - if (allocated(tmpBldPtMeshOri)) deallocate(tmpBldPtMeshOri) - if (allocated(tmpBldPtMeshVel)) deallocate(tmpBldPtMeshVel) - if (allocated(tmpBldPtMeshAcc)) deallocate(tmpBldPtMeshAcc) - if (allocated(tmpBldPtMeshFrc)) deallocate(tmpBldPtMeshFrc) - end subroutine FailCleanup - !> Validate and set some of the outputs (values must be stored before here as some might be changed) subroutine ValidateSetInputs(ErrStat3,ErrMsg3) integer(IntKi), intent( out) :: ErrStat3 !< temporary error status character(ErrMsgLen), intent( out) :: ErrMsg3 !< temporary error message + + ErrStat3 = ErrID_None + ErrMsg3 = "" + ! Interporder if ( InterpOrder < 1_IntKi .or. InterpOrder > 2_IntKi ) then - call SetErrStat(ErrID_Fatal,"InterpOrder passed into AeroDyn_Inflow_C_Init must be 1 (linear) or 2 (quadratic)",ErrStat3,ErrMsg3,RoutineName) + call SetErrStat(ErrID_Fatal,"InterpOrder passed into ADI_C_Init must be 1 (linear) or 2 (quadratic)",ErrStat3,ErrMsg3,RoutineName) return endif @@ -622,6 +750,7 @@ subroutine ValidateSetInputs(ErrStat3,ErrMsg3) endif endif + ! check fileFmt if ( WrOutputsData%fileFmt /= idFmtNone .and. WrOutputsData%fileFmt /= idFmtAscii .and. & WrOutputsData%fileFmt /= idFmtBinary .and. WrOutputsData%fileFmt /= idFmtBoth) then @@ -673,14 +802,17 @@ end subroutine SetupFileOutputs subroutine ShowPassedData() character(1) :: TmpFlag integer :: i,j - call WrScr("Interface debugging: Variables passed in through interface") + call WrSCr("") call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: Variables passed in through interface") + call WrScr(" ADI_C_Init") + call WrScr(" --------------------------------------------------------") call WrScr(" FileInfo") - TmpFlag="F"; if (ADinputFilePassed) TmpFlag="T" + TmpFlag="F"; if (ADinputFilePassed==1_c_int) TmpFlag="T" call WrScr(" ADinputFilePassed_C "//TmpFlag ) call WrScr(" ADinputFileString_C (ptr addr) "//trim(Num2LStr(LOC(ADinputFileString_C))) ) call WrScr(" ADinputFileStringLength_C "//trim(Num2LStr( ADinputFileStringLength_C )) ) - TmpFlag="F"; if (IfWinputFilePassed) TmpFlag="T" + TmpFlag="F"; if (IfWinputFilePassed==1_c_int) TmpFlag="T" call WrScr(" IfWinputFilePassed_C "//TmpFlag ) call WrScr(" IfWinputFileString_C (ptr addr)"//trim(Num2LStr(LOC(IfWinputFileString_C))) ) call WrScr(" IfWinputFileStringLength_C "//trim(Num2LStr( IfWinputFileStringLength_C )) ) @@ -703,208 +835,96 @@ subroutine ShowPassedData() call WrScr(" wrOuts_C "//trim(Num2LStr( wrOuts_C )) ) call WrScr(" DT_Outs_C "//trim(Num2LStr( DT_Outs_C )) ) call WrScr(" Flags") - TmpFlag="F"; if (storeHHVel) TmpFlag="T" + TmpFlag="F"; if (storeHHVel==1_c_int) TmpFlag="T" call WrScr(" storeHHVel "//TmpFlag ) call WrScr(" WrVTK_in "//trim(Num2LStr( WrVTK_in )) ) call WrScr(" WrVTK_inType "//trim(Num2LStr( WrVTK_inType )) ) - TmpFlag="F"; if (TransposeDCM_in) TmpFlag="T" - call WrScr(" TransposeDCM_in "//TmpFlag ) - call WrScr(" Init Data") - call WrNR(" Hub Position ") - call WrMatrix(HubPos_C,CU,'(3(ES15.7e2))') - call WrNR(" Hub Orientation ") - call WrMatrix(HubOri_C,CU,'(9(ES23.15e2))') - call WrNR(" Nacelle Position ") - call WrMatrix(NacPos_C,CU,'(3(ES15.7e2))') - call WrNR(" Nacelle Orientation ") - call WrMatrix(NacOri_C,CU,'(9(ES23.15e2))') - call WrScr(" NumBlades_C "//trim(Num2LStr( NumBlades_C )) ) - if (debugverbose > 1) then - call WrScr(" Root Positions") - do i=1,NumBlades_C - j=3*(i-1) - call WrMatrix(BldRootPos_C(j+1:j+3),CU,'(3(ES15.7e2))') - enddo - call WrScr(" Root Orientations") - do i=1,NumBlades_C - j=9*(i-1) - call WrMatrix(BldRootOri_C(j+1:j+9),CU,'(9(ES23.15e2))') - enddo - endif - call WrScr(" NumMeshPts_C "//trim(Num2LStr( NumMeshPts_C )) ) - if (debugverbose > 1) then - call WrScr(" Mesh Positions") - do i=1,NumMeshPts_C - j=3*(i-1) - call WrMatrix(InitMeshPos_C(j+1:j+3),CU,'(3(ES15.7e2))') - enddo - call WrScr(" Mesh Orientations") - do i=1,NumMeshPts_C - j=9*(i-1) - call WrMatrix(InitMeshOri_C(j+1:j+9),CU,'(9(ES23.15e2))') - enddo - endif call WrScr("-----------------------------------------------------------") end subroutine ShowPassedData !> This subroutine sets the interface meshes to map to the input motions to the AD !! meshes - subroutine SetMotionLoadsInterfaceMeshes(ErrStat3,ErrMsg3) - integer(IntKi), intent( out) :: ErrStat3 !< temporary error status - character(ErrMsgLen), intent( out) :: ErrMsg3 !< temporary error message - integer(IntKi) :: iNode - real(ReKi) :: InitPos(3) - real(R8Ki) :: Orient(3,3) - !------------------------------------------------------------- - ! Set the interface meshes for motion inputs and loads output - !------------------------------------------------------------- - ! Motion mesh for blades - call MeshCreate( BldPtMotionMesh , & - IOS = COMPONENT_INPUT , & - Nnodes = NumMeshPts , & - ErrStat = ErrStat3 , & - ErrMess = ErrMsg3 , & - TranslationDisp = .TRUE., Orientation = .TRUE., & - TranslationVel = .TRUE., RotationVel = .TRUE., & - TranslationAcc = .TRUE., RotationAcc = .FALSE. ) - if (ErrStat3 >= AbortErrLev) return - - do iNode=1,NumMeshPts - ! initial position and orientation of node - InitPos = tmpBldPtMeshPos(1:3,iNode) - if (TransposeDCM) then - Orient = transpose(tmpBldPtMeshOri(1:3,1:3,iNode)) - else - Orient = tmpBldPtMeshOri(1:3,1:3,iNode) - endif - call OrientRemap(Orient) - call MeshPositionNode( BldPtMotionMesh , & - iNode , & - InitPos , & ! position - ErrStat3, ErrMsg3 , & - Orient ) ! orientation - if (ErrStat3 >= AbortErrLev) return -!FIXME: if we need to switch to line2 instead of point, do that here. - call MeshConstructElement ( BldPtMotionMesh, ELEMENT_POINT, ErrStat3, ErrMsg3, iNode ) - if (ErrStat3 >= AbortErrLev) return + subroutine SetupMotionLoadsInterfaceMeshes() + integer(IntKi) :: iWT !< current rotor/turbine + integer(IntKi) :: iNode + integer(IntKi) :: maxBlades !< find out maximum number of blades on all turbine rotors + + maxBlades = 0_IntKi + do iWT=1,Sim%NumTurbines + maxBlades = max(maxBlades,Sim%WT(iWT)%NumBlades) enddo - call MeshCommit ( BldPtMotionMesh, ErrStat3, ErrMsg3 ) - if (ErrStat3 >= AbortErrLev) return - BldPtMotionMesh%RemapFlag = .TRUE. - - ! For checking the mesh, uncomment this. - ! note: CU is is output unit (platform dependent). - if (debugverbose >= 4) call MeshPrintInfo( CU, BldPtMotionMesh, MeshName='BldPtMotionMesh' ) - - -! !------------------------------------------------------------- -! ! Motion mesh for nacelle -- TODO: add this mesh for nacelle load transfers -! call MeshCreate( NacMotionMesh , & -! IOS = COMPONENT_INPUT , & -! Nnodes = 1 , & -! ErrStat = ErrStat3 , & -! ErrMess = ErrMsg3 , & -! TranslationDisp = .TRUE., Orientation = .TRUE., & -! TranslationVel = .TRUE., RotationVel = .TRUE., & -! TranslationAcc = .TRUE., RotationAcc = .FALSE. ) -! if (ErrStat3 >= AbortErrLev) return -! -! InitPos = real(NacPos_C( 1:3),ReKi) -! Orient = reshape( real(NacOri_C(1:9),ReKi), (/3,3/) ) -! call OrientRemap(Orient) -! call MeshPositionNode( NacMotionMesh , & -! 1 , & -! InitPos , & ! position -! ErrStat3, ErrMsg3 , & -! Orient ) ! orientation -! if (ErrStat3 >= AbortErrLev) return -! -! call MeshConstructElement ( NacMotionMesh, ELEMENT_POINT, ErrStat3, ErrMsg3, p1=1 ) -! if (ErrStat3 >= AbortErrLev) return -! -! call MeshCommit ( NacMotionMesh, ErrStat3, ErrMsg3 ) -! if (ErrStat3 >= AbortErrLev) return -! NacMotionMesh%RemapFlag = .TRUE. -! -! ! For checking the mesh, uncomment this. -! ! note: CU is is output unit (platform dependent). -! if (debugverbose >= 4) call MeshPrintInfo( CU, NacMotionMesh, MeshName='NacMotionMesh' ) -! + allocate(Map_BldPtMotion_2_AD_Blade( maxBlades,Sim%NumTurbines),STAT=ErrStat2); if (Failed0('Map_BldPtMotion_2_AD_Blade' )) return + allocate(Map_AD_BldLoad_P_2_BldPtLoad(maxBlades,Sim%NumTurbines),STAT=ErrStat2); if (Failed0('Map_AD_BldLoad_P_2_BldPtLoad')) return + + ! step through all turbine rotors + do iWT=1,Sim%NumTurbines + + !------------------------------------------------------------- + ! Load mesh for blades + CALL MeshCopy( SrcMesh = BldPtMotionMesh(iWT) ,& + DestMesh = BldPtLoadMesh(iWT) ,& + CtrlCode = MESH_SIBLING ,& + IOS = COMPONENT_OUTPUT ,& + ErrStat = ErrStat2 ,& + ErrMess = ErrMsg2 ,& + Force = .TRUE. ,& + Moment = .TRUE. ) + if(Failed()) return + BldPtLoadMesh(iWT)%RemapFlag = .FALSE. + + ! Temp mesh for load transfer + CALL MeshCopy( SrcMesh = BldPtLoadMesh(iWT) ,& + DestMesh = BldPtLoadMesh_tmp(iWT),& + CtrlCode = MESH_COUSIN ,& + IOS = COMPONENT_OUTPUT ,& + ErrStat = ErrStat2 ,& + ErrMess = ErrMsg2 ,& + Force = .TRUE. ,& + Moment = .TRUE. ) + if(Failed()) return + BldPtLoadMesh_tmp(iWT)%RemapFlag = .FALSE. + + + ! For checking the mesh + ! note: CU is is output unit (platform dependent). + if (debugverbose >= 4) call MeshPrintInfo( CU, BldPtLoadMesh(iWT), MeshName='BldPtLoadMesh'//trim(Num2LStr(iWT)) ) + + +! !------------------------------------------------------------- +! ! Load mesh for nacelle +! CALL MeshCopy( SrcMesh = NacMotionMesh(iWT) ,& +! DestMesh = NacLoadMesh(iWT) ,& +! CtrlCode = MESH_SIBLING ,& +! IOS = COMPONENT_OUTPUT ,& +! ErrStat = ErrStat2 ,& +! ErrMess = ErrMsg2 ,& +! Force = .TRUE. ,& +! Moment = .TRUE. ) +! if(Failed()) return +! NacLoadMesh(iWT)%RemapFlag = .FALSE. ! - !------------------------------------------------------------- - ! Load mesh for blades - CALL MeshCopy( SrcMesh = BldPtMotionMesh ,& - DestMesh = BldPtLoadMesh ,& - CtrlCode = MESH_SIBLING ,& - IOS = COMPONENT_OUTPUT ,& - ErrStat = ErrStat3 ,& - ErrMess = ErrMsg3 ,& - Force = .TRUE. ,& - Moment = .TRUE. ) - if (ErrStat3 >= AbortErrLev) return - BldPtLoadMesh%RemapFlag = .TRUE. - - ! Temp mesh for load transfer - CALL MeshCopy( SrcMesh = BldPtLoadMesh ,& - DestMesh = BldPtLoadMesh_tmp ,& - CtrlCode = MESH_COUSIN ,& - IOS = COMPONENT_OUTPUT ,& - ErrStat = ErrStat3 ,& - ErrMess = ErrMsg3 ,& - Force = .TRUE. ,& - Moment = .TRUE. ) - if (ErrStat3 >= AbortErrLev) return - BldPtLoadMesh_tmp%RemapFlag = .TRUE. +! ! For checking the mesh, uncomment this. +! ! note: CU is is output unit (platform dependent). +! if (debugverbose >= 4) call MeshPrintInfo( CU, NacLoadMesh(iWT), MeshName='NacLoadMesh'//trim(Num2LStr(iWT)) ) - ! For checking the mesh - ! note: CU is is output unit (platform dependent). - if (debugverbose >= 4) call MeshPrintInfo( CU, BldPtLoadMesh, MeshName='BldPtLoadMesh' ) - - -! !------------------------------------------------------------- -! ! Load mesh for nacelle -- TODO: add this mesh for nacelle load transfers -! CALL MeshCopy( SrcMesh = NacMotionMesh ,& -! DestMesh = NacLoadMesh ,& -! CtrlCode = MESH_SIBLING ,& -! IOS = COMPONENT_OUTPUT ,& -! ErrStat = ErrStat3 ,& -! ErrMess = ErrMsg3 ,& -! Force = .TRUE. ,& -! Moment = .TRUE. ) -! if (ErrStat3 >= AbortErrLev) return -! NacLoadMesh%RemapFlag = .TRUE. -! -! ! For checking the mesh, uncomment this. -! ! note: CU is is output unit (platform dependent). -! if (debugverbose >= 4) call MeshPrintInfo( CU, NacLoadMesh, MeshName='NacLoadMesh' ) - + !------------------------------------------------------------- + ! Set the mapping meshes + ! blades + do i=1,Sim%WT(iWT)%NumBlades + call MeshMapCreate( BldPtMotionMesh(iWT), ADI%u(1)%AD%rotors(iWT)%BladeMotion(i), Map_BldPtMotion_2_AD_Blade(i,iWT), ErrStat2, ErrMsg2 ); if(Failed()) return + call MeshMapCreate( ADI%y%AD%rotors(iWT)%BladeLoad(i), BldPtLoadMesh(iWT), Map_AD_BldLoad_P_2_BldPtLoad(i,iWT), ErrStat2, ErrMsg2 ); if(Failed()) return + enddo + ! nacelle -- TODO: add this mesh for nacelle load transfers +! if ( ADI%y%AD%rotors(iWT)%NacelleLoad%Committed ) then +! call MeshMapCreate( NacMotionMesh(iWT), ADI%u(1)%AD%rotors(iWT)%NacelleMotion, Map_NacPtMotion_2_AD_Nac(iWT), ErrStat2, ErrMsg2 ); if(Failed()) return +! call MeshMapCreate( ADI%y%AD%rotors(iWT)%NacelleLoad, NacLoadMesh(iWT), Map_AD_Nac_2_NacPtLoad(iWT), ErrStat2, ErrMsg2 ); if(Failed()) return +! endif - !------------------------------------------------------------- - ! Set the mapping meshes - ! blades - allocate(Map_BldPtMotion_2_AD_Blade(Sim%WT(1)%NumBlades),Map_AD_BldLoad_P_2_BldPtLoad(Sim%WT(1)%NumBlades),STAT=ErrStat3) - if (ErrStat3 /= 0) then - ErrStat3 = ErrID_Fatal - ErrMsg3 = "Could not allocate Map_BldPtMotion_2_AD_Blade" - return - endif - do i=1,Sim%WT(1)%NumBlades - call MeshMapCreate( BldPtMotionMesh, ADI%u(1)%AD%rotors(1)%BladeMotion(i), Map_BldPtMotion_2_AD_Blade(i), ErrStat3, ErrMsg3 ) - if (ErrStat3 >= AbortErrLev) return - call MeshMapCreate( ADI%y%AD%rotors(1)%BladeLoad(i), BldPtLoadMesh, Map_AD_BldLoad_P_2_BldPtLoad(i), ErrStat3, ErrMsg3 ) - if (ErrStat3 >= AbortErrLev) return - enddo - ! nacelle -- TODO: add this mesh for nacelle load transfers -! if ( y%AD%rotors(1)%NacelleLoad%Committed ) then -! call MeshMapCreate( NacMotionMesh, ADI%u(1)%AD%rotors(1)%NacelleMotion, Map_NacPtMotion_2_AD_Nac, ErrStat3, ErrMsg3 ) -! if (ErrStat3 >= AbortErrLev) return -! call MeshMapCreate( ADI%y%AD%rotors(1)%NacelleLoad, NacLoadMesh, Map_AD_Nac_2_NacPtLoad, ErrStat3, ErrMsg3 ) -! if (ErrStat3 >= AbortErrLev) return -! endif + enddo ! iWT - end subroutine SetMotionLoadsInterfaceMeshes + end subroutine SetupMotionLoadsInterfaceMeshes !------------------------------------------------------------- @@ -912,29 +932,26 @@ end subroutine SetMotionLoadsInterfaceMeshes !! If more than one input node was passed in, but only a single AD node !! exists, then give error that too many !! nodes passed. - subroutine CheckNodes(ErrStat3,ErrMsg3) - integer(IntKi), intent( out) :: ErrStat3 !< temporary error status - character(ErrMsgLen), intent( out) :: ErrMsg3 !< temporary error message - ErrStat3 = ErrID_None - ErrMsg3 = "" + subroutine CheckNodes(iWT) + integer(IntKi), intent(in ) :: iWT !< current rotor/turbine ! FIXME: this is a placeholder in case we think of some sanity checks to perform. ! - some check that nodes make some sense -- might be caught in meshmapping ! - some checks on hub/nacelle being near middle of the rotor? Not sure if that matters end subroutine CheckNodes -END SUBROUTINE AeroDyn_Inflow_C_Init +END SUBROUTINE ADI_C_Init !!=============================================================================================================== !!--------------------------------------------- AeroDyn ReInit--------------------------------------------------- !!=============================================================================================================== !!TODO: finish this routine so it is usable if we need re-init capability for coupling -!SUBROUTINE AeroDyn_Inflow_C_ReInit( DT_C, TMax_C, & -! ErrStat_C, ErrMsg_C) BIND (C, NAME='AeroDyn_Inflow_C_ReInit') +!SUBROUTINE ADI_C_ReInit( DT_C, TMax_C, & +! ErrStat_C, ErrMsg_C) BIND (C, NAME='ADI_C_ReInit') ! implicit none !#ifndef IMPLICIT_DLLEXPORT -!!DEC$ ATTRIBUTES DLLEXPORT :: AeroDyn_Inflow_C_ReInit -!!GCC$ ATTRIBUTES DLLEXPORT :: AeroDyn_Inflow_C_ReInit +!!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_ReInit +!!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_ReInit !#endif ! ! real(c_double), intent(in ) :: DT_C !< Timestep used with AD for stepping forward from t to t+dt. Must be constant. @@ -946,7 +963,7 @@ END SUBROUTINE AeroDyn_Inflow_C_Init ! character(ErrMsgLen) :: ErrMsg !< aggregated error message ! integer(IntKi) :: ErrStat2 !< temporary error status from a call ! character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call -! character(*), parameter :: RoutineName = 'AeroDyn_Inflow_C_ReInit' !< for error handling +! character(*), parameter :: RoutineName = 'ADI_C_ReInit' !< for error handling ! ! ! Initialize error handling ! ErrStat = ErrID_None @@ -969,103 +986,49 @@ END SUBROUTINE AeroDyn_Inflow_C_Init ! call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) ! endif ! end function Failed -!END SUBROUTINE AeroDyn_Inflow_C_ReInit +!END SUBROUTINE ADI_C_ReInit !=============================================================================================================== !--------------------------------------------- AeroDyn CalcOutput --------------------------------------------- !=============================================================================================================== - -SUBROUTINE AeroDyn_Inflow_C_CalcOutput(Time_C, & - HubPos_C, HubOri_C, HubVel_C, HubAcc_C, & - NacPos_C, NacOri_C, NacVel_C, NacAcc_C, & - BldRootPos_C, BldRootOri_C, BldRootVel_C, BldRootAcc_C, & - NumMeshPts_C, & - MeshPos_C, MeshOri_C, MeshVel_C, MeshAcc_C, & - MeshFrc_C, OutputChannelValues_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='AeroDyn_Inflow_C_CalcOutput') +!> This routine calculates the outputs at Time_C using the states and inputs provided. +!! NOTE: make sure to call ADI_C_SetRotorMotion before calling CalcOutput +SUBROUTINE ADI_C_CalcOutput(Time_C, & + OutputChannelValues_C, ErrStat_C, ErrMsg_C) BIND (C, NAME='ADI_C_CalcOutput') implicit none #ifndef IMPLICIT_DLLEXPORT -!DEC$ ATTRIBUTES DLLEXPORT :: AeroDyn_Inflow_C_CalcOutput -!GCC$ ATTRIBUTES DLLEXPORT :: AeroDyn_Inflow_C_CalcOutput +!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_CalcOutput +!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_CalcOutput #endif real(c_double), intent(in ) :: Time_C - real(c_float), intent(in ) :: HubPos_C( 3 ) !< Hub position - real(c_double), intent(in ) :: HubOri_C( 9 ) !< Hub orientation - real(c_float), intent(in ) :: HubVel_C( 6 ) !< Hub velocity - real(c_float), intent(in ) :: HubAcc_C( 6 ) !< Hub acceleration - real(c_float), intent(in ) :: NacPos_C( 3 ) !< Nacelle position - real(c_double), intent(in ) :: NacOri_C( 9 ) !< Nacelle orientation - real(c_float), intent(in ) :: NacVel_C( 6 ) !< Nacelle velocity - real(c_float), intent(in ) :: NacAcc_C( 6 ) !< Nacelle acceleration - real(c_float), intent(in ) :: BldRootPos_C( 3*Sim%WT(1)%NumBlades ) !< Blade root positions - real(c_double), intent(in ) :: BldRootOri_C( 9*Sim%WT(1)%NumBlades ) !< Blade root orientations - real(c_float), intent(in ) :: BldRootVel_C( 6*Sim%WT(1)%NumBlades ) !< Blade root velocities - real(c_float), intent(in ) :: BldRootAcc_C( 6*Sim%WT(1)%NumBlades ) !< Blade root accelerations - ! Blade mesh nodes - integer(c_int), intent(in ) :: NumMeshPts_C !< Number of mesh points we are transfering motions to and output loads to - real(c_float), intent(in ) :: MeshPos_C( 3*NumMeshPts_C ) !< A 3xNumMeshPts_C array [x,y,z] - real(c_double), intent(in ) :: MeshOri_C( 9*NumMeshPts_C ) !< A 9xNumMeshPts_C array [r11,r12,r13,r21,r22,r23,r31,r32,r33] - real(c_float), intent(in ) :: MeshVel_C( 6*NumMeshPts_C ) !< A 6xNumMeshPts_C array [x,y,z] - real(c_float), intent(in ) :: MeshAcc_C( 6*NumMeshPts_C ) !< A 6xNumMeshPts_C array [x,y,z] - real(c_float), intent( out) :: MeshFrc_C( 6*NumMeshPts_C ) !< A 6xNumMeshPts_C array [Fx,Fy,Fz,Mx,My,Mz] -- forces and moments (global) real(c_float), intent( out) :: OutputChannelValues_C(ADI%p%NumOuts) integer(c_int), intent( out) :: ErrStat_C character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) ! Local variables real(DbKi) :: Time - integer(IntKi) :: iNode integer(IntKi) :: ErrStat !< aggregated error status character(ErrMsgLen) :: ErrMsg !< aggregated error message integer(IntKi) :: ErrStat2 !< temporary error status from a call character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call - character(*), parameter :: RoutineName = 'AeroDyn_Inflow_C_CalcOutput' !< for error handling + character(*), parameter :: RoutineName = 'ADI_C_CalcOutput' !< for error handling ! Initialize error handling ErrStat = ErrID_None ErrMsg = "" - ! Sanity check -- number of node points cannot change - if ( NumMeshPts /= int(NumMeshPts_C, IntKi) ) then - ErrStat2 = ErrID_Fatal - ErrMsg2 = "Number of node points passed in changed. This must be constant throughout simulation" - if (Failed()) return - endif - - ! Convert the inputs from C to Fortrn Time = REAL(Time_C,DbKi) - ! Reshape mesh position, orientation, velocity, acceleration - tmpBldPtMeshPos(1:3,1:NumMeshPts) = reshape( real(MeshPos_C(1:3*NumMeshPts),ReKi), (/3, NumMeshPts/) ) - tmpBldPtMeshOri(1:3,1:3,1:NumMeshPts) = reshape( real(MeshOri_C(1:9*NumMeshPts),R8Ki), (/3,3,NumMeshPts/) ) - tmpBldPtMeshVel(1:6,1:NumMeshPts) = reshape( real(MeshVel_C(1:6*NumMeshPts),ReKi), (/6, NumMeshPts/) ) - tmpBldPtMeshAcc(1:6,1:NumMeshPts) = reshape( real(MeshAcc_C(1:6*NumMeshPts),ReKi), (/6, NumMeshPts/) ) - - - ! Transfer motions to input meshes - call Set_MotionMesh( ErrStat2, ErrMsg2 ); if (Failed()) return - call AD_SetInputMotion( ADI%u(1), & - HubPos_C, HubOri_C, HubVel_C, HubAcc_C, & - NacPos_C, NacOri_C, NacVel_C, NacAcc_C, & - BldRootPos_C, BldRootOri_C, BldRootVel_C, BldRootAcc_C, & - ErrStat2, ErrMsg2 ) ! transfer input motion mesh to u(1) meshes - if (Failed()) return - ! Call the main subroutine ADI_CalcOutput to get the resulting forces and moments at time T + call ADI_CopyInput (ADI_u, ADI%u(1), MESH_UPDATECOPY, Errstat2, ErrMsg2) ! copy new inputs over + if (Failed()) return CALL ADI_CalcOutput( Time, ADI%u(1), ADI%p, ADI%x(STATE_CURR), ADI%xd(STATE_CURR), ADI%z(STATE_CURR), ADI%OtherState(STATE_CURR), ADI%y, ADI%m, ErrStat2, ErrMsg2 ) if (Failed()) return - ! Transfer resulting load meshes to intermediate mesh - call AD_TransferLoads( ADI%u(1), ADI%y, ErrStat2, ErrMsg2 ) - if (Failed()) return - - ! Set output force/moment array - call Set_OutputLoadArray( ) - MeshFrc_C(1:6*NumMeshPts) = reshape( real(tmpBldPtMeshFrc(1:6,1:NumMeshPts), c_float), (/6*NumMeshPts/) ) - - ! Get the output channel info out of y - OutputChannelValues_C = REAL(ADI%y%WriteOutput, C_FLOAT) + ! Get the output channel info out of y + OutputChannelValues_C = REAL(ADI%y%WriteOutput, C_FLOAT) !------------------------------------------------------- ! write outputs @@ -1087,7 +1050,7 @@ logical function Failed() Failed = ErrStat >= AbortErrLev if (Failed) call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) end function Failed -END SUBROUTINE AeroDyn_Inflow_C_CalcOutput +END SUBROUTINE ADI_C_CalcOutput !=============================================================================================================== !--------------------------------------------- AeroDyn UpdateStates ------------------------------------------- @@ -1097,62 +1060,32 @@ END SUBROUTINE AeroDyn_Inflow_C_CalcOutput !! Since we don't really know if we are doing correction steps or not, we will track the previous state and !! reset to those if we are repeating a timestep (normally this would be handled by the OF glue code, but since !! the states are not passed across the interface, we must handle them here). -SUBROUTINE AeroDyn_Inflow_C_UpdateStates( Time_C, TimeNext_C, & - HubPos_C, HubOri_C, HubVel_C, HubAcc_C, & - NacPos_C, NacOri_C, NacVel_C, NacAcc_C, & - BldRootPos_C, BldRootOri_C, BldRootVel_C, BldRootAcc_C, & - NumMeshPts_C, & - MeshPos_C, MeshOri_C, MeshVel_C, MeshAcc_C, & - ErrStat_C, ErrMsg_C) BIND (C, NAME='AeroDyn_Inflow_C_UpdateStates') +!! NOTE: make sure to call ADI_C_SetRotorMotion before calling UpdateStates +SUBROUTINE ADI_C_UpdateStates( Time_C, TimeNext_C, & + ErrStat_C, ErrMsg_C) BIND (C, NAME='ADI_C_UpdateStates') implicit none #ifndef IMPLICIT_DLLEXPORT -!DEC$ ATTRIBUTES DLLEXPORT :: AeroDyn_Inflow_C_UpdateStates -!GCC$ ATTRIBUTES DLLEXPORT :: AeroDyn_Inflow_C_UpdateStates +!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_UpdateStates +!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_UpdateStates #endif real(c_double), intent(in ) :: Time_C real(c_double), intent(in ) :: TimeNext_C - real(c_float), intent(in ) :: HubPos_C( 3 ) !< Hub position - real(c_double), intent(in ) :: HubOri_C( 9 ) !< Hub orientation - real(c_float), intent(in ) :: HubVel_C( 6 ) !< Hub velocity - real(c_float), intent(in ) :: HubAcc_C( 6 ) !< Hub acceleration - real(c_float), intent(in ) :: NacPos_C( 3 ) !< Nacelle position - real(c_double), intent(in ) :: NacOri_C( 9 ) !< Nacelle orientation - real(c_float), intent(in ) :: NacVel_C( 6 ) !< Nacelle velocity - real(c_float), intent(in ) :: NacAcc_C( 6 ) !< Nacelle acceleration - real(c_float), intent(in ) :: BldRootPos_C( 3*Sim%WT(1)%NumBlades ) !< Blade root positions - real(c_double), intent(in ) :: BldRootOri_C( 9*Sim%WT(1)%NumBlades ) !< Blade root orientations - real(c_float), intent(in ) :: BldRootVel_C( 6*Sim%WT(1)%NumBlades ) !< Blade root velocities - real(c_float), intent(in ) :: BldRootAcc_C( 6*Sim%WT(1)%NumBlades ) !< Blade root accelerations - ! Blade mesh nodes - integer(c_int), intent(in ) :: NumMeshPts_C !< Number of mesh points we are transfering motions to and output loads to - real(c_float), intent(in ) :: MeshPos_C( 3*NumMeshPts_C ) !< A 3xNumMeshPts_C array [x,y,z] - real(c_double), intent(in ) :: MeshOri_C( 9*NumMeshPts_C ) !< A 9xNumMeshPts_C array [r11,r12,r13,r21,r22,r23,r31,r32,r33] - real(c_float), intent(in ) :: MeshVel_C( 6*NumMeshPts_C ) !< A 6xNumMeshPts_C array [x,y,z] - real(c_float), intent(in ) :: MeshAcc_C( 6*NumMeshPts_C ) !< A 6xNumMeshPts_C array [x,y,z] integer(c_int), intent( out) :: ErrStat_C character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) ! Local variables logical :: CorrectionStep ! if we are repeating a timestep in UpdateStates, don't update the inputs array - integer(IntKi) :: iNode integer(IntKi) :: ErrStat !< aggregated error status character(ErrMsgLen) :: ErrMsg !< aggregated error message integer(IntKi) :: ErrStat2 !< temporary error status from a call character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call - character(*), parameter :: RoutineName = 'AeroDyn_Inflow_C_UpdateStates' !< for error handling + character(*), parameter :: RoutineName = 'ADI_C_UpdateStates' !< for error handling ! Initialize error handling ErrStat = ErrID_None ErrMsg = "" CorrectionStep = .false. - ! Sanity check -- number of node points cannot change - if ( NumMeshPts /= int(NumMeshPts_C, IntKi) ) then - ErrStat2 = ErrID_Fatal - ErrMsg2 = "Number of node points passed in changed. This must be constant throughout simulation" - if (Failed()) return - endif - !------------------------------------------------------- ! Check the time for current timestep and next timestep @@ -1205,24 +1138,6 @@ SUBROUTINE AeroDyn_Inflow_C_UpdateStates( Time_C, TimeNext_C, & call ADI_CopyInput( ADI%u(INPUT_PRED), ADI%u(INPUT_CURR), MESH_UPDATECOPY, ErrStat2, ErrMsg2); if (Failed()) return endif - !------------------------------------------------------- - ! Set inputs for time T+dt -- u(INPUT_PRED) - !------------------------------------------------------- - ! Reshape mesh position, orientation, velocity, acceleration - tmpBldPtMeshPos(1:3,1:NumMeshPts) = reshape( real(MeshPos_C(1:3*NumMeshPts),ReKi), (/3, NumMeshPts/) ) - tmpBldPtMeshOri(1:3,1:3,1:NumMeshPts) = reshape( real(MeshOri_C(1:9*NumMeshPts),R8Ki), (/3,3,NumMeshPts/) ) - tmpBldPtMeshVel(1:6,1:NumMeshPts) = reshape( real(MeshVel_C(1:6*NumMeshPts),ReKi), (/6, NumMeshPts/) ) - tmpBldPtMeshAcc(1:6,1:NumMeshPts) = reshape( real(MeshAcc_C(1:6*NumMeshPts),ReKi), (/6, NumMeshPts/) ) - - ! Transfer motions to input meshes - call Set_MotionMesh( ErrStat2, ErrMsg2 ); if (Failed()) return - call AD_SetInputMotion( ADI%u(INPUT_PRED), & - HubPos_C, HubOri_C, HubVel_C, HubAcc_C, & - NacPos_C, NacOri_C, NacVel_C, NacAcc_C, & - BldRootPos_C, BldRootOri_C, BldRootVel_C, BldRootAcc_C, & - ErrStat2, ErrMsg2 ) ! transfer input motion mesh to u(1) meshes - if (Failed()) return - ! Set copy the current state over to the predicted state for sending to UpdateStates ! -- The STATE_PREDicted will get updated in the call. @@ -1233,6 +1148,11 @@ SUBROUTINE AeroDyn_Inflow_C_UpdateStates( Time_C, TimeNext_C, & CALL ADI_CopyOtherState (ADI%OtherState(STATE_CURR), ADI%OtherState(STATE_PRED), MESH_UPDATECOPY, Errstat2, ErrMsg2); if (Failed()) return + ! Copy newinputs for time u(INPUT_PRED) + call ADI_CopyInput (ADI_u, ADI%u(INPUT_PRED), MESH_UPDATECOPY, Errstat2, ErrMsg2) + if (Failed()) return + + ! Call the main subroutine ADI_UpdateStates to get the velocities CALL ADI_UpdateStates( ADI%InputTimes(INPUT_CURR), n_Global, ADI%u, ADI%InputTimes, ADI%p, ADI%x(STATE_PRED), ADI%xd(STATE_PRED), ADI%z(STATE_PRED), ADI%OtherState(STATE_PRED), ADI%m, ErrStat2, ErrMsg2 ) if (Failed()) return @@ -1264,18 +1184,18 @@ logical function Failed() Failed = ErrStat >= AbortErrLev if (Failed) call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) end function Failed -END SUBROUTINE AeroDyn_Inflow_C_UpdateStates +END SUBROUTINE ADI_C_UpdateStates !=============================================================================================================== !--------------------------------------------------- AeroDyn End----------------------------------------------- !=============================================================================================================== ! NOTE: the error handling in this routine is slightly different than the other routines -SUBROUTINE AeroDyn_Inflow_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='AeroDyn_Inflow_C_End') +SUBROUTINE ADI_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='ADI_C_End') implicit none #ifndef IMPLICIT_DLLEXPORT -!DEC$ ATTRIBUTES DLLEXPORT :: AeroDyn_Inflow_C_End -!GCC$ ATTRIBUTES DLLEXPORT :: AeroDyn_Inflow_C_End +!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_End +!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_End #endif integer(c_int), intent( out) :: ErrStat_C character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) @@ -1288,7 +1208,7 @@ SUBROUTINE AeroDyn_Inflow_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='AeroDyn_Inflo character(ErrMsgLen) :: ErrMsg !< aggregated error message integer :: ErrStat2 !< temporary error status from a call character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call - character(*), parameter :: RoutineName = 'AeroDyn_Inflow_C_End' !< for error handling + character(*), parameter :: RoutineName = 'ADI_C_End' !< for error handling ! Initialize error handling ErrStat = ErrID_None @@ -1309,19 +1229,12 @@ SUBROUTINE AeroDyn_Inflow_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='AeroDyn_Inflo else sWT = '' endif - call WrBinFAST(trim(WrOutputsData%Root)//trim(sWT)//'.outb', FileFmtID_ChanLen_In, 'AeroDyn_Inflow_C_Library', WrOutputsData%WriteOutputHdr, WrOutputsData%WriteOutputUnt, (/0.0_DbKi, Sim%dT/), WrOutputsData%storage(:,:,iWT), errStat2, errMsg2) + call WrBinFAST(trim(WrOutputsData%Root)//trim(sWT)//'.outb', FileFmtID_ChanLen_In, 'ADI_C_Library', WrOutputsData%WriteOutputHdr, WrOutputsData%WriteOutputUnt, (/0.0_DbKi, Sim%dT/), WrOutputsData%storage(:,:,iWT), errStat2, errMsg2) call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) enddo endif end if - ! clear out any globably allocated helper arrays - if (allocated(tmpBldPtMeshPos)) deallocate(tmpBldPtMeshPos) - if (allocated(tmpBldPtMeshOri)) deallocate(tmpBldPtMeshOri) - if (allocated(tmpBldPtMeshVel)) deallocate(tmpBldPtMeshVel) - if (allocated(tmpBldPtMeshAcc)) deallocate(tmpBldPtMeshAcc) - if (allocated(tmpBldPtMeshFrc)) deallocate(tmpBldPtMeshFrc) - ! Call the main subroutine ADI_End ! If u is not allocated, then we didn't get far at all in initialization, @@ -1362,75 +1275,580 @@ SUBROUTINE AeroDyn_Inflow_C_End(ErrStat_C,ErrMsg_C) BIND (C, NAME='AeroDyn_Inflo !if (allocated(InputTimes)) deallocate(InputTimes) ! Clear out mesh related data storage - call ClearMesh() + call ClearTmpStorage() + + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) +END SUBROUTINE ADI_C_End + + +!=============================================================================================================== +!--------------------------------------------- AeroDyn SetupRotor ---------------------------------------------- +!=============================================================================================================== +!> Setup the initial rotor root positions etc before initializing +subroutine ADI_C_SetupRotor(iWT_c, TurbineIsHAWT_c, TurbOrigin_C, & + HubPos_C, HubOri_C, & + NacPos_C, NacOri_C, & + NumBlades_C, BldRootPos_C, BldRootOri_C, & + NumMeshPts_C, InitMeshPos_C, InitMeshOri_C, & + ErrStat_C, ErrMsg_C) BIND (C, NAME='ADI_C_SetupRotor') + implicit none +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_SetupRotor +!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_SetupRotor +#endif + integer(c_int), intent(in ) :: iWT_c !< Wind turbine / rotor number + integer(c_int), intent(in ) :: TurbineIsHAWT_c !< true for HAWT, false for VAWT + real(c_float), intent(in ) :: TurbOrigin_C(3) !< turbine origin (tower base). Gets added to all meshes to shift turbine position. + ! Initial hub and blade root positions/orientations + real(c_float), intent(in ) :: HubPos_C( 3 ) !< Hub position + real(c_double), intent(in ) :: HubOri_C( 9 ) !< Hub orientation + real(c_float), intent(in ) :: NacPos_C( 3 ) !< Nacelle position + real(c_double), intent(in ) :: NacOri_C( 9 ) !< Nacelle orientation + integer(c_int), intent(in ) :: NumBlades_C !< Number of blades + real(c_float), intent(in ) :: BldRootPos_C( 3*NumBlades_C ) !< Blade root positions + real(c_double), intent(in ) :: BldRootOri_C( 9*NumBlades_C ) !< Blade root orientations + ! Initial nodes + integer(c_int), intent(in ) :: NumMeshPts_C !< Number of mesh points we are transfering motions to and output loads to + real(c_float), intent(in ) :: InitMeshPos_C( 3*NumMeshPts_C ) !< A 3xNumMeshPts_C array [x,y,z] + real(c_double), intent(in ) :: InitMeshOri_C( 9*NumMeshPts_C ) !< A 9xNumMeshPts_C array [r11,r12,r13,r21,r22,r23,r31,r32,r33] + integer(c_int), intent( out) :: ErrStat_C !< Error status + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) !< Error message (C_NULL_CHAR terminated) + + ! local vars + integer(IntKi) :: iWT !< current turbine + logical :: TurbineIsHAWT !< true for HAWT, false for VAWT + integer(IntKi) :: ErrStat !< aggregated error messagNumBlades_ee + character(ErrMsgLen) :: ErrMsg !< aggregated error message + integer(IntKi) :: ErrStat2 !< temporary error status from a call + character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call + integer(IntKi) :: i,j,k !< generic counters + character(*), parameter :: RoutineName = 'ADI_C_Init' !< for error handling + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + + ! For debugging the interface: + if (debugverbose > 0) then + call ShowPassedData() + endif + + ! turbine geometry + iWT = int(iWT_c, IntKi) + Sim%WT(iWT)%NumBlades = int(NumBlades_C, IntKi) + Sim%WT(iWT)%OriginInit(1:3) = real(TurbOrigin_C(1:3), ReKi) + + ! Aero calculation method -- AeroProjMod + ! APM_BEM_NoSweepPitchTwist - 1 - "Original AeroDyn model where momentum balance is done in the WithoutSweepPitchTwist system" + ! APM_BEM_Polar - 2 - "Use staggered polar grid for momentum balance in each annulus" + ! APM_LiftingLine - 3 - "Use the blade lifting line (i.e. the structural) orientation (currently for OLAF with VAWT)" + ! For now we will set (this may need to be changed later): + ! HAWT --> AeroProjMod = 1 + ! VAWT --> AeroProjMod = 3 + TurbineIsHAWT = TurbineIsHAWT_c==1_c_int + if (TurbineIsHAWT) then + InitInp%AD%rotors(iWT)%AeroProjMod = 1 + else + InitInp%AD%rotors(iWT)%AeroProjMod = 3 + endif + + + call AllocAry(InitInp%AD%rotors(iWT)%BladeRootPosition, 3, Sim%WT(iWT)%NumBlades, 'BldRootPos', errStat2, errMsg2 ); if (Failed()) return + call AllocAry(InitInp%AD%rotors(iWT)%BladeRootOrientation, 3, 3, Sim%WT(iWT)%NumBlades, 'BldRootOri', errStat2, errMsg2 ); if (Failed()) return + InitInp%AD%rotors(iWT)%originInit = Sim%WT(iWT)%OriginInit(1:3) + InitInp%AD%rotors(iWT)%HubPosition = real(HubPos_C(1:3),ReKi) + Sim%WT(iWT)%OriginInit(1:3) + InitInp%AD%rotors(iWT)%HubOrientation = reshape( real(HubOri_C(1:9),R8Ki), (/3,3/) ) + InitInp%AD%rotors(iWT)%NacellePosition = real(NacPos_C(1:3),ReKi) + Sim%WT(iWT)%OriginInit(1:3) + InitInp%AD%rotors(iWT)%NacelleOrientation = reshape( real(NacOri_C(1:9),R8Ki), (/3,3/) ) + InitInp%AD%rotors(iWT)%BladeRootPosition = reshape( real(BldRootPos_C(1:3*Sim%WT(iWT)%NumBlades),ReKi), (/ 3,Sim%WT(iWT)%NumBlades/) ) + InitInp%AD%rotors(iWT)%BladeRootOrientation = reshape( real(BldRootOri_C(1:9*Sim%WT(iWT)%NumBlades),R8Ki), (/3,3,Sim%WT(iWT)%NumBlades/) ) + do i=1,Sim%WT(iWT)%NumBlades + InitInp%AD%rotors(iWT)%BladeRootPosition(1:3,i) = InitInp%AD%rotors(iWT)%BladeRootPosition(1:3,i) + Sim%WT(iWT)%OriginInit(1:3) + enddo + if (TransposeDCM) then + InitInp%AD%rotors(iWT)%HubOrientation = transpose(InitInp%AD%rotors(iWT)%HubOrientation) + InitInp%AD%rotors(iWT)%NacelleOrientation = transpose(InitInp%AD%rotors(iWT)%NacelleOrientation) + do i=1,Sim%WT(iWT)%NumBlades + InitInp%AD%rotors(iWT)%BladeRootOrientation(1:3,1:3,i) = transpose(InitInp%AD%rotors(iWT)%BladeRootOrientation(1:3,1:3,i)) + enddo + endif + + ! Remap the orientation DCM just in case there is some issue with passed + call OrientRemap(InitInp%AD%rotors(iWT)%HubOrientation) + call OrientRemap(InitInp%AD%rotors(iWT)%NacelleOrientation) + do i=1,Sim%WT(iWT)%NumBlades + call OrientRemap(InitInp%AD%rotors(iWT)%BladeRootOrientation(1:3,1:3,i)) + enddo + + ! Number of blades and initial positions + ! - NumMeshPts is the number of interface Mesh points we are expecting on the python + ! side. Will validate this against what AD reads from the initialization info. + NumMeshPts(iWT) = int(NumMeshPts_C, IntKi) + if (NumMeshPts(iWT) < 1) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "At least one node point must be specified" + if (Failed()) return + endif + + call SetupMotionMesh() + + ! Set error status + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + +contains + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) then + call ClearTmpStorage() + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + endif + end function Failed + + ! check for failed where /= 0 is fatal + logical function Failed0(txt) + character(*), intent(in) :: txt + if (errStat /= 0) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Could not allocate "//trim(txt) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + endif + Failed0 = ErrStat >= AbortErrLev + if(Failed0) call ClearTmpStorage() + end function Failed0 + + !> This subroutine prints out all the variables that are passed in. Use this only + !! for debugging the interface on the Fortran side. + subroutine ShowPassedData() + character(1) :: TmpFlag + integer :: i,j + call WrSCr("") + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: Variables passed in through interface") + call WrScr(" ADI_C_SetupRotor -- rotor "//trim(Num2LStr(iWT_c))) + call WrScr(" --------------------------------------------------------") + call WrScr(" Turbine origin") + call WrMatrix(TurbOrigin_C,CU,'(3(ES15.7e2))') + call WrScr(" Init rotor positions/orientations (positions do not include Turbine origin offset)") + call WrNR(" Hub Position ") + call WrMatrix(HubPos_C,CU,'(3(ES15.7e2))') + call WrNR(" Hub Orientation ") + call WrMatrix(HubOri_C,CU,'(9(ES23.15e2))') + call WrNR(" Nacelle Position ") + call WrMatrix(NacPos_C,CU,'(3(ES15.7e2))') + call WrNR(" Nacelle Orientation ") + call WrMatrix(NacOri_C,CU,'(9(ES23.15e2))') + call WrScr(" NumBlades_C "//trim(Num2LStr(NumBlades_C)) ) + if (debugverbose > 1) then + call WrScr(" Root Positions") + do i=1,NumBlades_C + j=3*(i-1) + call WrMatrix(BldRootPos_C(j+1:j+3),CU,'(3(ES15.7e2))') + enddo + call WrScr(" Root Orientations") + do i=1,NumBlades_C + j=9*(i-1) + call WrMatrix(BldRootOri_C(j+1:j+9),CU,'(9(ES23.15e2))') + enddo + endif + call WrScr(" NumMeshPts_C "//trim(Num2LStr( NumMeshPts_C )) ) + if (debugverbose > 1) then + call WrScr(" Mesh Positions") + do i=1,NumMeshPts_C + j=3*(i-1) + call WrMatrix(InitMeshPos_C(j+1:j+3),CU,'(3(ES15.7e2))') + enddo + call WrScr(" Mesh Orientations") + do i=1,NumMeshPts_C + j=9*(i-1) + call WrMatrix(InitMeshOri_C(j+1:j+9),CU,'(9(ES23.15e2))') + enddo + endif + call WrScr("-----------------------------------------------------------") + end subroutine ShowPassedData + + subroutine SetupMotionMesh() + real(ReKi) :: InitPos(3) + real(R8Ki) :: Orient(3,3) + integer(IntKi) :: iNode + ! this may not be super efficient since we are allocating for every turbine, but this makes things a bit simpler + if ( allocated(tmpBldPtMeshPos) ) deallocate(tmpBldPtMeshPos) + if ( allocated(tmpBldPtMeshOri) ) deallocate(tmpBldPtMeshOri) + call AllocAry( tmpBldPtMeshPos, 3, NumMeshPts(iWT), "tmpBldPtMeshPos", ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( tmpBldPtMeshOri, 3, 3, NumMeshPts(iWT), "tmpBldPtMeshOri", ErrStat2, ErrMsg2 ); if (Failed()) return + + ! Reshape mesh position, orientation, velocity, acceleration + tmpBldPtMeshPos(1:3,1:NumMeshPts(iWT)) = reshape( real(InitMeshPos_C(1:3*NumMeshPts(iWT)),ReKi), (/3, NumMeshPts(iWT)/) ) + tmpBldPtMeshOri(1:3,1:3,1:NumMeshPts(iWT)) = reshape( real(InitMeshOri_C(1:9*NumMeshPts(iWT)),R8Ki), (/3,3,NumMeshPts(iWT)/) ) + + !------------------------------------------------------------- + ! Set the interface meshes for motion inputs and loads output + !------------------------------------------------------------- + ! Motion mesh for blades + call MeshCreate( BldPtMotionMesh(iWT) , & + IOS = COMPONENT_INPUT , & + Nnodes = NumMeshPts(iWT) , & + ErrStat = ErrStat2 , & + ErrMess = ErrMsg2 , & + TranslationDisp = .TRUE., Orientation = .TRUE., & + TranslationVel = .TRUE., RotationVel = .TRUE., & + TranslationAcc = .TRUE., RotationAcc = .FALSE. ) + if(Failed()) return + + do iNode=1,NumMeshPts(iWT) + ! initial position and orientation of node + InitPos = tmpBldPtMeshPos(1:3,iNode) + Sim%WT(iWT)%OriginInit(1:3) + if (TransposeDCM) then + Orient = transpose(tmpBldPtMeshOri(1:3,1:3,iNode)) + else + Orient = tmpBldPtMeshOri(1:3,1:3,iNode) + endif + call OrientRemap(Orient) + call MeshPositionNode( BldPtMotionMesh(iWT) , & + iNode , & + InitPos , & ! position + ErrStat2, ErrMsg2 , & + Orient ) ! orientation + if(Failed()) return + call MeshConstructElement ( BldPtMotionMesh(iWT), ELEMENT_POINT, ErrStat2, ErrMsg2, iNode ); if(Failed()) return + enddo + + call MeshCommit ( BldPtMotionMesh(iWT), ErrStat2, ErrMsg2 ); if(Failed()) return + BldPtMotionMesh(iWT)%RemapFlag = .FALSE. + + ! For checking the mesh, uncomment this. + ! note: CU is is output unit (platform dependent). + if (debugverbose >= 4) call MeshPrintInfo( CU, BldPtMotionMesh(iWT), MeshName='BldPtMotionMesh'//trim(Num2LStr(iWT)) ) + + +! !------------------------------------------------------------- +! ! Motion mesh for nacelle -- TODO: add this mesh for nacelle load transfers +! call MeshCreate( NacMotionMesh(iWT) , & +! IOS = COMPONENT_INPUT , & +! Nnodes = 1 , & +! ErrStat = ErrStat2 , & +! ErrMess = ErrMsg2 , & +! TranslationDisp = .TRUE., Orientation = .TRUE., & +! TranslationVel = .TRUE., RotationVel = .TRUE., & +! TranslationAcc = .TRUE., RotationAcc = .FALSE. ) +! if(Failed()) return +! +! InitPos = real(NacPos_C( 1:3),ReKi) + Sim%WT(iWT)%OriginInit(1:3) +! Orient = reshape( real(NacOri_C(1:9),ReKi), (/3,3/) ) +! call OrientRemap(Orient) +! call MeshPositionNode( NacMotionMesh(iWT) , & +! 1 , & +! InitPos , & ! position +! ErrStat2, ErrMsg2 , & +! Orient ) ! orientation +! if(Failed()) return +! +! call MeshConstructElement ( NacMotionMesh(iWT), ELEMENT_POINT, ErrStat2, ErrMsg2, p1=1 ); if(Failed()) return +! +! call MeshCommit ( NacMotionMesh(iWT), ErrStat2, ErrMsg2 ); if(Failed()) return +! NacMotionMesh(iWT)%RemapFlag = .FALSE. +! +! ! For checking the mesh, uncomment this. +! ! note: CU is is output unit (platform dependent). +! if (debugverbose >= 4) call MeshPrintInfo( CU, NacMotionMesh(iWT), MeshName='NacMotionMesh'//trim(Num2LStr(iWT)) ) + + ! clear the tmp stuff (will need to resize this later) + if ( allocated(tmpBldPtMeshPos) ) deallocate(tmpBldPtMeshPos) + if ( allocated(tmpBldPtMeshOri) ) deallocate(tmpBldPtMeshOri) + end subroutine SetupMotionMesh +end subroutine ADI_C_SetupRotor + +!=============================================================================================================== +!--------------------------------------------- AeroDyn SetRotorMotion ------------------------------------------ +!=============================================================================================================== +!> Set the motions for a single rotor. This must be called before ADI_C_CalcOutput +subroutine ADI_C_SetRotorMotion( iWT_c, & + HubPos_C, HubOri_C, HubVel_C, HubAcc_C, & + NacPos_C, NacOri_C, NacVel_C, NacAcc_C, & + BldRootPos_C, BldRootOri_C, BldRootVel_C, BldRootAcc_C, & + NumMeshPts_C, & + MeshPos_C, MeshOri_C, MeshVel_C, MeshAcc_C, & + ErrStat_C, ErrMsg_C) BIND (C, NAME='ADI_C_SetRotorMotion') + implicit none +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_SetRotorMotion +!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_SetRotorMotion +#endif + integer(c_int), intent(in ) :: iWT_c !< Wind turbine / rotor number + real(c_float), intent(in ) :: HubPos_C( 3 ) !< Hub position + real(c_double), intent(in ) :: HubOri_C( 9 ) !< Hub orientation + real(c_float), intent(in ) :: HubVel_C( 6 ) !< Hub velocity + real(c_float), intent(in ) :: HubAcc_C( 6 ) !< Hub acceleration + real(c_float), intent(in ) :: NacPos_C( 3 ) !< Nacelle position + real(c_double), intent(in ) :: NacOri_C( 9 ) !< Nacelle orientation + real(c_float), intent(in ) :: NacVel_C( 6 ) !< Nacelle velocity + real(c_float), intent(in ) :: NacAcc_C( 6 ) !< Nacelle acceleration + real(c_float), intent(in ) :: BldRootPos_C( 3*Sim%WT(iWT_c)%NumBlades ) !< Blade root positions + real(c_double), intent(in ) :: BldRootOri_C( 9*Sim%WT(iWT_c)%NumBlades ) !< Blade root orientations + real(c_float), intent(in ) :: BldRootVel_C( 6*Sim%WT(iWT_c)%NumBlades ) !< Blade root velocities + real(c_float), intent(in ) :: BldRootAcc_C( 6*Sim%WT(iWT_c)%NumBlades ) !< Blade root accelerations + ! Blade mesh nodes + integer(c_int), intent(in ) :: NumMeshPts_C !< Number of mesh points we are transfering motions to and output loads to + real(c_float), intent(in ) :: MeshPos_C( 3*NumMeshPts_C ) !< A 3xNumMeshPts_C array [x,y,z] + real(c_double), intent(in ) :: MeshOri_C( 9*NumMeshPts_C ) !< A 9xNumMeshPts_C array [r11,r12,r13,r21,r22,r23,r31,r32,r33] + real(c_float), intent(in ) :: MeshVel_C( 6*NumMeshPts_C ) !< A 6xNumMeshPts_C array [x,y,z] + real(c_float), intent(in ) :: MeshAcc_C( 6*NumMeshPts_C ) !< A 6xNumMeshPts_C array [x,y,z] + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + + ! Local variables + real(DbKi) :: Time + integer(IntKi) :: iNode + integer(IntKi) :: iWT !< current wind turbine / rotor + integer(IntKi) :: ErrStat !< aggregated error status + character(ErrMsgLen) :: ErrMsg !< aggregated error message + integer(IntKi) :: ErrStat2 !< temporary error status from a call + character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call + character(*), parameter :: RoutineName = 'ADI_C_SetRotorMotion' !< for error handling + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + ! For debugging the interface: + if (debugverbose > 0) then + call ShowPassedData() + endif + + ! current turbine number + iWT = int(iWT_c, IntKi) + + ! Sanity check -- number of node points cannot change + if ( NumMeshPts(iWT) /= int(NumMeshPts_C, IntKi) ) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Number of node points passed in changed. This must be constant throughout simulation" + if (Failed()) return + endif + + ! Reshape mesh position, orientation, velocity, acceleration + tmpBldPtMeshPos(1:3,1:NumMeshPts(iWT)) = reshape( real(MeshPos_C(1:3*NumMeshPts(iWT)),ReKi), (/3, NumMeshPts(iWT)/) ) + tmpBldPtMeshOri(1:3,1:3,1:NumMeshPts(iWT)) = reshape( real(MeshOri_C(1:9*NumMeshPts(iWT)),R8Ki), (/3,3,NumMeshPts(iWT)/) ) + tmpBldPtMeshVel(1:6,1:NumMeshPts(iWT)) = reshape( real(MeshVel_C(1:6*NumMeshPts(iWT)),ReKi), (/6, NumMeshPts(iWT)/) ) + tmpBldPtMeshAcc(1:6,1:NumMeshPts(iWT)) = reshape( real(MeshAcc_C(1:6*NumMeshPts(iWT)),ReKi), (/6, NumMeshPts(iWT)/) ) + + ! Transfer motions to input meshes + do iWT=1,Sim%NumTurbines + call Set_MotionMesh(iWT, ErrStat2, ErrMsg2); if (Failed()) return + call AD_SetInputMotion( iWT, ADI_u, & + HubPos_C, HubOri_C, HubVel_C, HubAcc_C, & + NacPos_C, NacOri_C, NacVel_C, NacAcc_C, & + BldRootPos_C, BldRootOri_C, BldRootVel_C, BldRootAcc_C, & + ErrStat2, ErrMsg2 ) ! transfer input motion mesh to u(1) meshes + if (Failed()) return + enddo + + ! Set error status call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + CONTAINS - !> Don't leave junk in memory. So destroy meshes and mappings. - subroutine ClearMesh() - ! Blade - call MeshDestroy( BldPtMotionMesh, ErrStat2, ErrMsg2 ) - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - call MeshDestroy( BldPtLoadMesh, ErrStat2, ErrMsg2 ) - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - ! Destroy mesh mappings - if (allocated(Map_BldPtMotion_2_AD_Blade)) then - do i=1,Sim%WT(1)%NumBlades - call NWTC_Library_Destroymeshmaptype( Map_BldPtMotion_2_AD_Blade(i), ErrStat2, ErrMsg2 ) - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + end function Failed + !> This subroutine prints out all the variables that are passed in. Use this only + !! for debugging the interface on the Fortran side. + subroutine ShowPassedData() + character(1) :: TmpFlag + integer :: i,j + call WrSCr("") + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: Variables passed in through interface") + call WrScr(" ADI_C_SetRotorMotion -- rotor "//trim(Num2LStr(iWT_c))) + call WrScr(" ("//trim(Num2LStr(Sim%WT(iWT_C)%numBlades))//" blades, "//trim(Num2LStr(NumMeshPts(iWT_C)))//" mesh nodes)") + call WrScr(" --------------------------------------------------------") + call WrScr(" rotor positions/orientations (positions do not include Turbine origin offset)") + call WrNR(" Hub Position ") + call WrMatrix(HubPos_C,CU,'(3(ES15.7e2))') + call WrNR(" Hub Orientation ") + call WrMatrix(HubOri_C,CU,'(9(ES23.15e2))') + call WrNR(" Hub Velocity ") + call WrMatrix(HubVel_C,CU,'(6(ES15.7e2))') + call WrNR(" Hub Acceleration ") + call WrMatrix(HubAcc_C,CU,'(6(ES15.7e2))') + + call WrNR(" Nacelle Position ") + call WrMatrix(NacPos_C,CU,'(3(ES15.7e2))') + call WrNR(" Nacelle Orientation ") + call WrMatrix(NacOri_C,CU,'(9(ES23.15e2))') + call WrNR(" Nacelle Velocity ") + call WrMatrix(NacVel_C,CU,'(6(ES15.7e2))') + call WrNR(" Nacelle Acceleration ") + call WrMatrix(NacAcc_C,CU,'(6(ES15.7e2))') + + if (debugverbose > 1) then + call WrScr(" Root Positions (positions do not include Turbine origin offset)") + do i=1,Sim%WT(iWT_c)%NumBlades + j=3*(i-1) + call WrMatrix(BldRootPos_C(j+1:j+3),CU,'(3(ES15.7e2))') + enddo + call WrScr(" Root Orientations") + do i=1,Sim%WT(iWT_c)%NumBlades + j=9*(i-1) + call WrMatrix(BldRootOri_C(j+1:j+9),CU,'(9(ES23.15e2))') + enddo + call WrScr(" Root Velocities") + do i=1,Sim%WT(iWT_c)%NumBlades + j=3*(i-1) + call WrMatrix(BldRootPos_C(j+1:j+3),CU,'(6(ES15.7e2))') + enddo + call WrScr(" Root Accelerations") + do i=1,Sim%WT(iWT_c)%NumBlades + j=3*(i-1) + call WrMatrix(BldRootAcc_C(j+1:j+3),CU,'(6(ES15.7e2))') enddo - deallocate(Map_BldPtMotion_2_AD_Blade) endif - if (allocated(Map_AD_BldLoad_P_2_BldPtLoad)) then - do i=1,Sim%WT(1)%NumBlades - call NWTC_Library_Destroymeshmaptype( Map_AD_BldLoad_P_2_BldPtLoad(i), ErrStat2, ErrMsg2 ) - call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call WrScr(" NumMeshPts_C "//trim(Num2LStr( NumMeshPts_C )) ) + if (debugverbose > 1) then + call WrScr(" Mesh Positions (positions do not include Turbine origin offset)") + do i=1,NumMeshPts_C + j=3*(i-1) + call WrMatrix(MeshPos_C(j+1:j+3),CU,'(3(ES15.7e2))') + enddo + call WrScr(" Mesh Orientations") + do i=1,NumMeshPts_C + j=9*(i-1) + call WrMatrix(MeshOri_C(j+1:j+9),CU,'(9(ES23.15e2))') + enddo + call WrScr(" Mesh Velocities") + do i=1,NumMeshPts_C + j=3*(i-1) + call WrMatrix(MeshVel_C(j+1:j+3),CU,'(6(ES15.7e2))') + enddo + call WrScr(" Mesh Accelerations") + do i=1,NumMeshPts_C + j=3*(i-1) + call WrMatrix(MeshAcc_C(j+1:j+3),CU,'(6(ES15.7e2))') enddo - deallocate(Map_AD_BldLoad_P_2_BldPtLoad) endif - ! Nacelle -- TODO: add this mesh for nacelle load transfers -! call MeshDestroy( NacMotionMesh, ErrStat2, ErrMsg2 ) -! call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) -! call MeshDestroy( NacLoadMesh, ErrStat2, ErrMsg2 ) -! call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) -! call NWTC_Library_Destroymeshmaptype( Map_AD_Nac_2_NacPtLoad , ErrStat2, ErrMsg2 ) -! call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) -! call NWTC_Library_Destroymeshmaptype( Map_NacPtMotion_2_AD_Nac , ErrStat2, ErrMsg2 ) -! call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + call WrScr("-----------------------------------------------------------") + end subroutine ShowPassedData + +end subroutine ADI_C_SetRotorMotion + +!=============================================================================================================== +!--------------------------------------------- AeroDyn GetRotorLoads ------------------------------------------- +!=============================================================================================================== +!> Get the loads from a single rotor. This must be called after ADI_C_CalcOutput +subroutine ADI_C_GetRotorLoads(iWT_C, & + NumMeshPts_C, MeshFrc_C, & + ErrStat_C, ErrMsg_C) BIND (C, NAME='ADI_C_GetRotorLoads') + implicit none +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: ADI_C_GetRotorLoads +!GCC$ ATTRIBUTES DLLEXPORT :: ADI_C_GetRotorLoads +#endif + integer(c_int), intent(in ) :: iWT_C !< Wind turbine / rotor number + integer(c_int), intent(in ) :: NumMeshPts_C !< Number of mesh points we are transfering motions to and output loads to + real(c_float), intent( out) :: MeshFrc_C( 6*NumMeshPts_C ) !< A 6xNumMeshPts_C array [Fx,Fy,Fz,Mx,My,Mz] -- forces and moments (global) + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + + ! Local variables + integer(IntKi) :: iWT !< current wind turbine / rotor + integer(IntKi) :: ErrStat !< aggregated error status + character(ErrMsgLen) :: ErrMsg !< aggregated error message + integer(IntKi) :: ErrStat2 !< temporary error status from a call + character(ErrMsgLen) :: ErrMsg2 !< temporary error message from a call + character(*), parameter :: RoutineName = 'ADI_C_SetRotorMotion' !< for error handling + + ! Initialize error handling + ErrStat = ErrID_None + ErrMsg = "" + + ! For debugging the interface: + if (debugverbose > 0) then + call ShowPassedData() + endif + + ! current turbine number + iWT = int(iWT_c, IntKi) + + ! Sanity check -- number of node points cannot change + if ( NumMeshPts(iWT) /= int(NumMeshPts_C, IntKi) ) then + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Number of node points passed in changed. This must be constant throughout simulation" + if (Failed()) return + endif + + ! Transfer resulting load meshes to intermediate mesh + call AD_TransferLoads( iWT, ADI%u(1), ADI%y, ErrStat2, ErrMsg2 ) + if (Failed()) return + + ! Set output force/moment array + call Set_OutputLoadArray(iWT) + MeshFrc_C(1:6*NumMeshPts(iWT)) = reshape( real(tmpBldPtMeshFrc(1:6,1:NumMeshPts(iWT)), c_float), (/6*NumMeshPts(iWT)/) ) + + ! Set error status + call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + +CONTAINS + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + if (Failed) call SetErr(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + end function Failed + !> This subroutine prints out all the variables that are passed in. Use this only + !! for debugging the interface on the Fortran side. + subroutine ShowPassedData() + character(1) :: TmpFlag + integer :: i,j + call WrSCr("") + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: Variables passed in through interface") + call WrScr(" ADI_C_GetRotorLoads -- rotor "//trim(Num2LStr(iWT_c))) + call WrScr(" --------------------------------------------------------") + call WrScr(" NumMeshPts_C "//trim(Num2LStr( NumMeshPts_C )) ) + call WrScr("-----------------------------------------------------------") + end subroutine ShowPassedData +end subroutine ADI_C_GetRotorLoads - end subroutine ClearMesh -END SUBROUTINE AeroDyn_Inflow_C_End +!=================================================================================================================================== +! Internal routines for setting meshes etc. +!=================================================================================================================================== + !> This routine is operating on module level data. Error handling here in case checks added -subroutine Set_MotionMesh(ErrStat3, ErrMsg3) +!! NOTE: the OriginInit is not included in the data passed in and must be added to the the position info here +subroutine Set_MotionMesh(iWT, ErrStat3, ErrMsg3) + integer(IntKi), intent(in ) :: iWT !< current rotor/turbine integer(IntKi), intent( out) :: ErrStat3 character(ErrMsgLen), intent( out) :: ErrMsg3 integer(IntKi) :: iNode ErrStat3 = 0_IntKi ErrMsg3 = '' ! Set mesh corresponding to input motions - do iNode=1,NumMeshPts - BldPtMotionMesh%TranslationDisp(1:3,iNode) = tmpBldPtMeshPos(1:3,iNode) - real(BldPtMotionMesh%Position(1:3,iNode), R8Ki) - BldPtMotionMesh%Orientation(1:3,1:3,iNode) = tmpBldPtMeshOri(1:3,1:3,iNode) - BldPtMotionMesh%TranslationVel( 1:3,iNode) = tmpBldPtMeshVel(1:3,iNode) - BldPtMotionMesh%RotationVel( 1:3,iNode) = tmpBldPtMeshVel(4:6,iNode) - BldPtMotionMesh%TranslationAcc( 1:3,iNode) = tmpBldPtMeshAcc(1:3,iNode) - !BldPtMotionMesh%RotationAcc( 1:3,iNode) = tmpBldPtMeshAcc(4:6,iNode) ! Rotational acc not included - call OrientRemap(BldPtMotionMesh%Orientation(1:3,1:3,iNode)) + do iNode=1,NumMeshPts(iWT) + BldPtMotionMesh(iWT)%TranslationDisp(1:3,iNode) = tmpBldPtMeshPos(1:3,iNode) + Sim%WT(iWT)%OriginInit(1:3) - real(BldPtMotionMesh(iWT)%Position(1:3,iNode), R8Ki) + BldPtMotionMesh(iWT)%Orientation(1:3,1:3,iNode) = tmpBldPtMeshOri(1:3,1:3,iNode) + BldPtMotionMesh(iWT)%TranslationVel( 1:3,iNode) = tmpBldPtMeshVel(1:3,iNode) + BldPtMotionMesh(iWT)%RotationVel( 1:3,iNode) = tmpBldPtMeshVel(4:6,iNode) + BldPtMotionMesh(iWT)%TranslationAcc( 1:3,iNode) = tmpBldPtMeshAcc(1:3,iNode) + !BldPtMotionMesh(iWT)%RotationAcc( 1:3,iNode) = tmpBldPtMeshAcc(4:6,iNode) ! Rotational acc not included + call OrientRemap(BldPtMotionMesh(iWT)%Orientation(1:3,1:3,iNode)) if (TransposeDCM) then - BldPtMotionMesh%Orientation(1:3,1:3,iNode) = transpose(BldPtMotionMesh%Orientation(1:3,1:3,iNode)) + BldPtMotionMesh(iWT)%Orientation(1:3,1:3,iNode) = transpose(BldPtMotionMesh(iWT)%Orientation(1:3,1:3,iNode)) endif enddo end subroutine Set_MotionMesh !> Map the motion of the intermediate input mesh over to the input meshes !! This routine is operating on module level data, hence few inputs -subroutine AD_SetInputMotion( u_local, & +!! NOTE: the OriginInit is not included in the data passed in and must be added to the the position info here +subroutine AD_SetInputMotion( iWT, u_local, & HubPos_C, HubOri_C, HubVel_C, HubAcc_C, & NacPos_C, NacOri_C, NacVel_C, NacAcc_C, & BldRootPos_C, BldRootOri_C, BldRootVel_C, BldRootAcc_C, & ErrStat, ErrMsg ) + integer(IntKi), intent(in ) :: iWT !< this turbine type(ADI_InputType), intent(inout) :: u_local ! Only one input (probably at T) real(c_float), intent(in ) :: HubPos_C( 3 ) !< Hub position real(c_double), intent(in ) :: HubOri_C( 9 ) !< Hub orientation @@ -1440,58 +1858,58 @@ subroutine AD_SetInputMotion( u_local, & real(c_double), intent(in ) :: NacOri_C( 9 ) !< Nacelle orientation real(c_float), intent(in ) :: NacVel_C( 6 ) !< Nacelle velocity real(c_float), intent(in ) :: NacAcc_C( 6 ) !< Nacelle acceleration - real(c_float), intent(in ) :: BldRootPos_C( 3*Sim%WT(1)%NumBlades ) !< Blade root positions - real(c_double), intent(in ) :: BldRootOri_C( 9*Sim%WT(1)%NumBlades ) !< Blade root orientations - real(c_float), intent(in ) :: BldRootVel_C( 6*Sim%WT(1)%NumBlades ) !< Blade root velocities - real(c_float), intent(in ) :: BldRootAcc_C( 6*Sim%WT(1)%NumBlades ) !< Blade root accelerations + real(c_float), intent(in ) :: BldRootPos_C( 3*Sim%WT(iWT)%NumBlades ) !< Blade root positions + real(c_double), intent(in ) :: BldRootOri_C( 9*Sim%WT(iWT)%NumBlades ) !< Blade root orientations + real(c_float), intent(in ) :: BldRootVel_C( 6*Sim%WT(iWT)%NumBlades ) !< Blade root velocities + real(c_float), intent(in ) :: BldRootAcc_C( 6*Sim%WT(iWT)%NumBlades ) !< Blade root accelerations integer(IntKi), intent( out) :: ErrStat character(ErrMsgLen), intent( out) :: ErrMsg integer(IntKi) :: i ErrStat = 0_IntKi ErrMsg = '' ! Hub -- NOTE: RotationalAcc not present in the mesh - if ( u_local%AD%rotors(1)%HubMotion%Committed ) then - u_local%AD%rotors(1)%HubMotion%TranslationDisp(1:3,1) = real(HubPos_C(1:3),R8Ki) - real(u_local%AD%rotors(1)%HubMotion%Position(1:3,1), R8Ki) - u_local%AD%rotors(1)%HubMotion%Orientation(1:3,1:3,1) = reshape( real(HubOri_C(1:9),R8Ki), (/3,3/) ) - u_local%AD%rotors(1)%HubMotion%TranslationVel(1:3,1) = real(HubVel_C(1:3), ReKi) - u_local%AD%rotors(1)%HubMotion%RotationVel(1:3,1) = real(HubVel_C(4:6), ReKi) - u_local%AD%rotors(1)%HubMotion%TranslationAcc(1:3,1) = real(HubAcc_C(1:3), ReKi) - call OrientRemap(u_local%AD%rotors(1)%HubMotion%Orientation(1:3,1:3,1)) + if ( u_local%AD%rotors(iWT)%HubMotion%Committed ) then + u_local%AD%rotors(iWT)%HubMotion%TranslationDisp(1:3,1) = real(HubPos_C(1:3),R8Ki) + Sim%WT(iWT)%OriginInit(1:3) - real(u_local%AD%rotors(iWT)%HubMotion%Position(1:3,1), R8Ki) + u_local%AD%rotors(iWT)%HubMotion%Orientation(1:3,1:3,1) = reshape( real(HubOri_C(1:9),R8Ki), (/3,3/) ) + u_local%AD%rotors(iWT)%HubMotion%TranslationVel(1:3,1) = real(HubVel_C(1:3), ReKi) + u_local%AD%rotors(iWT)%HubMotion%RotationVel(1:3,1) = real(HubVel_C(4:6), ReKi) + u_local%AD%rotors(iWT)%HubMotion%TranslationAcc(1:3,1) = real(HubAcc_C(1:3), ReKi) + call OrientRemap(u_local%AD%rotors(iWT)%HubMotion%Orientation(1:3,1:3,1)) if (TransposeDCM) then - u_local%AD%rotors(1)%HubMotion%Orientation(1:3,1:3,1) = transpose(u_local%AD%rotors(1)%HubMotion%Orientation(1:3,1:3,1)) + u_local%AD%rotors(iWT)%HubMotion%Orientation(1:3,1:3,1) = transpose(u_local%AD%rotors(iWT)%HubMotion%Orientation(1:3,1:3,1)) endif endif ! Nacelle -- NOTE: RotationalVel and RotationalAcc not present in the mesh - if ( u_local%AD%rotors(1)%NacelleMotion%Committed ) then - u_local%AD%rotors(1)%NacelleMotion%TranslationDisp(1:3,1) = real(NacPos_C(1:3),R8Ki) - real(u_local%AD%rotors(1)%NacelleMotion%Position(1:3,1), R8Ki) - u_local%AD%rotors(1)%NacelleMotion%Orientation(1:3,1:3,1) = reshape( real(NacOri_C(1:9),R8Ki), (/3,3/) ) - u_local%AD%rotors(1)%NacelleMotion%TranslationVel(1:3,1) = real(NacVel_C(1:3), ReKi) - u_local%AD%rotors(1)%NacelleMotion%TranslationAcc(1:3,1) = real(NacAcc_C(1:3), ReKi) - call OrientRemap(u_local%AD%rotors(1)%NacelleMotion%Orientation(1:3,1:3,1)) + if ( u_local%AD%rotors(iWT)%NacelleMotion%Committed ) then + u_local%AD%rotors(iWT)%NacelleMotion%TranslationDisp(1:3,1) = real(NacPos_C(1:3),R8Ki) + Sim%WT(iWT)%OriginInit(1:3) - real(u_local%AD%rotors(iWT)%NacelleMotion%Position(1:3,1), R8Ki) + u_local%AD%rotors(iWT)%NacelleMotion%Orientation(1:3,1:3,1) = reshape( real(NacOri_C(1:9),R8Ki), (/3,3/) ) + u_local%AD%rotors(iWT)%NacelleMotion%TranslationVel(1:3,1) = real(NacVel_C(1:3), ReKi) + u_local%AD%rotors(iWT)%NacelleMotion%TranslationAcc(1:3,1) = real(NacAcc_C(1:3), ReKi) + call OrientRemap(u_local%AD%rotors(iWT)%NacelleMotion%Orientation(1:3,1:3,1)) if (TransposeDCM) then - u_local%AD%rotors(1)%NacelleMotion%Orientation(1:3,1:3,1) = transpose(u_local%AD%rotors(1)%NacelleMotion%Orientation(1:3,1:3,1)) + u_local%AD%rotors(iWT)%NacelleMotion%Orientation(1:3,1:3,1) = transpose(u_local%AD%rotors(iWT)%NacelleMotion%Orientation(1:3,1:3,1)) endif endif ! Blade root - do i=0,Sim%WT(1)%numBlades-1 - if ( u_local%AD%rotors(1)%BladeRootMotion(i+1)%Committed ) then - u_local%AD%rotors(1)%BladeRootMotion(i+1)%TranslationDisp(1:3,1) = real(BldRootPos_C(3*i+1:3*i+3),R8Ki) - real(u_local%AD%rotors(1)%BladeRootMotion(i+1)%Position(1:3,1), R8Ki) - u_local%AD%rotors(1)%BladeRootMotion(i+1)%Orientation(1:3,1:3,1) = reshape( real(BldRootOri_C(9*i+1:9*i+9),R8Ki), (/3,3/) ) - u_local%AD%rotors(1)%BladeRootMotion(i+1)%TranslationVel(1:3,1) = real(BldRootVel_C(6*i+1:6*i+3), ReKi) - u_local%AD%rotors(1)%BladeRootMotion(i+1)%RotationVel(1:3,1) = real(BldRootVel_C(6*i+4:6*i+6), ReKi) - u_local%AD%rotors(1)%BladeRootMotion(i+1)%TranslationAcc(1:3,1) = real(BldRootAcc_C(6*i+1:6*i+3), ReKi) - u_local%AD%rotors(1)%BladeRootMotion(i+1)%RotationAcc(1:3,1) = real(BldRootAcc_C(6*i+4:6*i+6), ReKi) - call OrientRemap(u_local%AD%rotors(1)%BladeRootMotion(i+1)%Orientation(1:3,1:3,1)) + do i=0,Sim%WT(iWT)%numBlades-1 + if ( u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%Committed ) then + u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%TranslationDisp(1:3,1) = real(BldRootPos_C(3*i+1:3*i+3),R8Ki) + Sim%WT(iWT)%OriginInit(1:3) - real(u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%Position(1:3,1), R8Ki) + u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%Orientation(1:3,1:3,1) = reshape( real(BldRootOri_C(9*i+1:9*i+9),R8Ki), (/3,3/) ) + u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%TranslationVel(1:3,1) = real(BldRootVel_C(6*i+1:6*i+3), ReKi) + u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%RotationVel(1:3,1) = real(BldRootVel_C(6*i+4:6*i+6), ReKi) + u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%TranslationAcc(1:3,1) = real(BldRootAcc_C(6*i+1:6*i+3), ReKi) + u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%RotationAcc(1:3,1) = real(BldRootAcc_C(6*i+4:6*i+6), ReKi) + call OrientRemap(u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%Orientation(1:3,1:3,1)) if (TransposeDCM) then - u_local%AD%rotors(1)%BladeRootMotion(i+1)%Orientation(1:3,1:3,1) = transpose(u_local%AD%rotors(1)%BladeRootMotion(i+1)%Orientation(1:3,1:3,1)) + u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%Orientation(1:3,1:3,1) = transpose(u_local%AD%rotors(iWT)%BladeRootMotion(i+1)%Orientation(1:3,1:3,1)) endif endif enddo ! Blade mesh - do i=1,Sim%WT(1)%numBlades - if ( u_local%AD%rotors(1)%BladeMotion(i)%Committed ) then - call Transfer_Point_to_Line2( BldPtMotionMesh, u_local%AD%rotors(1)%BladeMotion(i), Map_BldPtMotion_2_AD_Blade(i), ErrStat, ErrMsg ) + do i=1,Sim%WT(iWT)%numBlades + if ( u_local%AD%rotors(iWT)%BladeMotion(i)%Committed ) then + call Transfer_Point_to_Line2( BldPtMotionMesh(iWT), u_local%AD%rotors(iWT)%BladeMotion(i), Map_BldPtMotion_2_AD_Blade(i,iWT), ErrStat, ErrMsg ) if (ErrStat >= AbortErrLev) return endif enddo @@ -1499,35 +1917,37 @@ end subroutine AD_SetInputMotion !> Map the loads of the output mesh to the intermediate output mesh. !! This routine is operating on module level data, hence few inputs -subroutine AD_TransferLoads( u_local, y_local, ErrStat3, ErrMsg3 ) - type(ADI_InputType), intent(in ) :: u_local ! Only one input (probably at T) +subroutine AD_TransferLoads( iWT, u_local, y_local, ErrStat3, ErrMsg3 ) + integer(IntKi), intent(in ) :: iWT !< this turbine + type(ADI_InputType), intent(in ) :: u_local ! Only one input (probably at T) type(ADI_OutputType), intent(in ) :: y_local ! Only one input (probably at T) integer(IntKi), intent( out) :: ErrStat3 character(ErrMsgLen), intent( out) :: ErrMsg3 integer(IntKi) :: i - BldPtLoadMesh%Force = 0.0_ReKi - BldPtLoadMesh%Moment = 0.0_ReKi - do i=1,Sim%WT(1)%NumBlades - if ( y_local%AD%rotors(1)%BladeLoad(i)%Committed ) then - if (debugverbose > 4) call MeshPrintInfo( CU, y_local%AD%rotors(1)%BladeLoad(i), MeshName='AD%rotors('//trim(Num2LStr(1))//')%BladeLoad('//trim(Num2LStr(i))//')' ) - call Transfer_Line2_to_Point( ADI%y%AD%rotors(1)%BladeLoad(i), BldPtLoadMesh_tmp, Map_AD_BldLoad_P_2_BldPtLoad(i), & - ErrStat3, ErrMsg3, u_local%AD%rotors(1)%BladeMotion(i), BldPtMotionMesh ) + BldPtLoadMesh(iWT)%Force = 0.0_ReKi + BldPtLoadMesh(iWT)%Moment = 0.0_ReKi + do i=1,Sim%WT(iWT)%NumBlades + if ( y_local%AD%rotors(iWT)%BladeLoad(i)%Committed ) then + if (debugverbose > 4) call MeshPrintInfo( CU, y_local%AD%rotors(iWT)%BladeLoad(i), MeshName='AD%rotors('//trim(Num2LStr(1))//')%BladeLoad('//trim(Num2LStr(i))//')' ) + call Transfer_Line2_to_Point( ADI%y%AD%rotors(iWT)%BladeLoad(i), BldPtLoadMesh_tmp(iWT), Map_AD_BldLoad_P_2_BldPtLoad(i,iWT), & + ErrStat3, ErrMsg3, u_local%AD%rotors(iWT)%BladeMotion(i), BldPtMotionMesh(iWT) ) if (ErrStat3 >= AbortErrLev) return - BldPtLoadMesh%Force = BldPtLoadMesh%Force + BldPtLoadMesh_tmp%Force - BldPtLoadMesh%Moment = BldPtLoadMesh%Moment + BldPtLoadMesh_tmp%Moment + BldPtLoadMesh(iWT)%Force = BldPtLoadMesh(iWT)%Force + BldPtLoadMesh_tmp(iWT)%Force + BldPtLoadMesh(iWT)%Moment = BldPtLoadMesh(iWT)%Moment + BldPtLoadMesh_tmp(iWT)%Moment endif enddo - if (debugverbose > 4) call MeshPrintInfo( CU, BldPtLoadMesh, MeshName='BldPtLoadMesh' ) + if (debugverbose > 4) call MeshPrintInfo( CU, BldPtLoadMesh(iWT), MeshName='BldPtLoadMesh'//trim(Num2LStr(iWT)) ) end subroutine AD_TransferLoads !> Transfer the loads from the load mesh to the temporary array for output !! This routine is operating on module level data, hence few inputs -subroutine Set_OutputLoadArray() +subroutine Set_OutputLoadArray(iWT) + integer(IntKi), intent(in ) :: iWT !< current rotor/turbine integer(IntKi) :: iNode ! Set mesh corresponding to input motions - do iNode=1,NumMeshPts - tmpBldPtMeshFrc(1:3,iNode) = BldPtLoadMesh%Force (1:3,iNode) - tmpBldPtMeshFrc(4:6,iNode) = BldPtLoadMesh%Moment(1:3,iNode) + do iNode=1,NumMeshPts(iWT) + tmpBldPtMeshFrc(1:3,iNode) = BldPtLoadMesh(iWT)%Force (1:3,iNode) + tmpBldPtMeshFrc(4:6,iNode) = BldPtLoadMesh(iWT)%Moment(1:3,iNode) enddo end subroutine Set_OutputLoadArray @@ -1598,12 +2018,12 @@ subroutine WrVTK_PointsRef(ErrStat3,ErrMsg3) ErrMsg3 = '' ! Blade point motion (structural mesh from driver) - call MeshWrVTKreference(RefPoint, BldPtMotionMesh, trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BldPtMotionMesh', ErrStat3, ErrMsg3) + call MeshWrVTKreference(RefPoint, BldPtMotionMesh(iWT), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BldPtMotionMesh', ErrStat3, ErrMsg3) if (ErrStat3 >= AbortErrLev) return ! Blade root motion (point only) if (allocated(rot_u(iWT)%BladeRootMotion)) then - do k=1,Sim%WT(1)%NumBlades + do k=1,Sim%WT(iWT)%NumBlades if (rot_u(iWT)%BladeRootMotion(k)%Committed) then call MeshWrVTKreference(RefPoint, rot_u(iWT)%BladeRootMotion(k), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BladeRootMotion'//trim(num2lstr(k)), ErrStat3, ErrMsg3 ) if (ErrStat3 >= AbortErrLev) return @@ -1637,7 +2057,7 @@ subroutine WrVTK_LinesRef(ErrStat3,ErrMsg3) ! Blades if (allocated(rot_u(iWT)%BladeMotion)) then - do k=1,Sim%WT(1)%NumBlades + do k=1,Sim%WT(iWT)%NumBlades if (rot_u(iWT)%BladeMotion(k)%Committed) then call MeshWrVTKreference(RefPoint, rot_u(iWT)%BladeMotion(k), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.Blade'//trim(num2lstr(k)), ErrStat3, ErrMsg3 ) if (ErrStat3 >= AbortErrLev) return @@ -1703,12 +2123,12 @@ subroutine WrVTK_Points(ErrStat3,ErrMsg3) ErrMsg3 = '' ! Blade point motion (structural mesh from driver) - call MeshWrVTK(RefPoint, BldPtMotionMesh, trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BldPtMotionMesh', n_Global, .true., ErrStat3, ErrMsg3, WrOutputsData%VTK_tWidth) + call MeshWrVTK(RefPoint, BldPtMotionMesh(iWT), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BldPtMotionMesh', n_Global, .true., ErrStat3, ErrMsg3, WrOutputsData%VTK_tWidth) if (ErrStat3 >= AbortErrLev) return ! Blade root motion (point only) if (allocated(rot_u(iWT)%BladeRootMotion)) then - do k=1,Sim%WT(1)%NumBlades + do k=1,Sim%WT(iWT)%NumBlades if (rot_u(iWT)%BladeRootMotion(k)%Committed) then call MeshWrVTK(RefPoint, rot_u(iWT)%BladeRootMotion(k), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.BladeRootMotion'//trim(num2lstr(k)), n_Global, .true., ErrStat3, ErrMsg3, WrOutputsData%VTK_tWidth) if (ErrStat3 >= AbortErrLev) return @@ -1765,7 +2185,7 @@ subroutine WrVTK_Surfaces(ErrStat3,ErrMsg3) ! Blades if (allocated(rot_u(iWT)%BladeMotion)) then - do k=1,Sim%WT(1)%NumBlades + do k=1,Sim%WT(iWT)%NumBlades if (rot_u(iWT)%BladeMotion(k)%Committed) then call MeshWrVTK_Ln2Surface (RefPoint, rot_u(iWT)%BladeMotion(k), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.Blade'//trim(num2lstr(k))//'Surface', & n_Global, OutputFields, errStat3, errMsg3, WrOutputsData%VTK_tWidth , verts=ADI%m%VTK_Surfaces(iWT)%BladeShape(k)%AirfoilCoords, & @@ -1797,7 +2217,7 @@ subroutine WrVTK_Lines(ErrStat3,ErrMsg3) ! Blades if (allocated(rot_u(iWT)%BladeMotion)) then - do k=1,Sim%WT(1)%NumBlades + do k=1,Sim%WT(iWT)%NumBlades if (rot_u(iWT)%BladeMotion(k)%Committed) then call MeshWrVTK(RefPoint, rot_u(iWT)%BladeMotion(k), trim(WrOutputsData%VTK_OutFileRoot)//trim(sWT)//'.Blade'//trim(num2lstr(k)), n_Global, .true., ErrStat3, ErrMsg3, WrOutputsData%VTK_tWidth) if (ErrStat3 >= AbortErrLev) return @@ -1852,5 +2272,92 @@ subroutine WrVTK_Ground (RefPoint, HalfLengths, FileRootName, errStat, errMsg) end subroutine WrVTK_Ground +!-------------------------------------------------------------------- +!> Set some temporary data storage arrays to simplify data conversion +subroutine SetTempStorage(ErrStat,ErrMsg) + INTEGER(IntKi), intent(out) :: errStat !< Indicates whether an error occurred (see NWTC_Library) + character(*), intent(out) :: errMsg !< Error message associated with the errStat + integer(IntKi) :: maxMeshPts, iWT + INTEGER(IntKi) :: errStat2 + CHARACTER(ErrMsgLen) :: errMsg2 + character(*), parameter :: RoutineName = 'SetTempStorage' !< for error handling + ErrStat = ErrID_None + ErrMsg = "" + if (.not. allocated(NumMeshPts)) then + ErrStat = ErrID_Fatal + ErrMSg = "Pre-Init has not been called yet" + return + endif + if (minval(NumMeshPts) < 0) then + ErrStat = ErrID_Fatal + ErrMSg = "ADI_C_SetupRotor haven't been called for all rotors" + return + endif + + ! Allocate temporary arrays to simplify data conversions + maxMeshPts=0_IntKi + do iWT=1,Sim%NumTurbines + maxMeshPts=max(maxMeshPts,NumMeshPts(iWT)) + enddo + if ( allocated(tmpBldPtMeshPos) ) deallocate(tmpBldPtMeshPos) + if ( allocated(tmpBldPtMeshOri) ) deallocate(tmpBldPtMeshOri) + if ( allocated(tmpBldPtMeshVel) ) deallocate(tmpBldPtMeshVel) + if ( allocated(tmpBldPtMeshAcc) ) deallocate(tmpBldPtMeshAcc) + if ( allocated(tmpBldPtMeshFrc) ) deallocate(tmpBldPtMeshFrc) + call AllocAry( tmpBldPtMeshPos, 3, maxMeshPts, "tmpBldPtMeshPos", ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( tmpBldPtMeshOri, 3, 3, maxMeshPts, "tmpBldPtMeshOri", ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( tmpBldPtMeshVel, 6, maxMeshPts, "tmpBldPtMeshVel", ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( tmpBldPtMeshAcc, 6, maxMeshPts, "tmpBldPtMeshAcc", ErrStat2, ErrMsg2 ); if (Failed()) return + call AllocAry( tmpBldPtMeshFrc, 6, maxMeshPts, "tmpBldPtMeshFrc", ErrStat2, ErrMsg2 ); if (Failed()) return + +contains + logical function Failed() + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + end function Failed +end subroutine SetTempStorage + + +subroutine ClearTmpStorage() + INTEGER(IntKi) :: errStat2 + CHARACTER(ErrMsgLen) :: errMsg2 + ! Arrays + if (allocated(tmpBldPtMeshPos)) deallocate(tmpBldPtMeshPos) + if (allocated(tmpBldPtMeshOri)) deallocate(tmpBldPtMeshOri) + if (allocated(tmpBldPtMeshVel)) deallocate(tmpBldPtMeshVel) + if (allocated(tmpBldPtMeshAcc)) deallocate(tmpBldPtMeshAcc) + if (allocated(tmpBldPtMeshFrc)) deallocate(tmpBldPtMeshFrc) + ! Meshes + if (allocated(BldPtMotionMesh )) call ClearMeshArr1(BldPtMotionMesh ) + if (allocated(BldPtLoadMesh )) call ClearMeshArr1(BldPtLoadMesh ) + if (allocated(BldPtLoadMesh_tmp)) call ClearMeshArr1(BldPtLoadMesh_tmp) + !if (allocated(NacMotionMesh )) call ClearMeshArr1(NacMotionMesh ) + !if (allocated(NacLoadMesh )) call ClearMeshArr1(NacLoadMesh ) + if (allocated(Map_BldPtMotion_2_AD_Blade )) call ClearMeshMapArr2(Map_BldPtMotion_2_AD_Blade ) + if (allocated(Map_AD_BldLoad_P_2_BldPtLoad)) call ClearMeshMapArr2(Map_AD_BldLoad_P_2_BldPtLoad) +contains + !> Don't leave junk in memory. So destroy meshes and mappings. + subroutine ClearMeshArr1(MeshName) + type(MeshType), allocatable :: MeshName(:) + integer :: i + do i=1,size(MeshName) + call MeshDestroy( MeshName(i), ErrStat2, ErrMsg2 ) ! ignore errors + enddo + deallocate(MeshName) + end subroutine ClearMeshArr1 + + subroutine ClearMeshMapArr2(MapName) + type(MeshMapType), allocatable :: MapName(:,:) + integer :: i,j + do j=1,size(MapName,2) + do i=1,size(MapName,1) + call NWTC_Library_Destroymeshmaptype( MapName(i,j), ErrStat2, ErrMsg2 ) + enddo + enddo + deallocate(MapName) + end subroutine ClearMeshMapArr2 + +end subroutine ClearTmpStorage + END MODULE AeroDyn_Inflow_C_BINDING diff --git a/modules/aerodyn/src/AeroDyn_Registry.txt b/modules/aerodyn/src/AeroDyn_Registry.txt index 86142df842..618b6540d6 100644 --- a/modules/aerodyn/src/AeroDyn_Registry.txt +++ b/modules/aerodyn/src/AeroDyn_Registry.txt @@ -77,6 +77,7 @@ typedef ^ ^ SiKi TowerRad # ..... Initialization data ....................................................................................................... # Define inputs that the initialization routine may need here: typedef AeroDyn/AD RotInitInputType IntKi NumBlades - - - "Number of blades on the turbine" - +typedef ^ RotInitInputType ReKi originInit {3} - 0 "X-Y-Z reference position for the turbine" m typedef ^ RotInitInputType ReKi HubPosition {3} - - "X-Y-Z reference position of hub" m typedef ^ RotInitInputType R8Ki HubOrientation {3}{3} - - "DCM reference orientation of hub" - typedef ^ RotInitInputType ReKi BladeRootPosition {:}{:} - - "X-Y-Z reference position of each blade root (3 x NumBlades)" m diff --git a/modules/aerodyn/src/AeroDyn_Types.f90 b/modules/aerodyn/src/AeroDyn_Types.f90 index 0bbcb29b27..3bd6d76a50 100644 --- a/modules/aerodyn/src/AeroDyn_Types.f90 +++ b/modules/aerodyn/src/AeroDyn_Types.f90 @@ -94,6 +94,7 @@ MODULE AeroDyn_Types ! ========= RotInitInputType ======= TYPE, PUBLIC :: RotInitInputType INTEGER(IntKi) :: NumBlades = 0_IntKi !< Number of blades on the turbine [-] + REAL(ReKi) , DIMENSION(1:3) :: originInit = 0.0_ReKi !< X-Y-Z reference position for the turbine [m] REAL(ReKi) , DIMENSION(1:3) :: HubPosition = 0.0_ReKi !< X-Y-Z reference position of hub [m] REAL(R8Ki) , DIMENSION(1:3,1:3) :: HubOrientation = 0.0_R8Ki !< DCM reference orientation of hub [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: BladeRootPosition !< X-Y-Z reference position of each blade root (3 x NumBlades) [m] @@ -834,6 +835,7 @@ subroutine AD_CopyRotInitInputType(SrcRotInitInputTypeData, DstRotInitInputTypeD ErrStat = ErrID_None ErrMsg = '' DstRotInitInputTypeData%NumBlades = SrcRotInitInputTypeData%NumBlades + DstRotInitInputTypeData%originInit = SrcRotInitInputTypeData%originInit DstRotInitInputTypeData%HubPosition = SrcRotInitInputTypeData%HubPosition DstRotInitInputTypeData%HubOrientation = SrcRotInitInputTypeData%HubOrientation if (allocated(SrcRotInitInputTypeData%BladeRootPosition)) then @@ -888,6 +890,7 @@ subroutine AD_PackRotInitInputType(Buf, Indata) character(*), parameter :: RoutineName = 'AD_PackRotInitInputType' if (Buf%ErrStat >= AbortErrLev) return call RegPack(Buf, InData%NumBlades) + call RegPack(Buf, InData%originInit) call RegPack(Buf, InData%HubPosition) call RegPack(Buf, InData%HubOrientation) call RegPack(Buf, allocated(InData%BladeRootPosition)) @@ -918,6 +921,8 @@ subroutine AD_UnPackRotInitInputType(Buf, OutData) if (Buf%ErrStat /= ErrID_None) return call RegUnpack(Buf, OutData%NumBlades) if (RegCheckErr(Buf, RoutineName)) return + call RegUnpack(Buf, OutData%originInit) + if (RegCheckErr(Buf, RoutineName)) return call RegUnpack(Buf, OutData%HubPosition) if (RegCheckErr(Buf, RoutineName)) return call RegUnpack(Buf, OutData%HubOrientation) diff --git a/modules/aerodyn/src/FVW.f90 b/modules/aerodyn/src/FVW.f90 index 616d96ffc7..505efe7034 100644 --- a/modules/aerodyn/src/FVW.f90 +++ b/modules/aerodyn/src/FVW.f90 @@ -584,7 +584,10 @@ subroutine FVW_UpdateStates( t, n, u, utimes, p, x, xd, z, OtherState, AFInfo, m m%ComputeWakeInduced = .FALSE. endif if (bReevaluation) then - call WrScr('[INFO] FVW: Update States: reevaluation at the same starting time') + if (m%InfoReEval) then + call WrScr('[INFO] FVW: Update States: reevaluation at the same starting time. This will not print on subsequent occurences.') + m%InfoReEval = .false. + endif call RollBackPreviousTimeStep() ! Cancel wake emission done in previous call m%ComputeWakeInduced = .TRUE. endif diff --git a/modules/aerodyn/src/FVW_Registry.txt b/modules/aerodyn/src/FVW_Registry.txt index 698cbf0634..4b7b8028e9 100644 --- a/modules/aerodyn/src/FVW_Registry.txt +++ b/modules/aerodyn/src/FVW_Registry.txt @@ -209,6 +209,7 @@ typedef ^ ^ ReKi typedef ^ ^ ReKi Uind :: - - "Induced velocities obtained at control points" - # Outputs typedef ^ ^ GridOutType GridOutputs {:} - - "Number of VTK grid to output" - +typedef ^ ^ Logical InfoReeval - .true. - "Give info about Reevaluation: gets set to false after first info statement" - # ........ Input ............ # Rotors diff --git a/modules/aerodyn/src/FVW_Types.f90 b/modules/aerodyn/src/FVW_Types.f90 index b1210cdd91..537f0b06b5 100644 --- a/modules/aerodyn/src/FVW_Types.f90 +++ b/modules/aerodyn/src/FVW_Types.f90 @@ -240,6 +240,7 @@ MODULE FVW_Types REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: CPs !< Control points used for wake rollup computation [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: Uind !< Induced velocities obtained at control points [-] TYPE(GridOutType) , DIMENSION(:), ALLOCATABLE :: GridOutputs !< Number of VTK grid to output [-] + LOGICAL :: InfoReeval = .true. !< Give info about Reevaluation: gets set to false after first info statement [-] END TYPE FVW_MiscVarType ! ======================= ! ========= Rot_InputType ======= @@ -3531,6 +3532,7 @@ subroutine FVW_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) if (ErrStat >= AbortErrLev) return end do end if + DstMiscData%InfoReeval = SrcMiscData%InfoReeval end subroutine subroutine FVW_DestroyMisc(MiscData, ErrStat, ErrMsg) @@ -3639,6 +3641,7 @@ subroutine FVW_PackMisc(Buf, Indata) call FVW_PackGridOutType(Buf, InData%GridOutputs(i1)) end do end if + call RegPack(Buf, InData%InfoReeval) if (RegCheckErr(Buf, RoutineName)) return end subroutine @@ -3750,6 +3753,8 @@ subroutine FVW_UnPackMisc(Buf, OutData) call FVW_UnpackGridOutType(Buf, OutData%GridOutputs(i1)) ! GridOutputs end do end if + call RegUnpack(Buf, OutData%InfoReeval) + if (RegCheckErr(Buf, RoutineName)) return end subroutine subroutine FVW_CopyRot_InputType(SrcRot_InputTypeData, DstRot_InputTypeData, CtrlCode, ErrStat, ErrMsg) diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index 69fdba80da..e4c0085143 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -367,7 +367,9 @@ ad_regression("ad_BAR_OLAF" "aerodyn;bem") ad_regression("ad_BAR_SineMotion" "aerodyn;bem") ad_regression("ad_BAR_SineMotion_UA4_DBEMT3" "aerodyn;bem") ad_regression("ad_BAR_RNAMotion" "aerodyn;bem") +ad_regression("ad_B1n2_OLAF" "aerodyn;OLAF") py_ad_regression("py_ad_5MW_OC4Semi_WSt_WavesWN" "aerodyn;bem;python") +py_ad_regression("py_ad_B1n2_OLAF" "aerodyn;OLAF;python") # UnsteadyAero ua_regression("ua_redfreq" "unsteadyaero") diff --git a/reg_tests/r-test b/reg_tests/r-test index a2420b965b..de5b8f4e47 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit a2420b965b48f83db07235efba0aa1a7563e0b87 +Subproject commit de5b8f4e47207b9589fbd740631ae5bda2bb5b2f