Sampling velocities upstream of the turbine for actuator disks (erf-m…

…odel#1751)

* Sampling velocities upstream of the turbine for actuator disks

* Correcting vtk file write

* Fill SMark only for actuator disks

* Updating AMReX submodule to latest

---------

Co-authored-by: Mahesh Natarajan <[email protected]>
Co-authored-by: Aaron M. Lattanzi <[email protected]>
Co-authored-by: Mahesh Natarajan <[email protected]>
Co-authored-by: Mahesh Natarajan <[email protected]>

Source/DataStructs/DataStruct.H

@@ -374,6 +374,18 @@ struct SolverChoice {
         pp.query("windfarm_loc_table",  windfarm_loc_table);
         pp.query("windfarm_spec_table", windfarm_spec_table);
 
+        // Sampling distance upstream of the turbine to find the
+        // incoming free stream velocity as a factor of the diameter of the
+        // turbine. ie. the sampling distance will be this number multiplied
+        // 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");
+        }
+
         // Test if time averaged data is to be output
         pp.query("time_avg_vel",time_avg_vel);
     }
@@ -578,5 +590,6 @@ struct SolverChoice {
 
     amrex::Real latitude_lo=-1e10, longitude_lo=-1e10;
     std::string windfarm_loc_table, windfarm_spec_table;
+    amrex::Real sampling_distance_by_D=-1.0;
 };
 #endif

Source/ERF.H

@@ -693,6 +693,10 @@ private:
     amrex::Vector<amrex::MultiFab> Nturb;
     amrex::Vector<amrex::MultiFab> vars_windfarm; // Fitch: Vabs, Vabsdt, dudt, dvdt, dTKEdt
                                                   // EWP:                 dudt, dvdt, dTKEdt
+
+    amrex::Vector<amrex::MultiFab> SMark; // A multifab that holds an integer corresponding
+                                                   // to the number of the wind turbine to sample
+                                                   // velocity upstream of the turbine
 #endif
 
     LandSurface lsm;
@@ -867,7 +871,7 @@ private:
                                                     "vorticity_x","vorticity_y","vorticity_z",
                                                     "magvel", "divU",
                                                     "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens",
-                                                    "eq_pot_temp", "num_turb",
+                                                    "eq_pot_temp", "num_turb", "SMark",
                                                     "dpdx", "dpdy", "pres_hse_x", "pres_hse_y",
                                                     "z_phys", "detJ" , "mapfac", "lat_m", "lon_m",
                                                     // Time averaged velocity

Source/ERF.cpp

@@ -116,6 +116,7 @@ ERF::ERF ()
 #ifdef ERF_USE_WINDFARM
     Nturb.resize(nlevs_max);
     vars_windfarm.resize(nlevs_max);
+    SMark.resize(nlevs_max);
 #endif
 
 #if defined(ERF_USE_RRTMGP)

Source/ERF_make_new_arrays.cpp

@@ -255,6 +255,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
     if (solverChoice.windfarm_type == WindFarmType::SimpleAD) {
         vars_windfarm[lev].define(ba, dm, 2, ngrow_state);// dudt, dvdt
         Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell
+        SMark[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell
     }
 #endif
 

Source/IO/Plotfile.cpp

@@ -447,6 +447,19 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             }
             mf_comp ++;
         }
+
+        if(containerHasElement(plot_var_names, "SMark") and solverChoice.windfarm_type == WindFarmType::SimpleAD) {
+             for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
+            {
+                const Box& bx = mfi.tilebox();
+                const Array4<Real>& derdat  = mf[lev].array(mfi);
+                const Array4<Real const>& SMark_array = SMark[lev].const_array(mfi);
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                    derdat(i, j, k, mf_comp) = SMark_array(i,j,k,0);
+                });
+            }
+            mf_comp ++;
+        }
 #endif
 
         int klo = geom[lev].Domain().smallEnd(2);

Source/Initialization/ERF_init_windfarm.cpp

@@ -33,8 +33,10 @@ ERF::init_windfarm (int lev)
     windfarm->write_turbine_locations_vtk();
 
     if(solverChoice.windfarm_type == WindFarmType::SimpleAD) {
-        windfarm->write_actuator_disks_vtk();
+        windfarm->fill_SMark_multifab(geom[lev], SMark[lev], solverChoice.sampling_distance_by_D);
+        windfarm->write_actuator_disks_vtk(geom[lev]);
     }
+
 }
 
 void

Source/WindFarmParametrization/InitWindFarm.cpp

