Munk: test modified for partial slip condition

SlimFunction: sediment erosion functions added

modules/shallowWater/tests/munk/test.py

@@ -82,7 +82,7 @@ bathGradient=functionConstant([0,0,0])
 claw.setBathymetryGradient(bathGradient)
 claw.setBathymetry(bath)
 #We enforce the boundary conditions (True means that slip is allowed along the boundaries)
-claw.addBoundaryCondition('Wall',claw.newBoundaryWall(True))
+claw.addBoundaryCondition('Wall',claw.newBoundaryWall(1.0))
 nu=functionConstant([10000])
 claw.setDiffusivity(nu)
 #The linear system and its solver

modules/slimFunction/slimFunction.cpp

@@ -1391,4 +1391,335 @@ void sedimentBottomFlux2D::operator()(functorCache &cm, fullMatrix<double> &val)
     }
 }
 
+sedimentBottomFlux2D2::sedimentBottomFlux2D2(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity, double dt, double d50, double d90, double rho_sed, double shields, int evaluation_point) {
+  _sed = sed;
+  _sol = sol;
+  _wind = wind;
+  _bath = bath;
+  _sedGround = sedGround;
+  _gravity = gravity;
+  _dt = dt;
+  _d50 = d50;
+  _d90 = d90;
+  _rho_sed = rho_sed;
+  _shields = shields;
+  _evaluation_point = evaluation_point;
+}
+
+void sedimentBottomFlux2D2::operator()(functorCache &cm, fullMatrix<double> &val) const {
+    const fullMatrix<double> &sed = cm.get(*_sed);
+    const fullMatrix<double> &sol = cm.get(*_sol);
+    const fullMatrix<double> &wind = cm.get(*_wind);
+    const fullMatrix<double> &bath = cm.get(*_bath);
+    const fullMatrix<double> &sedGround = cm.get(*_sedGround);
+    const fullMatrix<double> &gravity = cm.get(*_gravity);
+    val.resize(cm.nEvaluationPoint(),1,true);
+    double g = 9.81;
+    double K = 0.41;
+    double nu = 1e-6;
+    double s = _rho_sed/1000;
+    double print = _evaluation_point;
+
+    double negSed = 0;
+    double Sp = 0.7;
+    double M = 53.5*exp(-0.65*Sp);
+    double N = 5.65*exp(-2.5*Sp);
+    double n = 0.7+0.9*Sp;
+    double d_ast = _d50*pow((s-1.0)*g/pow(nu,2),(1.0/3.0));
+    printf("d_ast = %.4g \n", d_ast);
+    double ws = M*nu/(N*_d50)*pow((sqrt(1.0/4.0+pow((4.0*N/(3.0*pow(M,2.0))*pow(d_ast,3.0)),1.0)/n)-1.0/2.0),n);
+    printf("ws = %.4g \n", ws);
+    double k_s = 3.0*_d50;
+    printf("k_s = %.4g \n", k_s);
+    double u_astcr = _shields;
+	
+    for (int i = 0; i < cm.nEvaluationPoint(); i++) {
+      negSed = 0;
+      double uBot = hypot(sol(i,1), sol(i,2));
+      double H = bath(i,0)+sol(i,0);
+
+      double u_ast0 = ws*0.4;
+      double u_ast = uBot/(1.0/K*log(u_ast0*H/nu)+3.0);
+      int iter_u = 0;
+      while (abs(u_ast-u_ast0)>1e-5 && iter_u<30) {
+        u_ast0 = u_ast;
+        u_ast = uBot/(1.0/K*log(u_ast0*H/nu)+3.0);
+        iter_u++;
+      }
+      if (iter_u == 30) Msg::Fatal("Could not find the bottom velocity. Please refer to sedimentBottomFlux2D function");
+
+      double Up_ast = (uBot*pow(g,0.5))/(18.0*log(2.0*H/_d90));
+      double T = (Up_ast>u_astcr) ? (pow((Up_ast/u_astcr),2)-1.0) : 0;
+      double delta = std::min(1.0/2.0*k_s,0.01*H);
+      double c_ast_b = 0.015*(_d50*pow(T,1.5))/(delta*pow(d_ast,0.3));
+
+
+      double z_ast = ws/(K*u_ast);
+      double alpha = 3.25 + 0.55*log(z_ast);
+
+      
+      
+      if (i==print) printf("uBot = %.4g \n", uBot);
+      if (i==print) printf("H = %.4g \n", H);
+      if (i==print) printf("u_ast = %.4g \n", u_ast);
+      if (i==print) printf("T = %.4g \n", T);
+      if (i==print) printf("delta = %.4g \n", delta);
+      if (i==print) printf("c_ast_b = %.4g \n", c_ast_b);    
+      if (i==print) printf("z_ast = %.4g \n", z_ast);
+      if (i==print) printf("alpha = %.4g \n", alpha);
+      if (i==print) printf("Css = %.4g \n", sed(i,0));
+      
+
+      double settling = alpha*ws*sed(i,0);
+      if (i==print) printf("settling = %.4g \n", settling);
+      double erosionTotal = ws*c_ast_b;
+      if (i==print) printf("erosionTotal = %.4g \n", erosionTotal);
+      if (erosionTotal < 0) Msg::Fatal("Total erosion should always be positive! dt should be used for the max erosion term");
+      if (sed(i,0) < -1e-9) {
+          negSed++;
+          settling = 0;
+      }
+      val(i,0) = (erosionTotal - settling)/H;
+      val(i,0) = std::max(-sed(i,0)/_dt, val(i,0));		//No more settling than available
+      val(i,0) = std::min(sedGround(i,0)/(_dt*H), val(i,0)); 	//No more erosion than available
+      
+    }
+    if (negSed > 0) Msg::Warning("Wrong!! Negative sediment concentration");
+}
+
+sedimentBottomFlux2D3::sedimentBottomFlux2D3(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity, double dt, double d50, double S0, double rho_sed, double shields, int evaluation_point) {
+  _sed = sed;
+  _sol = sol;
+  _wind = wind;
+  _bath = bath;
+  _sedGround = sedGround;
+  _gravity = gravity;
+  _dt = dt;
+  _d50 = d50;
+  _S0 = S0;
+  _rho_sed = rho_sed;
+  _shields = shields;
+  _evaluation_point = evaluation_point;
+}
+
+void sedimentBottomFlux2D3::operator()(functorCache &cm, fullMatrix<double> &val) const {
+    const fullMatrix<double> &sed = cm.get(*_sed);
+    const fullMatrix<double> &sol = cm.get(*_sol);
+    const fullMatrix<double> &wind = cm.get(*_wind);
+    const fullMatrix<double> &bath = cm.get(*_bath);
+    const fullMatrix<double> &sedGround = cm.get(*_sedGround);
+    const fullMatrix<double> &gravity = cm.get(*_gravity);
+    val.resize(cm.nEvaluationPoint(),1,true);
+    double negSed = 0;
+    double g = 9.81;
+    double K = 0.41;
+    double nu = 1e-6;
+    double s = _rho_sed/1000;
+    double print = _evaluation_point;
+
+    double d_ast = _d50*pow((s-1.0)*g/pow(nu,2),(1.0/3.0));
+    double ws = ((_rho_sed-1000.)*g*pow(_d50,2))/(18.*nu*1000.);
+    double u_astcr = (d_ast>10) ? ws*0.4 : ws*4.0/d_ast;
+
+    printf("d_ast = %.4f \n", d_ast);
+    printf("ws = %.4f \n", ws);
+	
+    for (int i = 0; i < cm.nEvaluationPoint(); i++) {
+      negSed = 0;
+      double U = hypot(sol(i,1), sol(i,2));
+      double H = bath(i,0)+sol(i,0);
+
+      double Cf = g*H*_S0/pow(sol(i,1),2.);
+      double u_ast = pow(Cf,(1./2.))*U;
+      double z_ast = ws/(K*u_ast);
+      double theta_cr = 0.1046; //pow(u_astcr,2.)/((s-1.)*g*_d50);
+      double theta_s = pow(u_ast,2.)/((s-1.)*g*_d50);
+      double alpha = 3.25 + 0.55*log(z_ast);
+
+      double ca = alpha*sed(i,0);
+      double ca_ast = (theta_s > theta_cr) ? (0.331*pow(theta_s-theta_cr,1.75))/(1.+0.72*pow(theta_s-theta_cr,1.75)) : 0;
+      
+      if (i==print) printf("U = %.4f \n", U);
+      if (i==print) printf("H = %.4f \n", H);
+      if (i==print) printf("Cf = %.4f \n", Cf);
+      if (i==print) printf("u_ast = %.4f \n", u_ast);
+      if (i==print) printf("z_ast = %.4f \n", z_ast);
+      if (i==print) printf("theta_s = %.4f \n", theta_s);
+      if (i==print) printf("theta_cr = %.4f \n", theta_cr);
+      if (i==print) printf("alpha = %.4f \n", alpha);
+      if (i==print) printf("ca = %.4f \n", ca);
+      if (i==print) printf("ca_ast = %.4f \n", ca_ast);
+      if (i==print) printf("Css = %.4f \n", sed(i,0));
+      
+
+      double settling = ws*ca;
+      if (i==print) printf("settling = %.4g \n", settling);
+      double erosionTotal = ws*ca_ast;
+      if (i==print) printf("erosionTotal = %.4g \n", erosionTotal);
+      if (erosionTotal < 0) Msg::Fatal("Total erosion should always be positive! dt should be used for the max erosion term");
+      if (sed(i,0) < -1e-9) {
+          negSed++;
+          settling = 0;
+      }
+      val(i,0) = (erosionTotal - settling)/H;
+      val(i,0) = std::max(-sed(i,0)/_dt, val(i,0));		//No more settling than available
+      val(i,0) = std::min(sedGround(i,0)/(_dt*H), val(i,0)); 	//No more erosion than available
+      
+    }
+    if (negSed > 0) Msg::Warning("Wrong!! Negative sediment concentration");
+}
+
+sedimentBottomFlux2D4::sedimentBottomFlux2D4(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity, double dt, double d50, double n, double rho_sed, double tau_ast_c, int evaluation_point) {
+  _sed = sed;
+  _sol = sol;
+  _wind = wind;
+  _bath = bath;
+  _sedGround = sedGround;
+  _gravity = gravity;
+  _dt = dt;
+  _d50 = d50;
+  _n = n;
+  _rho_sed = rho_sed;
+  _tau_ast_c = tau_ast_c;
+  _evaluation_point = evaluation_point;
+}
+
+void sedimentBottomFlux2D4::operator()(functorCache &cm, fullMatrix<double> &val) const {
+    const fullMatrix<double> &sed = cm.get(*_sed);
+    const fullMatrix<double> &sol = cm.get(*_sol);
+    const fullMatrix<double> &wind = cm.get(*_wind);
+    const fullMatrix<double> &bath = cm.get(*_bath);
+    const fullMatrix<double> &sedGround = cm.get(*_sedGround);
+    const fullMatrix<double> &gravity = cm.get(*_gravity);
+    val.resize(cm.nEvaluationPoint(),1,true);
+    double negSed = 0;
+    double g = 9.81;
+    double K = 0.41;
+    double nu = 1e-6;
+    double print = _evaluation_point;
+
+    double ws = ((_rho_sed-1000.)*g*pow(_d50,2))/(18.*nu*1000.);
+
+    printf("ws = %.4f \n", ws);
+	
+    for (int i = 0; i < cm.nEvaluationPoint(); i++) {
+      negSed = 0;
+      double U = hypot(sol(i,1), sol(i,2));
+      double H = bath(i,0)+sol(i,0);
+
+      double u_ast0 = ws*0.4;
+      double u_ast = U/(1.0/K*log(u_ast0*H/nu)+3.0);
+      int iter_u = 0;
+      while (abs(u_ast-u_ast0)>1e-5 && iter_u<30) {
+        u_ast0 = u_ast;
+        u_ast = U/(1.0/K*log(u_ast0*H/nu)+3.0);
+        iter_u++;
+      }
+      if (iter_u == 30) Msg::Fatal("Could not find the bottom velocity. Please refer to sedimentBottomFlux2D function");
+
+      double z_ast = ws/(K*u_ast);
+      double alpha = 3.25 + 0.55*log(z_ast);
+      double Re = u_ast*_d50/nu;
+
+      double tau_b = 1000.*g*pow(_n,2)*sol(i,1)*U/pow(H,(1./3.));
+      double tau_ast_s = tau_b/((_rho_sed - 1000.)*g*_d50);
+      double C_ast = (tau_ast_s > _tau_ast_c) ? 0.65*0.0024*(tau_ast_s/_tau_ast_c - 1)/(1+0.0024*(tau_ast_s/_tau_ast_c - 1)) : 0;
+      
+      if (i==print) printf("U = %.4f \n", U);
+      if (i==print) printf("H = %.4f \n", H);
+      if (i==print) printf("u_ast = %.4f \n", u_ast);
+      if (i==print) printf("z_ast = %.4f \n", z_ast);
+      if (i==print) printf("alpha = %.4f \n", alpha);
+      if (i==print) printf("Re = %.4f \n", Re);
+      if (i==print) printf("tau_b = %.4f \n", tau_b);
+      if (i==print) printf("tau_ast_s = %.4f \n", tau_ast_s);
+      if (i==print) printf("tau_ast_c = %.4f \n", _tau_ast_c);
+      if (i==print) printf("C_ast = %.4f \n", C_ast);
+      if (i==print) printf("Css = %.4f \n", sed(i,0));
+      
+
+      double settling = ws*sed(i,0);
+      if (i==print) printf("settling = %.4g \n", settling);
+      double erosionTotal = ws*alpha*C_ast;
+      if (i==print) printf("erosionTotal = %.4g \n", erosionTotal);
+      if (erosionTotal < 0) Msg::Fatal("Total erosion should always be positive! dt should be used for the max erosion term");
+      if (sed(i,0) < -1e-9) {
+          negSed++;
+          settling = 0;
+      }
+      val(i,0) = (erosionTotal - settling)/H;
+      val(i,0) = std::max(-sed(i,0)/_dt, val(i,0));		//No more settling than available
+      val(i,0) = std::min(sedGround(i,0)/(_dt*H), val(i,0)); 	//No more erosion than available
+      
+    }
+    if (negSed > 0) Msg::Warning("Wrong!! Negative sediment concentration");
+}
+
+sedimentBottomFlux2D5::sedimentBottomFlux2D5(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity, double dt, double n, double tau_c, double M, double k1, double b, int evaluation_point) {
+  _sed = sed;
+  _sol = sol;
+  _wind = wind;
+  _bath = bath;
+  _sedGround = sedGround;
+  _gravity = gravity;
+  _dt = dt;
+  _n = n;
+  _tau_c = tau_c;
+  _M = M;
+  _k1 = k1;
+  _b = b;
+  _evaluation_point = evaluation_point;
+}
+
+void sedimentBottomFlux2D5::operator()(functorCache &cm, fullMatrix<double> &val) const {
+    const fullMatrix<double> &sed = cm.get(*_sed);
+    const fullMatrix<double> &sol = cm.get(*_sol);
+    const fullMatrix<double> &wind = cm.get(*_wind);
+    const fullMatrix<double> &bath = cm.get(*_bath);
+    const fullMatrix<double> &sedGround = cm.get(*_sedGround);
+    const fullMatrix<double> &gravity = cm.get(*_gravity);
+    val.resize(cm.nEvaluationPoint(),1,true);
+    double negSed = 0;
+    double g = 9.81;
+    double print = _evaluation_point;
+
+    double m = 1.;
+	
+    for (int i = 0; i < cm.nEvaluationPoint(); i++) {
+      negSed = 0;
+      double U = hypot(sol(i,1), sol(i,2));
+      double H = bath(i,0)+sol(i,0);
+      double ws = _k1*pow(sed(i,0),_b);
+
+      double tau_b = 1000.*g*pow(_n,2.)*sol(i,1)*U/pow(H,(1./3.));
+      double tau_d = _tau_c;
+      double P1 = (tau_b<=tau_d) ? (1-tau_b/tau_d) : 0;
+
+      
+      if (i==print) printf("U = %.4f \n", U);
+      if (i==print) printf("H = %.4f \n", H);
+      if (i==print) printf("ws = %.4f \n", ws);
+      if (i==print) printf("tau_b = %.4f \n", tau_b);
+      if (i==print) printf("P1 = %.4f \n", P1);
+      if (i==print) printf("Css = %.4f \n", sed(i,0));
+      
+
+      double settling = P1*ws*sed(i,0);
+      if (i==print) printf("settling = %.4g \n", settling);
+      double erosionTotal = (tau_b>_tau_c) ? _M*pow((tau_b/_tau_c - 1),m) : 0;
+      if (i==print) printf("erosionTotal = %.4g \n", erosionTotal);
+      if (erosionTotal < 0) Msg::Fatal("Total erosion should always be positive! dt should be used for the max erosion term");
+      if (sed(i,0) < -1e-9) {
+          negSed++;
+          settling = 0;
+      }
+      val(i,0) = (erosionTotal - settling)/H;
+      val(i,0) = std::max(-sed(i,0)/_dt, val(i,0));		//No more settling than available
+      val(i,0) = std::min(sedGround(i,0)/(_dt*H), val(i,0)); 	//No more erosion than available
+      
+    }
+    if (negSed > 0) Msg::Warning("Wrong!! Negative sediment concentration");
+}
+
+
 

