Skip to content

Commit

Permalink
Merge pull request #336 from hklion/drag
Browse files Browse the repository at this point in the history 
Add options for logarithmic and quadratic drag
  • Loading branch information
@hklion
hklion authored Jan 23, 2025
2 parents 2c19cee + db21f09 commit 1b99d10
Show file tree
Hide file tree
Showing 21 changed files with 677 additions and 17 deletions.
36 changes: 34 additions & 2 deletions Docs/sphinx_doc/Inputs.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -588,8 +588,40 @@ List of Parameters
+-----------------------------------+----------------------------------------+-------------------+----------------+
| **remora.Akt_bak** | Minimum/initial value of Akt | Real number | 1.0e-6 |
+-----------------------------------+----------------------------------------+-------------------+----------------+
| **remora.rdrag** | Bottom drag | Real number | 3.0e-4 |
+-----------------------------------+----------------------------------------+-------------------+----------------+

.. _list-of-parameters-drag:

List of drag-related parameters
--------------------------------

+-----------------------------------+------------------------------------------+-------------------+----------------+
| Parameter | Definition | Acceptable | Default |
| | | Values | |
+===================================+==========================================+===================+================+
| **remora.bottom_stress_type** | Bottom drag type | linear | linear |
| | | / quadratic | |
| | | / logarithmic | |
+-----------------------------------+------------------------------------------+-------------------+----------------+
| **remora.rdrag** | Linear drag coefficient (used if | Positive real | 3.0e-4 |
| | remora.bottom_stress_type = linear) | number | |
+-----------------------------------+------------------------------------------+-------------------+----------------+
| **remora.rdrag2** | Quadratic drag coefficient (used if | Positive real | 3.0e-3 |
| | remora.bottom_stress_type = quadratic) | number | |
+-----------------------------------+------------------------------------------+-------------------+----------------+
| **remora.Zob** | Bottom roughness [m] (used if | Positive real | 0.02 |
| | remora.bottom_stress_type = logarithmic | number | |
| | and for GLS | | |
+-----------------------------------+------------------------------------------+-------------------+----------------+
| **remora.Zos** | Surface roughness [m] (used for GLS) | Positive real | 0.02 |
| | | number | |
+-----------------------------------+------------------------------------------+-------------------+----------------+
| **remora.Cdb_min** | Minimum threshold for transfer | Positive real | 1.0e-6 |
| | coefficient of momentum | number | |
+-----------------------------------+------------------------------------------+-------------------+----------------+
| **remora.Cdb_max** | Maximum threshold for transfer | Positive real | 0.5 |
| | coefficient of momentum | number | |
+-----------------------------------+------------------------------------------+-------------------+----------------+


.. _list-of-parameters-gls:

Expand Down
27 changes: 24 additions & 3 deletions Source/Initialization/REMORA_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -324,6 +324,8 @@ void REMORA::resize_stuff(int lev)
vec_sustr.resize(lev+1);
vec_svstr.resize(lev+1);
vec_rdrag.resize(lev+1);
vec_rdrag2.resize(lev+1);
vec_ZoBot.resize(lev+1);
vec_bustr.resize(lev+1);
vec_bvstr.resize(lev+1);

Expand Down Expand Up @@ -451,8 +453,17 @@ void REMORA::init_stuff (int lev, const BoxArray& ba, const DistributionMapping&
vec_sustr[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0))); //2d, surface stress
vec_svstr[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0))); //2d

