Adding time to power output

erf-model · Oct 29, 2024 · 21f5569 · 21f5569
1 parent 1a95b9c
commit 21f5569
Show file tree

Hide file tree

Showing 14 changed files with 226 additions and 35 deletions.
diff --git a/Exec/GeneralActuatorDisk/ERF_prob.cpp b/Exec/GeneralActuatorDisk/ERF_prob.cpp
@@ -0,0 +1,164 @@
+#include "ERF_prob.H"
+#include "AMReX_Random.H"
+
+using namespace amrex;
+
+std::unique_ptr<ProblemBase>
+amrex_probinit(const amrex_real* problo, const amrex_real* probhi)
+{
+    return std::make_unique<Problem>(problo, probhi);
+}
+
+Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
+{
+  // Parse params
+  ParmParse pp("prob");
+  pp.query("rho_0", parms.rho_0);
+  pp.query("T_0", parms.T_0);
+  pp.query("A_0", parms.A_0);
+  pp.query("QKE_0", parms.QKE_0);
+
+  pp.query("U_0", parms.U_0);
+  pp.query("V_0", parms.V_0);
+  pp.query("W_0", parms.W_0);
+  pp.query("U_0_Pert_Mag", parms.U_0_Pert_Mag);
+  pp.query("V_0_Pert_Mag", parms.V_0_Pert_Mag);
+  pp.query("W_0_Pert_Mag", parms.W_0_Pert_Mag);
+  pp.query("T_0_Pert_Mag", parms.T_0_Pert_Mag);
+
+  pp.query("pert_deltaU", parms.pert_deltaU);
+  pp.query("pert_deltaV", parms.pert_deltaV);
+  pp.query("pert_periods_U", parms.pert_periods_U);
+  pp.query("pert_periods_V", parms.pert_periods_V);
+  pp.query("pert_ref_height", parms.pert_ref_height);
+  parms.aval = parms.pert_periods_U * 2.0 * PI / (probhi[1] - problo[1]);
+  parms.bval = parms.pert_periods_V * 2.0 * PI / (probhi[0] - problo[0]);
+  parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height;
+  parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height;
+
+  init_base_parms(parms.rho_0, parms.T_0);
+}
+
+void
+Problem::init_custom_pert(
+    const amrex::Box&  bx,
+    const amrex::Box& xbx,
+    const amrex::Box& ybx,
+    const amrex::Box& zbx,
+    amrex::Array4<amrex::Real const> const& /*state*/,
+    amrex::Array4<amrex::Real      > const& state_pert,
+    amrex::Array4<amrex::Real      > const& x_vel_pert,
+    amrex::Array4<amrex::Real      > const& y_vel_pert,
+    amrex::Array4<amrex::Real      > const& z_vel_pert,
+    amrex::Array4<amrex::Real      > const& /*r_hse*/,
+    amrex::Array4<amrex::Real      > const& /*p_hse*/,
+    amrex::Array4<amrex::Real const> const& /*z_nd*/,
+    amrex::Array4<amrex::Real const> const& /*z_cc*/,
+    amrex::GeometryData const& geomdata,
+    amrex::Array4<amrex::Real const> const& /*mf_m*/,
+    amrex::Array4<amrex::Real const> const& /*mf_u*/,
+    amrex::Array4<amrex::Real const> const& /*mf_v*/,
+    const SolverChoice& sc)
+{
+    const bool use_moisture = (sc.moisture_type != MoistureType::None);
+
+  ParallelForRNG(bx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept {
+    // Geometry
+    const Real* prob_lo = geomdata.ProbLo();
+    const Real* prob_hi = geomdata.ProbHi();
+    const Real* dx = geomdata.CellSize();
+    const Real x = prob_lo[0] + (i + 0.5) * dx[0];
+    const Real y = prob_lo[1] + (j + 0.5) * dx[1];
+    const Real z = prob_lo[2] + (k + 0.5) * dx[2];
+
+    // Define a point (xc,yc,zc) at the center of the domain
+    const Real xc = 0.5 * (prob_lo[0] + prob_hi[0]);
+    const Real yc = 0.5 * (prob_lo[1] + prob_hi[1]);
+    const Real zc = 0.5 * (prob_lo[2] + prob_hi[2]);
+
+    const Real r  = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc));
+
+    // Add temperature perturbations
+    if ((z <= parms_d.pert_ref_height) && (parms_d.T_0_Pert_Mag != 0.0)) {
+        Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
+        state_pert(i, j, k, RhoTheta_comp) = (rand_double*2.0 - 1.0)*parms_d.T_0_Pert_Mag;
+    }
+
+    // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain
+    state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0 * exp(-10.*r*r);
+
+    // Set an initial value for QKE
+    //state_pert(i, j, k, RhoQKE_comp) = parms_d.QKE_0;
+
+    if (use_moisture) {
+        state_pert(i, j, k, RhoQ1_comp) = 0.0;
+        state_pert(i, j, k, RhoQ2_comp) = 0.0;
+    }
+  });
+
+  // Set the x-velocity
+  ParallelForRNG(xbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept {
+    const Real* prob_lo = geomdata.ProbLo();
+    const Real* dx = geomdata.CellSize();
+    const Real y = prob_lo[1] + (j + 0.5) * dx[1];
+    const Real z = prob_lo[2] + (k + 0.5) * dx[2];
+
+    // Set the x-velocity
+    x_vel_pert(i, j, k) = parms_d.U_0;
+    if ((z <= parms_d.pert_ref_height) && (parms_d.U_0_Pert_Mag != 0.0))
+    {
+        Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
+        Real x_vel_prime = (rand_double*2.0 - 1.0)*parms_d.U_0_Pert_Mag;
+        x_vel_pert(i, j, k) += x_vel_prime;
+    }
+    if (parms_d.pert_deltaU != 0.0)
+    {
+        const amrex::Real yl = y - prob_lo[1];
+        const amrex::Real zl = z / parms_d.pert_ref_height;
+        const amrex::Real damp = std::exp(-0.5 * zl * zl);
+        x_vel_pert(i, j, k) += parms_d.ufac * damp * z * std::cos(parms_d.aval * yl);
+    }
+  });
+
+  // Set the y-velocity
+  ParallelForRNG(ybx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept {
+    const Real* prob_lo = geomdata.ProbLo();
+    const Real* dx = geomdata.CellSize();
+    const Real x = prob_lo[0] + (i + 0.5) * dx[0];
+    const Real z = prob_lo[2] + (k + 0.5) * dx[2];
+
+    // Set the y-velocity
+    y_vel_pert(i, j, k) = parms_d.V_0;
+    if ((z <= parms_d.pert_ref_height) && (parms_d.V_0_Pert_Mag != 0.0))
+    {
+        Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
+        Real y_vel_prime = (rand_double*2.0 - 1.0)*parms_d.V_0_Pert_Mag;
+        y_vel_pert(i, j, k) += y_vel_prime;
+    }
+    if (parms_d.pert_deltaV != 0.0)
+    {
+        const amrex::Real xl = x - prob_lo[0];
+        const amrex::Real zl = z / parms_d.pert_ref_height;
+        const amrex::Real damp = std::exp(-0.5 * zl * zl);
+        y_vel_pert(i, j, k) += parms_d.vfac * damp * z * std::cos(parms_d.bval * xl);
+    }
+  });
+
+  // Set the z-velocity
+  ParallelForRNG(zbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept {
+    const int dom_lo_z = geomdata.Domain().smallEnd()[2];
+    const int dom_hi_z = geomdata.Domain().bigEnd()[2];
+
+    // Set the z-velocity
+    if (k == dom_lo_z || k == dom_hi_z+1)
+    {
+        z_vel_pert(i, j, k) = 0.0;
+    }
+    else if (parms_d.W_0_Pert_Mag != 0.0)
+    {
+        Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
+        Real z_vel_prime = (rand_double*2.0 - 1.0)*parms_d.W_0_Pert_Mag;
+        z_vel_pert(i, j, k) = parms_d.W_0 + z_vel_prime;
+    }
+  });
+}
diff --git a/Source/ERF.H b/Source/ERF.H
@@ -435,7 +435,8 @@ public:
                        amrex::MultiFab& U_old, amrex::MultiFab& V_old, amrex::MultiFab& W_old,
                        amrex::MultiFab& mf_vars_windfarm,
                        const amrex::MultiFab& mf_Nturb,