modules/slimFunction/slimFunction.h

@@ -394,6 +394,41 @@ class sedimentBottomFlux2D : public functor {
   sedimentBottomFlux2D(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity, double factorW, double wMax, double u0, double v0, double w0, double a01, double a02, double a1, double n, double omega2, double eros0Fact);
 };
 
+class sedimentBottomFlux2D2 : public functor {
+  const functor *_sed, *_sol, *_wind, *_bath, *_sedGround, *_gravity;
+  double _dt, _d50, _d90, _rho_sed, _shields;
+  int _evaluation_point;
+ public :
+  void operator()(functorCache &cm, fullMatrix<double> &val) const ;
+  sedimentBottomFlux2D2(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity, double dt, double d50, double d90, double rho_sed, double shields, int evaluation_point=1);
+};
+
+class sedimentBottomFlux2D3 : public functor {
+  const functor *_sed, *_sol, *_wind, *_bath, *_sedGround, *_gravity;
+  double _dt, _d50, _S0, _rho_sed, _shields;
+  int _evaluation_point;
+ public :
+  void operator()(functorCache &cm, fullMatrix<double> &val) const ;
+  sedimentBottomFlux2D3(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity, double dt, double d50, double S0, double rho_sed, double shields, int evaluation_point=1);
+};
+
+class sedimentBottomFlux2D4 : public functor {
+  const functor *_sed, *_sol, *_wind, *_bath, *_sedGround, *_gravity;
+  double _dt, _d50, _n, _rho_sed, _tau_ast_c;
+  int _evaluation_point;
+ public :
+  void operator()(functorCache &cm, fullMatrix<double> &val) const ;
+  sedimentBottomFlux2D4(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity, double dt, double d50, double n, double rho_sed, double tau_ast_c = 0.034, int evaluation_point = 1);
+};
+
+class sedimentBottomFlux2D5 : public functor {
+  const functor *_sed, *_sol, *_wind, *_bath, *_sedGround, *_gravity;
+  double _dt, _n, _tau_c, _M, _k1, _b;
+  int _evaluation_point;
+ public :
+  void operator()(functorCache &cm, fullMatrix<double> &val) const ;
+  sedimentBottomFlux2D5(const functor *sed, const functor *sol, const functor *wind, const functor *bath, const functor *sedGround, const functor *gravity, double dt, double n, double tau_c = 0.1, double M = 4e-5, double k1 = 0.01, double b = 0.5, int evaluation_point = 1);
+};
 
 #endif