//2d, linear drag coefficient [m/s], defined at rho, somehow related to rdrg
vec_rdrag[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
if (solverChoice.bottom_stress_type == BottomStressType::linear) {
//2d, linear drag coefficient [m/s], defined at rho, somehow related to rdrg
vec_rdrag[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
} else if (solverChoice.bottom_stress_type == BottomStressType::quadratic) {
vec_rdrag2[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
}

if (solverChoice.bottom_stress_type == BottomStressType::logarithmic ||
solverChoice.vert_mixing_type == VertMixingType::GLS) {
vec_ZoBot[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
}

vec_bustr[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0))); //2d, bottom stress
vec_bvstr[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0)));
Expand Down Expand Up @@ -524,7 +535,17 @@ void REMORA::init_stuff (int lev, const BoxArray& ba, const DistributionMapping&
// when init_type = real. However, this does not appear to be necessary so removing

// Set initial linear drag coefficient
vec_rdrag[lev]->setVal(solverChoice.rdrag);
if (solverChoice.bottom_stress_type == BottomStressType::linear) {
vec_rdrag[lev]->setVal(solverChoice.rdrag);
} else if (solverChoice.bottom_stress_type == BottomStressType::quadratic) {
vec_rdrag2[lev]->setVal(solverChoice.rdrag2);
}

if (solverChoice.bottom_stress_type == BottomStressType::logarithmic ||
solverChoice.vert_mixing_type == VertMixingType::GLS) {
vec_ZoBot[lev]->setVal(solverChoice.Zob);
}


// ********************************************************************************************
// Create the REMORAFillPatcher object
Expand Down
4 changes: 4 additions & 0 deletions Source/REMORA.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -263,6 +263,10 @@ public:

/** Linear drag coefficient [m/s], defined at rho points */
amrex::Vector<std::unique_ptr<amrex::MultiFab>> vec_rdrag;
/** Quadratic drag coefficient [unitless], defined at rho points */
amrex::Vector<std::unique_ptr<amrex::MultiFab>> vec_rdrag2;
/** Bottom roughness length [m], defined at rho points */
amrex::Vector<std::unique_ptr<amrex::MultiFab>> vec_ZoBot;

/** Bottom stress in the u direction */
amrex::Vector<std::unique_ptr<amrex::MultiFab>> vec_bustr;
Expand Down
31 changes: 29 additions & 2 deletions Source/REMORA_DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,10 @@ enum class EOSType {
linear, nonlinear
};

enum class BottomStressType {
linear, quadratic, logarithmic
};

struct SolverChoice {
public:
void init_params()
Expand Down Expand Up @@ -84,8 +88,11 @@ struct SolverChoice {
amrex::Error("UV advection scheme unknown.");

pp.query("rdrag", rdrag);
pp.query("rdrag2", rdrag2);
pp.query("Zos", Zos);
pp.query("Zob", Zob);
pp.query("Cdb_max", Cdb_max);
pp.query("Cdb_min", Cdb_min);

// Order of spatial discretization
pp.query("spatial_order", spatial_order);
Expand Down Expand Up @@ -190,6 +197,18 @@ struct SolverChoice {
}
}

static std::string bottom_stress_type_string = "linear";
pp.query("bottom_stress_type", bottom_stress_type_string);
if (amrex::toLower(bottom_stress_type_string) == "linear") {
bottom_stress_type = BottomStressType::linear;
} else if (amrex::toLower(bottom_stress_type_string) == "quadratic") {
bottom_stress_type = BottomStressType::quadratic;
} else if (amrex::toLower(bottom_stress_type_string) == "logarithmic") {
bottom_stress_type = BottomStressType::logarithmic;
} else {
amrex::Abort("Don't know this bottom_stress_type");
}

amrex::Real tnu2_salt = amrex::Real(0.0);
amrex::Real tnu2_temp = amrex::Real(0.0);
amrex::Real tnu2_scalar = amrex::Real(0.0);
Expand Down Expand Up @@ -401,6 +420,9 @@ struct SolverChoice {
// EOS type
EOSType eos_type;

// Bottom stress type
BottomStressType bottom_stress_type;

// Mixing type and parameters
VertMixingType vert_mixing_type;
HorizMixingType horiz_mixing_type;
Expand All @@ -411,13 +433,18 @@ struct SolverChoice {
amrex::Real theta_b = amrex::Real(0.0);
amrex::Real tcline = amrex::Real(150.0);

// Linear drag coefficient
// Linear drag coefficient [m/s]
amrex::Real rdrag = amrex::Real(3e-4);
// Quadratic drag coefficient [dimensionless]
amrex::Real rdrag2 = amrex::Real(3e-3);

// Momentum stress scales
// Momentum stress scales [m]
amrex::Real Zob = amrex::Real(2e-2);
amrex::Real Zos = amrex::Real(2e-2);

amrex::Real Cdb_max = amrex::Real(0.5);
amrex::Real Cdb_min = amrex::Real(1e-6);

// Linear equation of state parameters
amrex::Real R0 = amrex::Real(1028); // background density value (Kg/m3) used in Linear Equation of State
amrex::Real S0 = amrex::Real(35.0); // background salinity (nondimensional) constant
Expand Down
5 changes: 4 additions & 1 deletion Source/TimeIntegration/REMORA_gls.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -343,7 +343,6 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke,
}

Real Zos_min = std::max(solverChoice.Zos, 0.0001_rt);
Real Zob_min = std::max(solverChoice.Zob, 0.0001_rt);
Real Zos_eff = Zos_min;
Real Gadv = 1.0_rt/3.0_rt;
Real eps = 1.0e-10_rt;
Expand Down Expand Up @@ -484,6 +483,8 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke,
Array4<Real> const& msku = vec_msku[lev]->array(mfi);
Array4<Real> const& mskv = vec_mskv[lev]->array(mfi);

Array4<Real> const& ZoBot = vec_ZoBot[lev]->array(mfi);

FArrayBox fab_tmp_buoy(bx_growloxy,1, amrex::The_Async_Arena()); fab_tmp_buoy.template setVal<RunOn::Device>(0.);
FArrayBox fab_tmp_shear(bx_growloxy,1, amrex::The_Async_Arena()); fab_tmp_shear.template setVal<RunOn::Device>(0.);

Expand Down Expand Up @@ -720,6 +721,7 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke,
// Compute production and dissipation terms.
ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
Real Zob_min = std::max(ZoBot(i,j,0), 0.0001_rt);
//----------------------------------------------------------------------
// Time-step dissipation and vertical diffusion terms implicitly.
//----------------------------------------------------------------------
Expand Down Expand Up @@ -885,6 +887,7 @@ REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke,
});
ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
Real Zob_min = std::max(ZoBot(i,j,0), 0.0001_rt);
Akv(i,j,N+1)=Akv_bak+L_sft*Zos_eff*gls_cmu0*
std::sqrt(tke(i,j,N+1,nnew));
Akv(i,j,0)=Akv_bak+vonKar*Zob_min*gls_cmu0*
Expand Down
67 changes: 58 additions & 9 deletions Source/TimeIntegration/REMORA_setup_step.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ REMORA::setup_step (int lev, Real time, Real dt_lev)
MultiFab mf_AK(ba,dm,1,IntVect(NGROW,NGROW,0)); //2d missing j coordinate
MultiFab mf_DC(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
MultiFab mf_Hzk(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
MultiFab mf_logdrg_tmp(ba,dm,1,IntVect(NGROW,NGROW,0));

MultiFab* mf_z_r = vec_z_r[lev].get();
MultiFab* mf_z_w = vec_z_w[lev].get();
Expand All @@ -71,6 +72,8 @@ REMORA::setup_step (int lev, Real time, Real dt_lev)
std::unique_ptr<MultiFab>& mf_sustr = vec_sustr[lev];
std::unique_ptr<MultiFab>& mf_svstr = vec_svstr[lev];
std::unique_ptr<MultiFab>& mf_rdrag = vec_rdrag[lev];
std::unique_ptr<MultiFab>& mf_rdrag2 = vec_rdrag2[lev];
std::unique_ptr<MultiFab>& mf_ZoBot = vec_ZoBot[lev];
std::unique_ptr<MultiFab>& mf_bustr = vec_bustr[lev];
std::unique_ptr<MultiFab>& mf_bvstr = vec_bvstr[lev];

Expand Down Expand Up @@ -112,29 +115,75 @@ REMORA::setup_step (int lev, Real time, Real dt_lev)

auto N = Geom(lev).Domain().size()[2]-1; // Number of vertical "levs" aka, NZ

const Real vonKar = solverChoice.vonKar;
const Real Cdb_min = solverChoice.Cdb_min;
const Real Cdb_max = solverChoice.Cdb_max;

for ( MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
Array4<Real const> const& rdrag = mf_rdrag->const_array(mfi);
Array4<Real > const& bustr = mf_bustr->array(mfi);
Array4<Real > const& bvstr = mf_bvstr->array(mfi);
Array4<Real const> const& uold = U_old.const_array(mfi);
Array4<Real const> const& vold = V_old.const_array(mfi);
Array4<Real > const& logdrg_tmp = mf_logdrg_tmp.array(mfi);
Array4<Real const> const& z_r = mf_z_r->const_array(mfi);
Array4<Real const> const& z_w = mf_z_w->const_array(mfi);

Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
Box gbx2D = gbx2;
gbx2D.makeSlab(2,0);
Box ubx1 = mfi.grownnodaltilebox(0,IntVect(NGROW-1,NGROW-1,0));
Box ubx1D = ubx1;
ubx1D.makeSlab(2,0);
Box vbx1 = mfi.grownnodaltilebox(1,IntVect(NGROW-1,NGROW-1,0));
Box vbx1D = vbx1;
vbx1D.makeSlab(2,0);
// Set bottom stress as defined in set_vbx.F
ParallelFor(ubx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
bustr(i,j,0) = 0.5_rt * (rdrag(i-1,j,0)+rdrag(i,j,0))*(uold(i,j,0));
});
ParallelFor(vbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
bvstr(i,j,0) = 0.5_rt * (rdrag(i,j-1,0)+rdrag(i,j,0))*(vold(i,j,0));
});
if (solverChoice.bottom_stress_type == BottomStressType::linear) {
Array4<Real const> const& rdrag = mf_rdrag->const_array(mfi);
ParallelFor(ubx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
bustr(i,j,0) = 0.5_rt * (rdrag(i-1,j,0)+rdrag(i,j,0))*(uold(i,j,0));
});
ParallelFor(vbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
bvstr(i,j,0) = 0.5_rt * (rdrag(i,j-1,0)+rdrag(i,j,0))*(vold(i,j,0));
});
} else if (solverChoice.bottom_stress_type == BottomStressType::quadratic) {
Array4<Real const> const& rdrag2 = mf_rdrag2->const_array(mfi);
ParallelFor(ubx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
Real avg_v = 0.25_rt * (vold(i,j,0) + vold(i,j+1,0) + vold(i-1,j,0) + vold(i-1,j+1,0));
Real vel_mag = std::sqrt(uold(i,j,0)*uold(i,j,0) + avg_v * avg_v);
bustr(i,j,0) = 0.5_rt * (rdrag2(i-1,j,0) + rdrag2(i,j,0)) * uold(i,j,0) * vel_mag;
});
ParallelFor(vbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
Real avg_u = 0.25_rt * (uold(i,j,0) + uold(i+1,j,0) + uold(i,j-1,0) + uold(i+1,j-1,0));
Real vel_mag = std::sqrt(avg_u * avg_u + vold(i,j,0) * vold(i,j,0));
bvstr(i,j,0) = 0.5_rt * (rdrag2(i,j-1,0) + rdrag2(i,j,0)) * vold(i,j,0) * vel_mag;
});
} else if (solverChoice.bottom_stress_type == BottomStressType::logarithmic) {
Array4<Real const> const& ZoBot = mf_ZoBot->const_array(mfi);
ParallelFor(gbx2D, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
Real logz = 1.0_rt / (std::log((z_r(i,j,0) - z_w(i,j,0)) / ZoBot(i,j,0)));
Real cff = vonKar * vonKar * logz * logz;
logdrg_tmp(i,j,0) = std::min(Cdb_max,std::max(Cdb_min,cff));
});
ParallelFor(ubx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
Real avg_v = 0.25_rt * (vold(i,j,0) + vold(i,j+1,0) + vold(i-1,j,0) + vold(i-1,j+1,0));
Real vel_mag = std::sqrt(uold(i,j,0)*uold(i,j,0) + avg_v * avg_v);
bustr(i,j,0) = 0.5_rt * (logdrg_tmp(i-1,j,0)+logdrg_tmp(i,j,0)) * uold(i,j,0) * vel_mag;
});
ParallelFor(vbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
Real avg_u = 0.25_rt * (uold(i,j,0) + uold(i+1,j,0) + uold(i,j-1,0) + uold(i+1,j-1,0));
Real vel_mag = std::sqrt(avg_u * avg_u + vold(i,j,0) * vold(i,j,0));
bvstr(i,j,0) = 0.5_rt * (logdrg_tmp(i,j-1,0) + logdrg_tmp(i,j,0)) * vold(i,j,0) * vel_mag;
});
}
}
FillPatch(lev, time, *vec_bustr[lev].get(), GetVecOfPtrs(vec_bustr), BCVars::u2d_simple_bc, BdyVars::null,0,true,false);
FillPatch(lev, time, *vec_bvstr[lev].get(), GetVecOfPtrs(vec_bvstr), BCVars::v2d_simple_bc, BdyVars::null,0,true,false);
Expand Down
4 changes: 4 additions & 0 deletions Tests/CTestList.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,8 @@ if(WIN32)
add_test_r(Upwelling "Upwelling/*/upwelling.exe" "plt00010")
add_test_r(Upwelling_GLS "Upwelling/*/upwelling.exe" "plt00010")
add_test_r(Upwelling_NLEOS "Upwelling/*/upwelling.exe" "plt00010")
add_test_r(Upwelling_qdrag "Upwelling/*/upwelling.exe" "plt00010")
add_test_r(Upwelling_logdrag "Upwelling/*/upwelling.exe" "plt00010")
add_test_r(Channel_Test "Channel_Test/*/channel_test.exe" "plt00010")
add_test_r(DoubleGyre "DoubleGyre/*/doublegyre.exe" "plt00010")
else()
Expand All @@ -112,6 +114,8 @@ else()
add_test_r(Upwelling "Upwelling/upwelling" "plt00010")
add_test_r(Upwelling_GLS "Upwelling/upwelling" "plt00010")
add_test_r(Upwelling_NLEOS "Upwelling/upwelling" "plt00010")
add_test_r(Upwelling_qdrag "Upwelling/upwelling" "plt00010")
add_test_r(Upwelling_logdrag "Upwelling/upwelling" "plt00010")
add_test_r(Channel_Test "Channel_Test/channel_test" "plt00010")
add_test_r(DoubleGyre "DoubleGyre/doublegyre" "plt00010")
endif()
Expand Down
30 changes: 30 additions & 0 deletions Tests/REMORA_Gold_Files/Upwelling_logdrag/Header
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
HyperCLaw-V1.1
5
temp
salt
x_velocity
y_velocity
z_velocity
3
3000
0
0 0 -150
41000 80000 0