-                       const amrex::MultiFab& mf_SMark);
+                       const amrex::MultiFab& mf_SMark,
+                       const amrex::Real& time);
 #endif
 
 #ifdef ERF_USE_EB

diff --git a/Source/Initialization/ERF_init_windfarm.cpp b/Source/Initialization/ERF_init_windfarm.cpp
@@ -59,8 +59,9 @@ ERF::advance_windfarm (const Geometry& a_geom,
                        MultiFab& W_old,
                        MultiFab& mf_vars_windfarm,
                        const MultiFab& mf_Nturb,
-                       const MultiFab& mf_SMark)
+                       const MultiFab& mf_SMark,
+                       const Real& time)
 {
         windfarm->advance(a_geom, dt_advance, cons_in, mf_vars_windfarm,
-                          U_old, V_old, W_old, mf_Nturb, mf_SMark);
+                          U_old, V_old, W_old, mf_Nturb, mf_SMark, time);
 }
diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp
@@ -110,7 +110,7 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/)
 #if defined(ERF_USE_WINDFARM)
     if (solverChoice.windfarm_type != WindFarmType::None) {
         advance_windfarm(Geom(lev), dt_lev, S_old,
-                         U_old, V_old, W_old, vars_windfarm[lev], Nturb[lev], SMark[lev]);
+                         U_old, V_old, W_old, vars_windfarm[lev], Nturb[lev], SMark[lev], time);
     }
 
 #endif

