From a4b7da7b38b6db71c3b57f69f4e48fd4dbb6713c Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Fri, 27 Sep 2024 09:04:44 -0700 Subject: [PATCH] Some changes to the simplified actuator disk model (#1830) * Generalize simplified actuator disk model * Removing unwatned files --------- Co-authored-by: Mahesh Natarajan Co-authored-by: Mahesh Natarajan Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> --- .../ParaViewExportSTLMacro_qc.py | 94 ---------------- .../ParaViewExportSTLMacro_qr.py | 93 ---------------- .../PostProcessing/SuperCell3D_Blender2p80.py | 101 ------------------ Source/DataStructs/ERF_DataStruct.H | 16 +-- Source/ERF.H | 4 +- .../Fitch/ERF_AdvanceFitch.cpp | 8 +- .../ERF_AdvanceSimpleAD.cpp | 29 ++++- 7 files changed, 41 insertions(+), 304 deletions(-) delete mode 100644 Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qc.py delete mode 100644 Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qr.py delete mode 100644 Exec/SuperCell_3D/PostProcessing/SuperCell3D_Blender2p80.py diff --git a/Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qc.py b/Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qc.py deleted file mode 100644 index 36b7644e0..000000000 --- a/Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qc.py +++ /dev/null @@ -1,94 +0,0 @@ -# trace generated using paraview version 5.12.0 -#import paraview -#paraview.compatibility.major = 5 -#paraview.compatibility.minor = 12 - -#### import the simple module from the paraview -from paraview.simple import * -from numpy import * -#### disable automatic camera reset on 'Show' -paraview.simple._DisableFirstRenderCameraReset() - -n = 167 -# create a new 'AMReX/BoxLib Grid Reader' -plot_filesseries = AMReXBoxLibGridReader(registrationName='plot_files.series', FileNames=['/pscratch/sd/n/nataraj2/ERF/SuperCell_3D/SolnFiles/4th_Standard_alpha66_YOpenXPer/plot_files.series']) - - -# get animation scene -animationScene1 = GetAnimationScene() - -# update animation scene based on data timesteps -animationScene1.UpdateAnimationUsingDataTimeSteps() - -# Properties modified on plot_filesseries -plot_filesseries.CellArrayStatus = ['qc', 'qrain', 'rain_accum'] - -# get active view -renderView1 = GetActiveViewOrCreate('RenderView') - -# show data in view -plot_filesseriesDisplay = Show(plot_filesseries, renderView1, 'AMRRepresentation') - -# trace defaults for the display properties. -plot_filesseriesDisplay.Representation = 'Outline' - -# reset view to fit data -renderView1.ResetCamera() - -# get the material library -materialLibrary1 = GetMaterialLibrary() - -# update the view to ensure updated data information -renderView1.Update() - -# create a new 'Cell Data to Point Data' -cellDatatoPointData1 = CellDatatoPointData(registrationName='CellDatatoPointData1', Input=plot_filesseries) - -# show data in view -cellDatatoPointData1Display = Show(cellDatatoPointData1, renderView1, 'AMRRepresentation') - -# trace defaults for the display properties. -cellDatatoPointData1Display.Representation = 'Outline' - -# hide data in view -Hide(plot_filesseries, renderView1) - -# update the view to ensure updated data information -renderView1.Update() - -# create a new 'Contour' -contour1 = Contour(registrationName='Contour1', Input=cellDatatoPointData1) - -# Properties modified on contour1 -contour1.Isosurfaces = [1e-05] - -# show data in view -contour1Display = Show(contour1, renderView1, 'GeometryRepresentation') - -# trace defaults for the display properties. -contour1Display.Representation = 'Surface' - -# update the view to ensure updated data information -renderView1.Update() - -# hide data in view -Hide(cellDatatoPointData1, renderView1) - -for i in arange(0, n, 1): - - print("Doing scene %d of %d"%(i,n)) - - - if(i<=9): - filename = "/pscratch/sd/n/nataraj2/ERF/SuperCell_3D/SolnFiles/4th_Standard_alpha66_YOpenXPer/supercell_3d_00%d.stl"%i - elif(i<=99): - filename = "/pscratch/sd/n/nataraj2/ERF/SuperCell_3D/SolnFiles/4th_Standard_alpha66_YOpenXPer/supercell_3d_0%d.stl"%i - elif(i<=999): - filename = "/pscratch/sd/n/nataraj2/ERF/SuperCell_3D/SolnFiles/4th_Standard_alpha66_YOpenXPer/supercell_3d_%d.stl"%i - - print("filename is %s", filename) - -# save data - SaveData(filename, proxy=contour1) - - animationScene1.GoToNext() diff --git a/Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qr.py b/Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qr.py deleted file mode 100644 index 763681d1e..000000000 --- a/Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qr.py +++ /dev/null @@ -1,93 +0,0 @@ -# trace generated using paraview version 5.12.0 -#import paraview -#paraview.compatibility.major = 5 -#paraview.compatibility.minor = 12 - -#### import the simple module from the paraview -from paraview.simple import * -from numpy import * -#### disable automatic camera reset on 'Show' -paraview.simple._DisableFirstRenderCameraReset() - -n = 167 -# create a new 'AMReX/BoxLib Grid Reader' -plot_filesseries = AMReXBoxLibGridReader(registrationName='plot_files.series', FileNames=['/pscratch/sd/n/nataraj2/ERF/SuperCell_3D/SolnFiles/4th_Standard_alpha66_YOpenXPer/plot_files.series']) - - -# get animation scene -animationScene1 = GetAnimationScene() - -# update animation scene based on data timesteps -animationScene1.UpdateAnimationUsingDataTimeSteps() - -# Properties modified on plot_filesseries -plot_filesseries.CellArrayStatus = ['qrain', 'rain_accum'] - -# get active view -renderView1 = GetActiveViewOrCreate('RenderView') - -# show data in view -plot_filesseriesDisplay = Show(plot_filesseries, renderView1, 'AMRRepresentation') - -# trace defaults for the display properties. -plot_filesseriesDisplay.Representation = 'Outline' - -# reset view to fit data -renderView1.ResetCamera() - -# get the material library -materialLibrary1 = GetMaterialLibrary() - -# update the view to ensure updated data information -renderView1.Update() - -# create a new 'Cell Data to Point Data' -cellDatatoPointData1 = CellDatatoPointData(registrationName='CellDatatoPointData1', Input=plot_filesseries) - -# show data in view -cellDatatoPointData1Display = Show(cellDatatoPointData1, renderView1, 'AMRRepresentation') - -# trace defaults for the display properties. -cellDatatoPointData1Display.Representation = 'Outline' - -# hide data in view -Hide(plot_filesseries, renderView1) - -# update the view to ensure updated data information -renderView1.Update() - -# create a new 'Contour' -contour1 = Contour(registrationName='Contour1', Input=cellDatatoPointData1) - -# Properties modified on contour1 -contour1.Isosurfaces = [1e-04] - -# show data in view -contour1Display = Show(contour1, renderView1, 'GeometryRepresentation') - -# trace defaults for the display properties. -contour1Display.Representation = 'Surface' - -# update the view to ensure updated data information -renderView1.Update() - -# hide data in view -Hide(cellDatatoPointData1, renderView1) - -for i in arange(0, n, 1): - - print("Doing scene %d of %d"%(i,n)) - - - if(i<=9): - filename = "/pscratch/sd/n/nataraj2/ERF/SuperCell_3D/SolnFiles/4th_Standard_alpha66_YOpenXPer/supercell_3d_00%d.stl"%i - elif(i<=99): - filename = "/pscratch/sd/n/nataraj2/ERF/SuperCell_3D/SolnFiles/4th_Standard_alpha66_YOpenXPer/supercell_3d_0%d.stl"%i - elif(i<=999): - filename = "/pscratch/sd/n/nataraj2/ERF/SuperCell_3D/SolnFiles/4th_Standard_alpha66_YOpenXPer/supercell_3d_%d.stl"%i - - print("filename is %s", filename) - -# save data - SaveData(filename, proxy=contour1) - diff --git a/Exec/SuperCell_3D/PostProcessing/SuperCell3D_Blender2p80.py b/Exec/SuperCell_3D/PostProcessing/SuperCell3D_Blender2p80.py deleted file mode 100644 index bf0b95478..000000000 --- a/Exec/SuperCell_3D/PostProcessing/SuperCell3D_Blender2p80.py +++ /dev/null @@ -1,101 +0,0 @@ -import bpy - -bpy.context.scene.render.engine = 'CYCLES' - -bpy.ops.object.delete(use_global=False) -#bpy.ops.wm.open_mainfile(filepath="supercell_3d_qc_only_trying_for_python.blend") -bpy.ops.wm.open_mainfile(filepath="supercell_3d_qc_only_straight_cam_for_python.blend") - -#bpy.context.scene.cycles.samples = 1024 # Adjust as needed -#bpy.context.scene.render.resolution_x = 8000 # Width -#bpy.context.scene.render.resolution_y = 4000 # Height -bpy.context.scene.render.resolution_percentage = 50 # Full resolution -bpy.context.scene.render.image_settings.file_format='PNG' -bpy.context.scene.render.image_settings.color_depth = '8' # 8 or 16 bits per channel -bpy.context.scene.render.image_settings.compression = 0 # Set to 15% compression -bpy.context.scene.render.image_settings.color_mode = 'RGBA' - -bpy.context.scene.render.use_border = True -bpy.context.scene.render.border_min_x = 0.0 -bpy.context.scene.render.border_max_x = 1.0 -bpy.context.scene.render.border_min_y = 0.15 -bpy.context.scene.render.border_max_y = 0.85 - -for i in range(5,167,1): - - bpy.ops.object.text_add(location=(0, 0, 0)) - text_obj = bpy.context.object - text_obj.data.body = "Time (hr): %f"%(i/60.0) - text_obj.location = (-30, 0, 15) - text_obj.data.size = 1.0 - text_obj.scale = (2.0, 2.0, 2.0) # Scale uniformly in all directions - text_obj.rotation_euler = (90, 0, 0) # Replace with desired rotation - text_obj.data.materials.clear() - - if(i<=9): - qc_filename = 'supercell_3d_00%d0.stl'%i - elif(i<=99): - qc_filename = 'supercell_3d_0%d0.stl'%i - elif(i<=999): - qc_filename = 'supercell_3d_%d0.stl'%i - - file_loc = './STLFiles/qc/%s'%qc_filename - - if(i<=9): - img_filename = 'supercell_3d_00%d.jpg'%i - elif(i<=99): - img_filename = 'supercell_3d_0%d.jpg'%i - elif(i<=999): - img_filename = 'supercell_3d_%d.jpg'%i - - print("filename: ",file_loc); - #drop=bpy.ops.import_mesh.stl(filepath=file_loc) - drop=bpy.ops.wm.stl_import(filepath=file_loc) - mat_drop = bpy.data.materials.get("Drop") - if mat_drop is None: - # create material - mat_drop = bpy.data.materials.new(name="Drop") - - mat_drop.use_nodes=True - nodes=mat_drop.node_tree.nodes - for node in nodes: - nodes.remove(node) - - diffuseshader = nodes.new(type='ShaderNodeBsdfDiffuse') - diffuseshader.inputs[1].default_value=0.0 - diffuseshader.inputs['Color'].default_value = (0.846, 0.846, 0.846, 1.0) # RGBA - - - #glossyshader = nodes.new(type='ShaderNodeBsdfGlossy') - #glossyshader.inputs[1].default_value=0.0 - mixshader = nodes.new(type='ShaderNodeMixShader') - mixshader.inputs['Fac'].default_value = 0.1 - dropoutput = nodes.new(type='ShaderNodeOutputMaterial') - - links = mat_drop.node_tree.links - link = links.new(diffuseshader.outputs[0], mixshader.inputs[1]) - #link = links.new(glossyshader.outputs[0], mixshader.inputs[2]) - link = links.new(mixshader.outputs[0], dropoutput.inputs[0]) - - # Get material - ob = bpy.context.active_object - if ob.data.materials: - # assign to 1st material slot - ob.data.materials[1] = mat_drop - else: - # no slots - ob.data.materials.append(mat_drop) - bpy.data.objects[7].scale = (0.001, 0.001, 0.001) - - for area in bpy.context.screen.areas: - if area.type == 'VIEW_3D': - area.spaces[0].shading.type = 'RENDERED' - - bpy.context.scene.render.filepath = "./Images/%s"%img_filename - bpy.ops.render.render(write_still=True) - - bpy.ops.object.select_all(action='DESELECT') - bpy.data.objects[7].select_set(True) - bpy.ops.object.delete() - text_obj.select_set(True) # Select the object - bpy.ops.object.delete() # Delete the selected object diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index e78e668fe..a4e4616f9 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -402,17 +402,17 @@ struct SolverChoice { // by the diameter of the turbine pp.query("sampling_distance_by_D", sampling_distance_by_D); if(windfarm_type==WindFarmType::SimpleAD and sampling_distance_by_D < 0.0) { - amrex::Abort("To use simplified actuator disks, you need to provide a variable " - "erf.sampling_distance_by_D in the inputs which specifies the upstream " - "distance as a factor of the turbine diameter at which the incoming free stream " - "velocity will be computed at"); + amrex::Abort("To use simplified actuator disks, you need to provide a variable" + " erf.sampling_distance_by_D in the inputs which specifies the upstream" + " distance as a factor of the turbine diameter at which the incoming free stream" + " velocity will be computed at."); } pp.query("turb_disk_angle_from_x", turb_disk_angle); if(windfarm_type==WindFarmType::SimpleAD and turb_disk_angle < 0.0) { - amrex::Abort("To use simplified actuator disks, you need to provide a variable " - " erf.turb_disk_angle_from_x in the inputs which is the angle of the face of the " - " turbine disk from the x-axis. A turbine facing an oncoming flow in the x-direction " - " will have turb_disk_angle value of 90 deg "); + amrex::Abort("To use simplified actuator disks, you need to provide a variable" + " erf.turb_disk_angle_from_x in the inputs which is the angle of the face of the" + " turbine disk from the x-axis. A turbine facing an oncoming flow in the x-direction" + " will have turb_disk_angle value of 90 deg."); } pp.query("windfarm_x_shift",windfarm_x_shift); diff --git a/Source/ERF.H b/Source/ERF.H index 884412fe9..386649e5f 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -411,7 +411,9 @@ public: const amrex::Real& dt_advance, amrex::MultiFab& cons_in, amrex::MultiFab& U_old, amrex::MultiFab& V_old, amrex::MultiFab& W_old, - amrex::MultiFab& mf_vars_windfarm, const amrex::MultiFab& mf_Nturb); + amrex::MultiFab& mf_vars_windfarm, + const amrex::MultiFab& mf_Nturb, + const amrex::MultiFab& mf_SMark); #endif #ifdef ERF_USE_EB diff --git a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp index 41d0d45ff..7770ea32d 100644 --- a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp +++ b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp @@ -124,12 +124,12 @@ Fitch::source_terms_cellcentered (const Geometry& geom, mf_vars_fitch.setVal(0.0); Real d_hub_height = hub_height; Real d_rotor_rad = rotor_rad; - Gpu::DeviceVector d_wind_speed(wind_speed.size()); - Gpu::DeviceVector d_thrust_coeff(thrust_coeff.size()); + Gpu::DeviceVector d_wind_speed(wind_speed.size()); + Gpu::DeviceVector d_thrust_coeff(thrust_coeff.size()); // Copy data from host vectors to device vectors - Gpu::copy(Gpu::hostToDevice, wind_speed.begin(), wind_speed.end(), d_wind_speed.begin()); - Gpu::copy(Gpu::hostToDevice, thrust_coeff.begin(), thrust_coeff.end(), d_thrust_coeff.begin()); + Gpu::copy(Gpu::hostToDevice, wind_speed.begin(), wind_speed.end(), d_wind_speed.begin()); + Gpu::copy(Gpu::hostToDevice, thrust_coeff.begin(), thrust_coeff.end(), d_thrust_coeff.begin()); for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) { diff --git a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp index 486d174fa..8aa569a03 100644 --- a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp +++ b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp @@ -1,5 +1,6 @@ #include #include +#include using namespace amrex; @@ -173,6 +174,17 @@ SimpleAD::source_terms_cellcentered (const Geometry& geom, Real ny = -std::sin(turb_disk_angle); Real d_turb_disk_angle = turb_disk_angle; + Gpu::DeviceVector d_wind_speed(wind_speed.size()); + Gpu::DeviceVector d_thrust_coeff(thrust_coeff.size()); + + // Copy data from host vectors to device vectors + Gpu::copy(Gpu::hostToDevice, wind_speed.begin(), wind_speed.end(), d_wind_speed.begin()); + Gpu::copy(Gpu::hostToDevice, thrust_coeff.begin(), thrust_coeff.end(), d_thrust_coeff.begin()); + + const Real* wind_speed_d = d_wind_speed.dataPtr(); + const Real* thrust_coeff_d = d_thrust_coeff.dataPtr(); + const int n_spec_table = d_wind_speed.size(); + for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& gbx = mfi.growntilebox(1); @@ -193,11 +205,22 @@ SimpleAD::source_terms_cellcentered (const Geometry& geom, Real avg_vel = d_freestream_velocity_ptr[it]/(d_disk_cell_count_ptr[it] + 1e-10); Real phi = d_freestream_phi_ptr[it]/(d_disk_cell_count_ptr[it] + 1e-10); + Real C_T = interpolate_1d(wind_speed_d, thrust_coeff_d, avg_vel, n_spec_table); + Real a; + if(C_T <= 1) { + a = 0.5 - 0.5*std::pow(1.0-C_T,0.5); + } Real Uinfty_dot_nhat = avg_vel*(std::cos(phi)*nx + std::sin(phi)*ny); if(SMark_array(ii,jj,kk,1) == static_cast(it)) { - check_int++; - source_x = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*0.25*(1.0-0.25)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi); - source_y = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*0.25*(1.0-0.25)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi); + check_int++; + if(C_T <= 1) { + source_x = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*a*(1.0-a)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi); + source_y = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*a*(1.0-a)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi); + } + else { + source_x = -0.5*C_T*std::pow(Uinfty_dot_nhat, 2.0)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi); + source_y = -0.5*C_T*std::pow(Uinfty_dot_nhat, 2.0)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi); + } } } if(check_int > 1){