last PR forgot clang-format, python code cleanup

src/pks/flow/constitutive_relations/wrm/rel_perm_frzBC_evaluator.cc

@@ -158,7 +158,7 @@ RelPermFrzBCEvaluator::Evaluate_(const State& S, const std::vector<CompositeVect
     int index = (*wrms_->first)[c];
     double sat_res = wrms_->second[index]->residualSaturation();
     double coef = BrooksCoreyFrzCoef::frzcoef(sat_c[0][c], sat_gas_c[0][c], omega_);
-    res_c[0][c] = std::max(wrms_->second[index]->k_relative(1. - sat_gas_c[0][c]) * coef,  min_val_);
+    res_c[0][c] = std::max(wrms_->second[index]->k_relative(1. - sat_gas_c[0][c]) * coef, min_val_);
   }
 
   // -- Potentially evaluate the model on boundary faces as well.
@@ -187,7 +187,8 @@ RelPermFrzBCEvaluator::Evaluate_(const State& S, const std::vector<CompositeVect
       double krel;
 
       double coef_b = BrooksCoreyFrzCoef::frzcoef(sat_bf[0][bf], sat_gas_bf[0][bf], omega_);
-      double coef_c =BrooksCoreyFrzCoef::frzcoef(sat_c[0][cells[0]], sat_gas_c[0][cells[0]], omega_);
+      double coef_c =
+        BrooksCoreyFrzCoef::frzcoef(sat_c[0][cells[0]], sat_gas_c[0][cells[0]], omega_);
 
       if (boundary_krel_ == BoundaryRelPerm::HARMONIC_MEAN) {
         double krelb =

src/pks/flow/constitutive_relations/wrm/rel_perm_frzBC_evaluator.hh

@@ -124,4 +124,3 @@ class RelPermFrzBCEvaluator : public EvaluatorSecondaryMonotypeCV {
 
 } // namespace Flow
 } // namespace Amanzi
-

src/pks/flow/constitutive_relations/wrm/wrm_brooks_corey.cc

@@ -27,8 +27,8 @@ WRMBrooksCorey::WRMBrooksCorey(Teuchos::ParameterList& plist) : plist_(plist)
 };
 
 