diff --git a/Source/WindFarmParametrization/ERF_WindFarm.H b/Source/WindFarmParametrization/ERF_WindFarm.H
@@ -90,10 +90,11 @@ public:
                   amrex::MultiFab& V_old,
                   amrex::MultiFab& W_old,
                   const amrex::MultiFab& mf_Nturb,
-                  const amrex::MultiFab& mf_SMark) override
+                  const amrex::MultiFab& mf_SMark,
+                  const amrex::Real& time) override
     {
         m_windfarm_model[0]->advance(a_geom, dt_advance, cons_in, mf_vars_windfarm,
-                                     U_old, V_old, W_old, mf_Nturb, mf_SMark);
+                                     U_old, V_old, W_old, mf_Nturb, mf_SMark, time);
     }
 
     void set_turb_spec(const amrex::Real& a_rotor_rad, const amrex::Real& a_hub_height,

diff --git a/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp b/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp
@@ -14,9 +14,11 @@ EWP::advance (const Geometry& geom,
               MultiFab& V_old,
               MultiFab& W_old,
               const MultiFab& mf_Nturb,
-              const MultiFab& mf_SMark)
+              const MultiFab& mf_SMark,
+              const Real& time)
  {
     AMREX_ALWAYS_ASSERT(mf_SMark.nComp() > 0);
+    AMREX_ALWAYS_ASSERT(time > -1.0);
     source_terms_cellcentered(geom, cons_in, mf_vars_ewp, U_old, V_old, W_old, mf_Nturb);
     update(dt_advance, cons_in, U_old, V_old, mf_vars_ewp);
 }

diff --git a/Source/WindFarmParametrization/EWP/ERF_EWP.H b/Source/WindFarmParametrization/EWP/ERF_EWP.H
@@ -22,7 +22,8 @@ public:
                   amrex::MultiFab& V_old,
                   amrex::MultiFab& W_old,
                   const amrex::MultiFab& mf_Nturb,
