C++ Soil::Cw2方法代码示例

double
SoilWater::MaxExfiltration (const Geometry& geo, const size_t edge,
                            const Soil& soil, const double T) const
{
  const size_t n = geo.edge_other (edge, Geometry::cell_above);
  const double h0 = h (n);
  const double K0 = soil.K (n, h0, h_ice (n), T);
  if (max_exfiltration_gradient > 0.0)
    return K0 * max_exfiltration_gradient;
  const double Cw2 = soil.Cw2 (n, h0);
  const double Theta0 = Theta (n);
  const double Theta_surf = soil.Theta_res (n);
  const double delta_Theta = Theta0 - Theta_surf;
  const double z0 = geo.cell_z (n);
  // Darcy formulated for Theta between middle of node and soil surface.
  return - (K0 / Cw2) * (delta_Theta / z0);
}

//.........这里部分代码省略.........
	  //Initialize diffusive matrix
	  Solver::Matrix diff (cell_size);
	  // diff = ublas::zero_matrix<double> (cell_size, cell_size);
	  diffusion (geo, Kedge, diff);

	  //Initialize gravitational matrix
	  ublas::vector<double> grav (cell_size); //ublass compatibility
	  grav = ublas::zero_vector<double> (cell_size);
	  gravitation (geo, Kedge, grav);

	  // Boundary matrices and vectors
	  ublas::banded_matrix<double>  Dm_mat (cell_size, cell_size, 
                                                0, 0); // Dir bc
	  Dm_mat = ublas::zero_matrix<double> (cell_size, cell_size);
	  ublas::vector<double>  Dm_vec (cell_size); // Dir bc
	  Dm_vec = ublas::zero_vector<double> (cell_size);
	  ublas::vector<double> Gm (cell_size); // Dir bc
	  Gm = ublas::zero_vector<double> (cell_size);
	  ublas::vector<double> B (cell_size); // Neu bc 
	  B = ublas::zero_vector<double> (cell_size);

	  lowerboundary (geo, groundwater, active_lysimeter, h,
                         Kedge,
                         dq, Dm_mat, Dm_vec, Gm, B, msg);
	  upperboundary (geo, soil, T, surface, state, remaining_water, h,
                         Kedge,
                         dq, Dm_mat, Dm_vec, Gm, B, ddt, debug, msg, dt);
          Darcy (geo, Kedge, h, dq); //for calculating drain fluxes 


	  //Initialize water capacity  matrix
	  ublas::banded_matrix<double> Cw (cell_size, cell_size, 0, 0);
	  for (size_t c = 0; c < cell_size; c++)
	    Cw (c, c) = soil.Cw2 (c, h[c]);
	  
          std::vector<double> h_std (cell_size);
          //ublas vector -> std vector 
          std::copy(h.begin (), h.end (), h_std.begin ());

#ifdef TEST_OM_DEN_ER_BRUGT
          for (size_t cell = 0; cell != cell_size ; ++cell) 
            {				
              S_macro (cell) = (S_matrix[cell] + S_drain[cell]) 
                * geo.cell_volume (cell);
            }
#endif

	  //Initialize sum matrix
	  Solver::Matrix summat (cell_size);  
	  noalias (summat) = diff + Dm_mat;

	  //Initialize sum vector
	  ublas::vector<double> sumvec (cell_size);  
	  sumvec = grav + B + Gm + Dm_vec - S_vol
#ifdef TEST_OM_DEN_ER_BRUGT
            - S_macro
#endif
            ; 

	  // QCw is shorthand for Qmatrix * Cw
	  Solver::Matrix Q_Cw (cell_size);
	  noalias (Q_Cw) = prod (Qmat, Cw);

	  //Initialize A-matrix
	  Solver::Matrix A (cell_size);  
	  noalias (A) = (1.0 / ddt) * Q_Cw - summat;  

void
SoilWater::tick_after (const Geometry& geo,
                       const Soil& soil, const SoilHeat& soil_heat, 
                       const bool initial, 
                       Treelog& msg)
{
  TREELOG_SUBMODEL (msg, "SoilWater");

  // We need old K for primary/secondary flux division.
  std::vector<double> K_old = K_cell_;

  // Update cells.
  const size_t cell_size = geo.cell_size ();
  daisy_assert (K_cell_.size () == cell_size);
  daisy_assert (Cw2_.size () == cell_size);
  daisy_assert (h_.size () == cell_size);
  daisy_assert (h_ice_.size () == cell_size);
  daisy_assert (K_old.size () == cell_size);
  daisy_assert (Theta_.size () == cell_size);
  daisy_assert (Theta_primary_.size () == cell_size);
  daisy_assert (Theta_secondary_.size () == cell_size);
  daisy_assert (Theta_tertiary_.size () == cell_size);

  double z_low = geo.top ();
  table_low = NAN;
  double z_high = geo.bottom ();
  table_high = NAN;

  for (size_t c = 0; c < cell_size; c++)
    {
      // Groundwater table.
      const double z = geo.cell_top (c);
      const double h = h_[c];
      const double table = z + h;
      if (h < 0)
        {
          if (approximate (z, z_low))
            {
              if (!std::isnormal (table_low)
                  || table < table_low)
                table_low = table;
            }
          else if (z < z_low)
            table_low = table;
        }
      else if (approximate (z, z_high))
        {
          if (!std::isnormal (table_high)
              || table > table_high)
            table_high = table;
        }
      else if (z > z_high)
        table_high = table;

      // Conductivity.
      K_cell_[c] = soil.K (c, h_[c], h_ice_[c], soil_heat.T (c));
      
      // Specific water capacity.
      Cw2_[c] = soil.Cw2 (c, h_[c]);

      // Primary and secondary water.
      if (Theta_[c] <= 0.0)
        {
          std::ostringstream tmp;
          tmp << "Theta[" << c << "] = " << Theta_[c];
          daisy_bug (tmp.str ());
          Theta_[c] = 1e-9;
        }
      const double h_lim = soil.h_secondary (c);
      if (h_lim >= 0.0 || h_[c] <= h_lim)
        {
          // No active secondary domain.
          Theta_primary_[c] = Theta_[c];
          Theta_secondary_[c] = 0.0;
        }
      else 
        {
          // Secondary domain activated.
          const double Theta_lim = soil.Theta (c, h_lim, h_ice_[c]);
          daisy_assert (Theta_lim > 0.0);
          if (Theta_[c] >= Theta_lim)
            {
              Theta_primary_[c] = Theta_lim;
              Theta_secondary_[c] = Theta_[c] - Theta_lim;
            }
          else
            { 
              std::ostringstream tmp;
              tmp << "h[" << c << "] = " << h_[c] 
                  << "; Theta[" << c << "] = " << Theta_[c] 
                  << "\nh_lim = " << h_lim << "; Theta_lim = " << Theta_lim
                  << "\nStrenge h > h_lim, yet Theta <= Theta_lim";
              msg.bug (tmp.str ());
              Theta_primary_[c] = Theta_[c];
              Theta_secondary_[c] = 0.0;
            }
        }
    }

  if (!std::isnormal (table_high))
//.........这里部分代码省略.........