-void 
-WRMBrooksCorey::InitializeFromPlist_() 
+void
+WRMBrooksCorey::InitializeFromPlist_()
 {
   lambda_ = plist_.get<double>("Brooks Corey lambda [-]");
   b_ = 1. / lambda_; // clapp hornberger b
@@ -45,7 +45,7 @@ WRMBrooksCorey::InitializeFromPlist_()
  * The original curve is regulized on interval (s0, 1) using the
  * Hermite interpolant of order 3. Formulas (3.11)-(3.12).
  ****************************************************************** */
-double 
+double
 WRMBrooksCorey::k_relative(double s)
 {
   if (s <= s0_) {
@@ -62,7 +62,7 @@ WRMBrooksCorey::k_relative(double s)
 /* ******************************************************************
  * D Relative permeability / D saturation
  ****************************************************************** */
-double 
+double
 WRMBrooksCorey::d_k_relative(double s)
 {
   if (s <= s0_) {
@@ -71,7 +71,7 @@ WRMBrooksCorey::d_k_relative(double s)
     return dkdse / (1. - sr_);
   } else if (s == 1.) {
     return 0.0;
-  } else{
+  } else {
     return fit_kr_.Derivative(s);
   }
 }
@@ -80,7 +80,7 @@ WRMBrooksCorey::d_k_relative(double s)
 /* ******************************************************************
  * Saturation formula (3.5)-(3.8).
  ****************************************************************** */
-double 
+double
 WRMBrooksCorey::saturation(double pc)
 {
   if (pc <= p_sat_) {
@@ -94,10 +94,10 @@ WRMBrooksCorey::saturation(double pc)
 /* ******************************************************************
  * Derivative of the saturation formula w.r.t. capillary pressure.
  ****************************************************************** */
-double 
+double
 WRMBrooksCorey::d_saturation(double pc)
 {
-  if ( pc <= p_sat_) {
+  if (pc <= p_sat_) {
     return 0.;
   } else {
     return -(1. - sr_) * lambda_ * pow(p_sat_ / pc, lambda_) / pc;
@@ -108,7 +108,7 @@ WRMBrooksCorey::d_saturation(double pc)
 /* ******************************************************************
  * Pressure as a function of saturation, formula (3.9).
  ****************************************************************** */
-double 
+double
 WRMBrooksCorey::capillaryPressure(double s)
 {
   double se = (s - sr_) / (1. - sr_);
@@ -121,7 +121,7 @@ WRMBrooksCorey::capillaryPressure(double s)
 /* ******************************************************************
  * Derivative of pressure formulat w.r.t. saturation.
  ****************************************************************** */
-double 
+double
 WRMBrooksCorey::d_capillaryPressure(double s)
 {
   double se = (s - sr_) / (1. - sr_);

src/pks/flow/constitutive_relations/wrm/wrm_brooks_corey.hh

@@ -73,18 +73,17 @@ class WRMBrooksCorey : public WRM {
 
   Teuchos::ParameterList& plist_;
 
-  double lambda_, p_sat_;  // Brooks and Corey parameters: lambda, p_sat
-  double sr_;  // residual saturation
-  double b_;  // frequently used constant, clapp-hornberger b = 1/lambda
+  double lambda_, p_sat_; // Brooks and Corey parameters: lambda, p_sat
+  double sr_;             // residual saturation
+  double b_;              // frequently used constant, clapp-hornberger b = 1/lambda
 
   double s0_; // regularization threshold in saturation
   Amanzi::Utils::Spline fit_kr_;
 
-  static Utils::RegisteredFactory<WRM,WRMBrooksCorey> factory_;
-
+  static Utils::RegisteredFactory<WRM, WRMBrooksCorey> factory_;
 };
 
-}  // namespace
-}  // namespace
+} // namespace Flow
+} // namespace Amanzi
 
 #endif

src/pks/flow/constitutive_relations/wrm/wrm_brooks_corey_reg.hh

@@ -13,7 +13,7 @@
 namespace Amanzi {
 namespace Flow {
 
-Utils::RegisteredFactory<WRM,WRMBrooksCorey> WRMBrooksCorey::factory_("Brooks-Corey");
+Utils::RegisteredFactory<WRM, WRMBrooksCorey> WRMBrooksCorey::factory_("Brooks-Corey");
 
 } // namespace Flow
 } // namespace Amanzi

src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_defs.hh

@@ -60,10 +60,10 @@ struct ModelParams {
       evap_transition_width(100.),   // transition on evaporation from surface to
                                      // evaporation from subsurface [m],
                                      // THIS IS DEPRECATED
-      KB(0.),                        // a paramter denoting the log ratio of momentum roughness length
-                                     // to vapor roughness length for vapor flux or 
-                                     // to thermal roughness length for sensible heat flux.
-      // Da0_a, Da0_b, Cd0_c, Cd0_d are fitting parameters in the formulization of the log ratio 
+      KB(0.), // a paramter denoting the log ratio of momentum roughness length
+              // to vapor roughness length for vapor flux or
+              // to thermal roughness length for sensible heat flux.
+      // Da0_a, Da0_b, Cd0_c, Cd0_d are fitting parameters in the formulization of the log ratio
       // of momentum roughness length to the vapor roughness length for vapor flux [-]. Setting
       // Da0_a = Cd0_c = Cd0_d = 0 means momentum roughness length = vapor roughness length.
       Da0_a(0.),

src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.cc

@@ -292,13 +292,13 @@ UpdateEnergyBalanceWithSnow_Inner(const GroundProperties& surf,
   double vapor_pressure_skin = SaturatedVaporPressure(snow.temp);
   // KB here is to formulize the log ration between momentum roughness length and vapor roughness length
   // to add the effect of surface microtopography. Details see Gao et al. 2021 (WRR) Eq.(17) and similar
-  // studies by Mölder & Lindroth 2001 (Agricultural and Forest Meteorology) but in heat transport. 
-  // Theories see Chapter 4 of book [Brutsaert, W. (1982). Evaporation into the atmosphere: Theory, 
-  // history, and applications]. Da0_a, Da0_b, Cd0_c, Cd0_d here are four fitting coefficients in the 
-  // formulization. There are several different values proposed by different studies (see Brutsaert 1982). 
-  // In Gao et al. (2021), these four coefficients are determined by lab experiments. 
-  // If set Da0_a, Cd0_c, Cd0_d to 0, it means that vapor roughness length is assumed equal to momentum 
-  // roughness length. Currently, this approach still needs testing or other formulization approaches 
+  // studies by Mölder & Lindroth 2001 (Agricultural and Forest Meteorology) but in heat transport.
+  // Theories see Chapter 4 of book [Brutsaert, W. (1982). Evaporation into the atmosphere: Theory,
+  // history, and applications]. Da0_a, Da0_b, Cd0_c, Cd0_d here are four fitting coefficients in the
+  // formulization. There are several different values proposed by different studies (see Brutsaert 1982).
+  // In Gao et al. (2021), these four coefficients are determined by lab experiments.
+  // If set Da0_a, Cd0_c, Cd0_d to 0, it means that vapor roughness length is assumed equal to momentum
+  // roughness length. Currently, this approach still needs testing or other formulization approaches
   // may also be added later for testing. So keep Da0_a = Cd0_c = Cd0_d = 0 for users.
   double u_star =
     met.Us * c_von_Karman /
@@ -307,7 +307,8 @@ UpdateEnergyBalanceWithSnow_Inner(const GroundProperties& surf,
                CalcRoughnessFactor(snow.height, surf.roughness, snow.roughness) /
                params.dynamic_viscosity_air;
   double KB =
-    (params.Da0_a * std::pow(Re0, params.Da0_b) - (params.Cd0_c * std::log(Re0) + params.Cd0_d)) * c_von_Karman;
+    (params.Da0_a * std::pow(Re0, params.Da0_b) - (params.Cd0_c * std::log(Re0) + params.Cd0_d)) *
+    c_von_Karman;
   double Dhe_latent = WindFactor(
     met.Us, met.Z_Us, CalcRoughnessFactor(snow.height, surf.roughness, snow.roughness), KB);
   eb.fQe = LatentHeat(Dhe_latent * Sqig,
@@ -386,7 +387,8 @@ UpdateEnergyBalanceWithoutSnow(const GroundProperties& surf,
   double u_star = met.Us * c_von_Karman / std::log(met.Z_Us / surf.roughness);
   double Re0 = params.density_air * u_star * surf.roughness / params.dynamic_viscosity_air;
   double KB =
-    (params.Da0_a * std::pow(Re0, params.Da0_b) - (params.Cd0_c * std::log(Re0) + params.Cd0_d)) * c_von_Karman;
+    (params.Da0_a * std::pow(Re0, params.Da0_b) - (params.Cd0_c * std::log(Re0) + params.Cd0_d)) *
+    c_von_Karman;
   double Dhe_latent = WindFactor(met.Us, met.Z_Us, surf.roughness, KB);
   double Rsoil = EvaporativeResistanceGround(surf, met, params, vapor_pressure_skin);
   double coef = 1.0 / (Rsoil + 1.0 / (Dhe_latent * Sqig));

tools/utils/convert_parameters_vg2bc.py

@@ -1,3 +1,8 @@
+"""Computes Brooks-Corey water retention curve parameters given a set
+of van Genuchten parameters, using the method of Lenhard et al. (1989)
+or Ma et al. (1999) method 2.
+"""
+
 from argparse import ArgumentParser
 from plot_wrm import *
 import numpy as np
@@ -37,9 +42,10 @@ def plot_vg_and_bc(alpha, n, sr, smoothing_interval_sat=0.05):
     ax[0].plot(sat, kr_bc)
     ax[0].set_xlabel('Saturation')
     ax[0].set_ylabel('Relative permeability')
-    ax[1].plot(sat, pc_vg, label='van Genuchten')
-    ax[1].plot(sat, pc_bc, label='Brooks-Corey, \nlambda={:.2f}, \
-               \nclapp_horn_b={:.2f}, \np_sat={:.2f} Pa'.format(1/clapp_horn_b, clapp_horn_b, p_sat))
+    ax[1].plot(sat, pc_vg, label='van Genuchten, \n  alpha={:.2g}, \
+               \n  n={:.2f}'.format(alpha, n))
+    ax[1].plot(sat, pc_bc, label='Brooks-Corey, \n  lambda={:.2f}, \
+               \n  clapp_horn_b={:.2f}, \n  p_sat={:.2f} Pa'.format(1/clapp_horn_b, clapp_horn_b, p_sat))
 
     ax[1].set_yscale('log')
     ax[1].set_xlabel('Saturation')
@@ -48,11 +54,11 @@ def plot_vg_and_bc(alpha, n, sr, smoothing_interval_sat=0.05):
 
 
 def get_args():
-    parser = ArgumentParser()
+    parser = ArgumentParser(description=__doc__)
     parser.add_argument('alpha', type=float, help='van Genuchten alpha [Pa^-1]')
     parser.add_argument('n', type=float, help='van Genuchten n')
     parser.add_argument('--sr', type=float, help='residual saturation')
-    parser.add_argument('--smooth_sat', type=float, help='smoothing_interval_sat')
+    parser.add_argument('--smooth_sat', type=float, default=0, help='smoothing_interval_sat')
     parser.add_argument('--plot', action='store_true', help='plot van Genuchten vs. Brooks-Corey')
     args = parser.parse_args()
     return args
@@ -67,6 +73,3 @@ def get_args():
     if args.plot:
         plot_vg_and_bc(args.alpha, args.n, args.sr, args.smooth_sat)
         plt.show()
-
-
-