-                  const amrex::MultiFab& mf_SMark) override;
+                  const amrex::MultiFab& mf_SMark,
+                  const amrex::Real& time) override;
 
     void source_terms_cellcentered (const amrex::Geometry& geom,
                                     const amrex::MultiFab& cons_in,

diff --git a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp
@@ -53,10 +53,12 @@ Fitch::advance (const Geometry& geom,
                 MultiFab& V_old,
                 MultiFab& W_old,
                 const MultiFab& mf_Nturb,
-                const MultiFab& mf_SMark)
+                const MultiFab& mf_SMark,
+                const Real& time)
 {
     AMREX_ALWAYS_ASSERT(W_old.nComp() > 0);
     AMREX_ALWAYS_ASSERT(mf_SMark.nComp() > 0);
+    AMREX_ALWAYS_ASSERT(time > -1.0);
     source_terms_cellcentered(geom, cons_in, mf_vars_fitch, U_old, V_old, W_old, mf_Nturb);
     update(dt_advance, cons_in, U_old, V_old, mf_vars_fitch);
 }

diff --git a/Source/WindFarmParametrization/Fitch/ERF_Fitch.H b/Source/WindFarmParametrization/Fitch/ERF_Fitch.H
@@ -22,7 +22,8 @@ public:
                   amrex::MultiFab& V_old,
                   amrex::MultiFab& W_old,
                   const amrex::MultiFab& mf_Nturb,
-                  const amrex::MultiFab& mf_SMark) override;
+                  const amrex::MultiFab& mf_SMark,
+                  const amrex::Real& time) override;
 
     void source_terms_cellcentered (const amrex::Geometry& geom,
                                     const amrex::MultiFab& cons_in,

diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp
@@ -14,11 +14,13 @@ GeneralAD::advance (const Geometry& geom,
                   MultiFab& V_old,
                   MultiFab& W_old,
                   const MultiFab& mf_Nturb,
-                  const MultiFab& mf_SMark)
+                  const MultiFab& mf_SMark,
+                  const Real& time)
 {
     AMREX_ALWAYS_ASSERT(W_old.nComp() > 0);
     AMREX_ALWAYS_ASSERT(mf_Nturb.nComp() > 0);
     AMREX_ALWAYS_ASSERT(mf_vars_generalAD.nComp() > 0);
+    AMREX_ALWAYS_ASSERT(time > -1.0);
     compute_freestream_velocity(cons_in, U_old, V_old, mf_SMark);
     source_terms_cellcentered(geom, cons_in, mf_SMark, mf_vars_generalAD);
     update(dt_advance, cons_in, U_old, V_old, W_old, mf_vars_generalAD);

diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H
@@ -21,7 +21,8 @@ public:
                   amrex::MultiFab& V_old,
                   amrex::MultiFab& W_old,
                   const amrex::MultiFab& mf_Nturb,
-                  const amrex::MultiFab& mf_SMark) override;
+                  const amrex::MultiFab& mf_SMark,
+                  const amrex::Real& time) override;
 
     void compute_freestream_velocity (const amrex::MultiFab& cons_in,
                                      const amrex::MultiFab& U_old,

diff --git a/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H b/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H
@@ -22,7 +22,8 @@ public:
                   amrex::MultiFab& V_old,
                   amrex::MultiFab& W_old,
                   const amrex::MultiFab& mf_Nturb,
-                  const amrex::MultiFab& mf_SMark) = 0;
+                  const amrex::MultiFab& mf_SMark,
+                  const amrex::Real& time) = 0;
 
     virtual void set_turb_spec(const amrex::Real&  rotor_rad, const amrex::Real& hub_height,
                                const amrex::Real& thrust_coeff_standing, const amrex::Vector<amrex::Real>& wind_speed,

diff --git a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp
@@ -13,13 +13,44 @@ SimpleAD::advance (const Geometry& geom,
                   MultiFab& V_old,
                   MultiFab& W_old,
                   const MultiFab& mf_Nturb,
-                  const MultiFab& mf_SMark)
+                  const MultiFab& mf_SMark,
+                  const Real& time)
 {
     AMREX_ALWAYS_ASSERT(W_old.nComp() > 0);
     AMREX_ALWAYS_ASSERT(mf_Nturb.nComp() > 0);
     compute_freestream_velocity(cons_in, U_old, V_old, mf_SMark);
     source_terms_cellcentered(geom, cons_in, mf_SMark, mf_vars_simpleAD);
     update(dt_advance, cons_in, U_old, V_old, mf_vars_simpleAD);
+    compute_power_output(time);
+}
+
+void
+SimpleAD::compute_power_output(const Real& time)
+{
+     get_turb_loc(xloc, yloc);
+     get_turb_spec(rotor_rad, hub_height, thrust_coeff_standing,
+                  wind_speed, thrust_coeff, power);
+
+     const int n_spec_table = wind_speed.size();
+  // Compute power based on the look-up table
+
+    if (ParallelDescriptor::IOProcessor()){
+        static std::ofstream file("power_output.txt", std::ios::app);
+        // Check if the file opened successfully
+        if (!file.is_open()) {
+            std::cerr << "Error opening file!" << std::endl;
+            Abort("Could not open file to write power output in ERF_AdvanceSimpleAD.cpp");
+        }
+        Real total_power = 0.0;
+        for(int it=0; it<xloc.size(); it++){
+            Real avg_vel = freestream_velocity[it]/(disk_cell_count[it] + 1e-10);
+            Real turb_power = interpolate_1d(wind_speed.data(), power.data(), avg_vel, n_spec_table);
+            total_power = total_power + turb_power;
+            //printf("avg vel and power is %d %0.15g, %0.15g\n", it, avg_vel, turb_power);
+        }
+        file << time << " " << total_power << "\n";
+        file.flush();
+    }
 }
 
 void