@@ -245,6 +245,58 @@ WindFarm::fill_Nturb_multifab(const Geometry& geom,
 }
 
 
+void
+WindFarm::fill_SMark_multifab(const Geometry& geom,
+                              MultiFab& mf_SMark,
+                              const Real& sampling_distance_by_D)
+{
+    amrex::Gpu::DeviceVector<Real> d_xloc(xloc.size());
+    amrex::Gpu::DeviceVector<Real> d_yloc(yloc.size());
+    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, xloc.begin(), xloc.end(), d_xloc.begin());
+    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, yloc.begin(), yloc.end(), d_yloc.begin());
+
+    Real d_rotor_rad = rotor_rad;
+    Real d_hub_height = hub_height;
+    Real d_sampling_distance = sampling_distance_by_D*2.0*rotor_rad;
+
+    Real* d_xloc_ptr     = d_xloc.data();
+    Real* d_yloc_ptr     = d_yloc.data();
+
+    mf_SMark.setVal(0);
+
+    int i_lo = geom.Domain().smallEnd(0); int i_hi = geom.Domain().bigEnd(0);
+    int j_lo = geom.Domain().smallEnd(1); int j_hi = geom.Domain().bigEnd(1);
+    int k_lo = geom.Domain().smallEnd(2); int k_hi = geom.Domain().bigEnd(2);
+    auto dx = geom.CellSizeArray();
+    auto ProbLoArr = geom.ProbLoArray();
+    int num_turb = xloc.size();
+
+     // Initialize wind farm
+    for ( MFIter mfi(mf_SMark,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+        const Box& bx     = mfi.tilebox();
+        auto  SMark_array = mf_SMark.array(mfi);
+        ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+            int ii = amrex::min(amrex::max(i, i_lo), i_hi);
+            int jj = amrex::min(amrex::max(j, j_lo), j_hi);
+            int kk = amrex::min(amrex::max(k, k_lo), k_hi);
+
+            Real x1 = ProbLoArr[0] + ii*dx[0];
+            Real x2 = ProbLoArr[0] + (ii+1)*dx[0];
+
+            Real y = ProbLoArr[1] + (jj+0.5) * dx[1];
+            Real z = ProbLoArr[2] + (kk+0.5) * dx[2];
+
+            for(int it=0; it<num_turb; it++){
+                if(d_xloc_ptr[it]-d_sampling_distance+1e-12 > x1 and d_xloc_ptr[it]-d_sampling_distance+1e-12 < x2) {
+                   if(std::pow((y-d_yloc_ptr[it])*(y-d_yloc_ptr[it]) + (z-d_hub_height)*(z-d_hub_height),0.5) < d_rotor_rad) {
+                       SMark_array(i,j,k,0) = it;
+                    }
+                }
+            }
+        });
+    }
+}
+
 void
 WindFarm::write_turbine_locations_vtk()
 {
@@ -265,37 +317,75 @@ WindFarm::write_turbine_locations_vtk()
 
 
 void
-WindFarm::write_actuator_disks_vtk()
+WindFarm::write_actuator_disks_vtk(const Geometry& geom)
 {
+
     if (ParallelDescriptor::IOProcessor()){
-        FILE* file_actuator_disks;
-        file_actuator_disks = fopen("actuator_disks.vtk","w");
-        fprintf(file_actuator_disks, "%s\n","# vtk DataFile Version 3.0");
-        fprintf(file_actuator_disks, "%s\n","Actuator Disks");
-        fprintf(file_actuator_disks, "%s\n","ASCII");
-        fprintf(file_actuator_disks, "%s\n","DATASET POLYDATA");
+        FILE *file_actuator_disks_all, *file_actuator_disks_in_dom;
+        file_actuator_disks_all = fopen("actuator_disks_all.vtk","w");
+        fprintf(file_actuator_disks_all, "%s\n","# vtk DataFile Version 3.0");
+        fprintf(file_actuator_disks_all, "%s\n","Actuator Disks");
+        fprintf(file_actuator_disks_all, "%s\n","ASCII");
+        fprintf(file_actuator_disks_all, "%s\n","DATASET POLYDATA");
+
+        file_actuator_disks_in_dom = fopen("actuator_disks_in_dom.vtk","w");
+        fprintf(file_actuator_disks_in_dom, "%s\n","# vtk DataFile Version 3.0");
+        fprintf(file_actuator_disks_in_dom, "%s\n","Actuator Disks");
+        fprintf(file_actuator_disks_in_dom, "%s\n","ASCII");
+        fprintf(file_actuator_disks_in_dom, "%s\n","DATASET POLYDATA");
+
         int npts = 100;
-        fprintf(file_actuator_disks, "%s %ld %s\n", "POINTS", xloc.size()*npts, "float");
+        fprintf(file_actuator_disks_all, "%s %ld %s\n", "POINTS", xloc.size()*npts, "float");
+        auto ProbLoArr = geom.ProbLoArray();
+        auto ProbHiArr = geom.ProbHiArray();
+        int num_turb_in_dom = 0;
+
+        // Find the number of turbines inside the specified computational domain
+
         for(int it=0; it<xloc.size(); it++){
+            Real x = xloc[it];
+            Real y = yloc[it];
+            if(x > ProbLoArr[0] and x < ProbHiArr[0] and y > ProbLoArr[1] and y < ProbHiArr[1]) {
+                num_turb_in_dom++;
+            }
+        }
+        fprintf(file_actuator_disks_in_dom, "%s %ld %s\n", "POINTS", num_turb_in_dom*npts, "float");
+
+        for(int it=0; it<xloc.size(); it++){
+            Real x;
+            x = xloc[it]+1e-12;
             for(int pt=0;pt<100;pt++){
-                Real x, y, z;
+                Real y, z;
                 Real theta = 2.0*M_PI/npts*pt;
-                x = xloc[it]+1e-12;
                 y = yloc[it]+rotor_rad*cos(theta);
                 z = hub_height+rotor_rad*sin(theta);
-                fprintf(file_actuator_disks, "%0.15g %0.15g %0.15g\n", x, y, z);
+                fprintf(file_actuator_disks_all, "%0.15g %0.15g %0.15g\n", x, y, z);
+                if(xloc[it] > ProbLoArr[0] and xloc[it] < ProbHiArr[0] and yloc[it] > ProbLoArr[1] and yloc[it] < ProbHiArr[1]) {
+                    fprintf(file_actuator_disks_in_dom, "%0.15g %0.15g %0.15g\n", x, y, z);
+                }
             }
         }
-        fprintf(file_actuator_disks, "%s %ld %ld\n", "LINES", xloc.size()*(npts-1), static_cast<long int>(xloc.size()*(npts-1)*3));
+        fprintf(file_actuator_disks_all, "%s %ld %ld\n", "LINES", xloc.size()*(npts-1), static_cast<long int>(xloc.size()*(npts-1)*3));
+        fprintf(file_actuator_disks_in_dom, "%s %ld %ld\n", "LINES", num_turb_in_dom*(npts-1), static_cast<long int>(num_turb_in_dom*(npts-1)*3));
         for(int it=0; it<xloc.size(); it++){
             for(int pt=0;pt<99;pt++){
-                fprintf(file_actuator_disks, "%ld %ld %ld\n",
+                fprintf(file_actuator_disks_all, "%ld %ld %ld\n",
                                              static_cast<long int>(2),
                                              static_cast<long int>(it*npts+pt),
                                              static_cast<long int>(it*npts+pt+1));
             }
         }
-        fclose(file_actuator_disks);
+         for(int it=0; it<num_turb_in_dom; it++){
+            for(int pt=0;pt<99;pt++){
+                fprintf(file_actuator_disks_in_dom, "%ld %ld %ld\n",
+                                             static_cast<long int>(2),
+                                             static_cast<long int>(it*npts+pt),
+                                             static_cast<long int>(it*npts+pt+1));
+            }
+        }
+
+        fclose(file_actuator_disks_all);
+        fclose(file_actuator_disks_in_dom);
     }
 }
 

Source/WindFarmParametrization/WindFarm.H

@@ -55,11 +55,16 @@ public:
 
     void read_windfarm_spec_table(const std::string windfarm_spec_table);
 
-    void fill_Nturb_multifab(const amrex::Geometry& geom, amrex::MultiFab& mf_Nturb);
+    void fill_Nturb_multifab(const amrex::Geometry& geom,
+                             amrex::MultiFab& mf_Nturb);
+
+    void fill_SMark_multifab(const amrex::Geometry& geom,
+                             amrex::MultiFab& mf_SMark,
+                             const amrex::Real& sampling_distance_by_D);
 
     void write_turbine_locations_vtk();
 
-    void write_actuator_disks_vtk();
+    void write_actuator_disks_vtk(const amrex::Geometry& geom);
 
     void advance (const amrex::Geometry& a_geom,
                   const amrex::Real& dt_advance,