From 9b09b5dcfbe4d1ce1b3bcee504cd651f188dc436 Mon Sep 17 00:00:00 2001
From: HenrikAsmuth <henrik.asmuth@geo.uu.se>
Date: Fri, 2 Sep 2022 16:28:42 +0200
Subject: [PATCH] Add geostrophic wind option to BoundaryLayer

---
 apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp | 19 ++++++++++++++-----
 1 file changed, 14 insertions(+), 5 deletions(-)

diff --git a/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp b/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
index 48af64b9f..c397092e0 100644
--- a/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
+++ b/apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
@@ -118,6 +118,8 @@ void multipleLevel(const std::string& configPath)
 
     const uint nodes_per_H = config.contains("nz")? config.getValue<uint>("nz"): 64;
 
+    const bool driveWithGeostrophicWind = config.contains("driveWithGeostrophicWind")? config.getValue<bool>("driveWithGeostrophicWind"): false;
+
     // all in s
     const float tStartOut   = config.getValue<real>("tStartOut");
     const float tOut        = config.getValue<real>("tOut");
@@ -155,7 +157,7 @@ void multipleLevel(const std::string& configPath)
 
     para->setPrintFiles(true);
 
-    para->setForcing(pressureGradientLB, 0, 0);
+    if(!driveWithGeostrophicWind) para->setForcing(pressureGradientLB, 0, 0);
     para->setVelocityLB(velocityLB);
     para->setViscosityLB(viscosityLB);
     para->setVelocityRatio( dx / dt );
@@ -200,10 +202,17 @@ void multipleLevel(const std::string& configPath)
     para->setHasWallModelMonitor(true);
     bcFactory.setStressBoundaryCondition(BoundaryConditionFactory::StressBC::StressPressureBounceBack);
 
-
-    // gridBuilder->setVelocityBoundaryCondition(SideType::PZ, 0.0, 0.0, 0.0);
-    gridBuilder->setSlipBoundaryCondition(SideType::PZ,  0.0,  0.0, 0.0);
-    bcFactory.setSlipBoundaryCondition(BoundaryConditionFactory::SlipBC::SlipBounceBack); //SlipCompressibleTurbulentViscosity
+    if(driveWithGeostrophicWind) 
+    {
+        real u_geostrophic = u_star/0.4*log(H/z0) * dt / dx ;
+        gridBuilder->setVelocityBoundaryCondition(SideType::PZ, u_geostrophic, 0.0, 0.0);
+        bcFactory.setVelocityBoundaryCondition(BoundaryConditionFactory::VelocityBC::VelocityCompressible);
+    }
+    else //drive with pressure gradient
+    {
+        gridBuilder->setSlipBoundaryCondition(SideType::PZ,  0.0,  0.0, 0.0);
+        bcFactory.setSlipBoundaryCondition(BoundaryConditionFactory::SlipBC::SlipBounceBack); 
+    }
 
     real cPi = 3.1415926535897932384626433832795;
     para->setInitialCondition([&](real coordX, real coordY, real coordZ, real &rho, real &vx, real &vy, real &vz) {
-- 
GitLab