@@ -225,24 +256,4 @@ SimpleAD::source_terms_cellcentered (const Geometry& geom,
             simpleAD_array(i,j,k,1) = source_y;
          });
     }
-
-  // Compute power based on the look-up table
-
-    if (ParallelDescriptor::IOProcessor()){
-        static std::ofstream file("power_output.txt", std::ios::app);
-        // Check if the file opened successfully
-        if (!file.is_open()) {
-            std::cerr << "Error opening file!" << std::endl;
-            Abort("Could not open file to write power output in ERF_AdvanceSimpleAD.cpp");
-        }
-        Real total_power = 0.0;
-        for(int it=0; it<xloc.size(); it++){
-            Real avg_vel = freestream_velocity[it]/(disk_cell_count[it] + 1e-10);
-            Real turb_power = interpolate_1d(wind_speed.data(), power.data(), avg_vel, n_spec_table);
-            total_power = total_power + turb_power;
-            //printf("avg vel and power is %d %0.15g, %0.15g\n", it, avg_vel, turb_power);
-        }
-        file << total_power << "\n";
-        file.flush();
-    }
 }
diff --git a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_SimpleAD.H b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_SimpleAD.H
@@ -21,7 +21,8 @@ public:
                   amrex::MultiFab& V_old,
                   amrex::MultiFab& W_old,
                   const amrex::MultiFab& mf_Nturb,
-                  const amrex::MultiFab& mf_SMark) override;
+                  const amrex::MultiFab& mf_SMark,
+                  const amrex::Real& time) override;
 
     void compute_freestream_velocity(const amrex::MultiFab& cons_in,
                                      const amrex::MultiFab& U_old,
@@ -39,6 +40,8 @@ public:
                  amrex::MultiFab& V_old,
                  const amrex::MultiFab& mf_vars);
 
+    void compute_power_output(const amrex::Real& time);
+
 protected:
     amrex::Vector<amrex::Real> xloc, yloc;
     amrex::Real turb_disk_angle;