((0,0,0) (40,79,15) (0,0,0))
10
1000 1000 9.375
0
0
0 1 3000
10
0 41000
0 80000
-150 0
Level_0/Cell
1
3
amrexvec_nu_x
amrexvec_nu_y
amrexvec_nu_z
Level_0/Nu_nd
Binary file added Tests/REMORA_Gold_Files/Upwelling_logdrag/Level_0/Cell_D_00000
View file
Binary file not shown.
16 changes: 16 additions & 0 deletions Tests/REMORA_Gold_Files/Upwelling_logdrag/Level_0/Cell_H
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
1
1
5
0
(1 0
((0,0,0) (40,79,15) (0,0,0))
)
1
FabOnDisk: Cell_D_00000 0

1,5
1.45204851068578318e+01,3.49999999999999361e+01,-8.62178518000840041e-04,-6.34471600115517966e-05,0.00000000000000000e+00,

1,5
2.19344849766820396e+01,3.50000000000000568e+01,6.29280714783082826e-06,5.61017252629464303e-05,0.00000000000000000e+00,

Binary file added Tests/REMORA_Gold_Files/Upwelling_logdrag/Level_0/Nu_nd_D_00000
View file
Binary file not shown.
16 changes: 16 additions & 0 deletions Tests/REMORA_Gold_Files/Upwelling_logdrag/Level_0/Nu_nd_H
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
1
1
3
0
(1 0
((0,0,0) (41,80,16) (1,1,1))
)
1
FabOnDisk: Nu_nd_D_00000 0

1,3
0.00000000000000000e+00,0.00000000000000000e+00,-1.25844679229360145e-05,

1,3
0.00000000000000000e+00,0.00000000000000000e+00,1.26471716899404043e+02,

Loading

0 comments on commit 1b99d10

Please sign in to comment.