C++ Treelog::debug方法代码示例

 static void propagate_debug (Treelog& msg, int nest, const std::string& text)
 {
   switch (nest)
     {
     case is_unknown:
     case is_debug:
     case is_plain:
     case is_warning:
     case is_error:
     case is_bug:
       msg.debug (text);
       break;
     case is_close:
       msg.close ();
       break;
     case is_touch:
       msg.touch ();
       break;
     case is_flush:
       msg.flush ();
       break;
     default:
       msg.open (text);
     }
 }

void 
Horizon::Implementation::initialize (Hydraulic& hydraulic,
                                     const Texture& texture, 
                                     const double quarts,
                                     int som_size,
                                     const bool top_soil,
                                     const double center_z,
                                     Treelog& msg)
{
  if (CEC < 0.0)
    {
      CEC = default_CEC (texture);
      std::ostringstream tmp;
      tmp << "(CEC " << CEC 
          << " [cmolc/kg]) ; Estimated from clay and humus.";
      msg.debug (tmp.str ());
    }

  hydraulic.initialize (texture, dry_bulk_density, top_soil, CEC, center_z, msg);
  if (som_size > 0)
    {
      // Fill out SOM_fractions and SOM_C_per_N.
      if (SOM_C_per_N.size () > 0 && SOM_C_per_N.size () < som_size + 0U)
        SOM_C_per_N.insert (SOM_C_per_N.end (),
                            som_size - SOM_C_per_N.size (), 
                            SOM_C_per_N.back ());
      if (SOM_fractions.size () > 0 && SOM_fractions.size () < som_size + 0U)
        SOM_fractions.insert (SOM_fractions.end (),
                              som_size - SOM_fractions.size (), 
                              0.0);
    }

  // Did we specify 'dry_bulk_density'?  Else calculate it now.
  if (dry_bulk_density < 0.0)
    {
      dry_bulk_density = texture.rho_soil_particles () 
        * (1.0 - hydraulic.porosity ());
      std::ostringstream tmp;
      tmp << "(dry_bulk_density " << dry_bulk_density 
          << " [g/cm^3]) ; Estimated from porosity and texture.";
      msg.debug (tmp.str ());
    }

  hor_heat.initialize (hydraulic, texture, quarts, msg);
  
}

  void initialize (const Texture& texture,
                   double rho_b, const bool top_soil, const double CEC,
                   const double center_z, Treelog& msg)
  {
    TREELOG_MODEL (msg);
    std::ostringstream tmp;

    // Find Theta_sat.
    if (Theta_sat < 0.0)
      {
        if (rho_b < 0.0)
          {
            msg.error ("You must specify either dry bulk density or porosity");
            rho_b = 1.5;
            tmp << "Forcing rho_b = "  << rho_b << " g/cm^3\n";
          }
        Theta_sat = 1.0 - rho_b / texture.rho_soil_particles ();
        tmp << "(Theta_sat " << Theta_sat << " [])\n";
        daisy_assert (Theta_sat < 1.0);
      }
    if (Theta_sat <= Theta_fc)
      {
        msg.error ("Field capacity must be below saturation point");
        Theta_sat = (1.0 + 4.0 * Theta_fc) / 5.0;
        tmp << "Forcing Theta_sat = " << Theta_sat << " []\n";
      }

    // Find Theta_wp.
    if (Theta_wp < 0.0)
      {
        const double clay_lim // USDA Clay
          = texture.fraction_of_minerals_smaller_than ( 2.0 /* [um] */);
        const double silt_lim // USDA Silt 
          = texture.fraction_of_minerals_smaller_than (50.0 /* [um] */);
        daisy_assert (clay_lim >= 0.0);
        daisy_assert (silt_lim >= clay_lim);
        daisy_assert (silt_lim <= 1.0);
        const double mineral = texture.mineral ();
        const double clay = mineral * clay_lim * 100 /* [%] */;
        const double silt = mineral * (silt_lim - clay_lim) * 100 /* [%] */;
        const double humus = texture.humus * 100 /* [%] */;
        // Madsen and Platou (1983).
        Theta_wp = 0.758 * humus + 0.520 * clay + 0.075 * silt + 0.42;
        Theta_wp /= 100.0;      // [%] -> []
      }

    b = find_b (Theta_wp, Theta_fc);
    h_b = find_h_b (Theta_wp, Theta_fc, Theta_sat, b);
    tmp << "(b " << b << " [])\n"
        << "(h_b " << h_b << " [cm])";
    msg.debug (tmp.str ());

    // Must be called last (K_init depends on the other parameters).
    Hydraulic::initialize (texture, rho_b, top_soil, CEC, center_z, msg);
  }    

void
SoilWater::drain (const std::vector<double>& v, Treelog& msg)
{
  forward_sink (v, v);

  daisy_assert (S_sum_.size () == v.size ());
  daisy_assert (S_drain_.size () == v.size ());
  daisy_assert (S_soil_drain_.size () == v.size ());
  for (unsigned i = 0; i < v.size (); i++)
    {
      if (v[i] < -1e-8)
        {
          std::ostringstream tmp;
          tmp << "draining " << v[i] << " [h^-1] from cell " << i;
          msg.debug (tmp.str ());
        }
      S_sum_[i] += v[i];
      S_drain_[i] += v[i];
      S_soil_drain_[i] += v[i];
    }
}

//.........这里部分代码省略.........
  Theta_error = Theta;

  // Start time loop 
  double time_left = dt;	// How much of the large time step left.
  double ddt = dt;		// We start with small == large time step.
  int number_of_time_step_reductions = 0;
  int iterations_with_this_time_step = 0;
  

  int n_small_time_steps = 0;
  
  while (time_left > 0.0)
    {
      if (ddt > time_left)
	ddt = time_left;

      std::ostringstream tmp_ddt;
      tmp_ddt << "Time t = " << (dt - time_left) 
              << "; ddt = " << ddt
              << "; steps " << n_small_time_steps 
              << "; time left = " << time_left;
      Treelog::Open nest (msg, tmp_ddt.str ());

      if (n_small_time_steps > 0
          && (n_small_time_steps%msg_number_of_small_time_steps) == 0)
        {
          msg.touch ();
          msg.flush ();
        }
      
      n_small_time_steps++;
      if (n_small_time_steps > max_number_of_small_time_steps) 
        {
          msg.debug ("Too many small timesteps");
          throw "Too many small timesteps";
        }
      
      // Initialization for each small time step.

      if (debug > 0)
	{
	  std::ostringstream tmp;
	  tmp << "h = " << h << "\n";
	  tmp << "Theta = " << Theta;
	  msg.message (tmp.str ());
	}

      int iterations_used = 0;
      h_previous = h;
      Theta_previous = Theta;

      if (debug == 5)
	{
	  std::ostringstream tmp;
	  tmp << "Remaining water at start: " << remaining_water;
	  msg.message (tmp.str ());
	}

      ublas::vector<double> h_conv;

      for (size_t cell = 0; cell != cell_size ; ++cell)
        active_lysimeter[cell] = h (cell) > h_lysimeter (cell);

      for (size_t edge = 0; edge != edge_size ; ++edge)
        {
          Kold[edge] = find_K_edge (soil, geo, edge, h, h_ice, h_previous, T);

void 
Horizon::initialize_base (const bool top_soil,
                          const int som_size, const double center_z, 
                          const Texture& texture, Treelog& msg)
{ 
  TREELOG_MODEL (msg);
  const double clay_lim = texture_below ( 2.0 /* [um] USDA Clay */);
  fast_clay = texture.mineral () * clay_lim;
  fast_humus = texture.humus;
  impl->initialize (*hydraulic, texture, quartz () * texture.mineral (),
                    som_size, top_soil, center_z, msg); 
  impl->secondary->initialize (msg);

  std::ostringstream tmp;
  const double h_lim = secondary_domain ().h_lim ();
  if (h_lim >= 0.0)
    impl->primary_sorption_fraction = 1.0;
  else
    {
      const double h_sat = 0.0;
      const double Theta_sat = hydraulic->Theta (h_sat);

      // Integrate h over Theta.
      PLF h_int;
      const double h_min = Hydraulic::r2h (impl->r_pore_min);
      const double Theta_min =hydraulic->Theta (h_min);
      h_int.add (Theta_min, 0.0);
      static const int intervals = 100;
      double delta = (Theta_sat - Theta_min) / (intervals + 0.0);
      double sum = 0.0;
      for (double Theta = Theta_min + delta; Theta < Theta_sat; Theta += delta)
        {
          double my_h = hydraulic->h (Theta);
          sum += my_h * delta;
          h_int.add (Theta, sum);
        }

      const double h_wp = -15000.0;
      const double Theta_wp = hydraulic->Theta (h_wp);
      const double Theta_lim =  hydraulic->Theta (h_lim);
      tmp << "A saturated secondary domain contain " 
          << 100.0 * (Theta_sat - Theta_lim) / (Theta_sat - Theta_wp)
          << " % of plant available water\n";
      impl->primary_sorption_fraction = h_int (Theta_lim) / h_int (Theta_sat);
      tmp << "Primary domain contains "
          << 100.0 * impl->primary_sorption_fraction
          << " % of the available sorption sites\n";
    }
  tmp << "h\th\tTheta\tK\n"
      << "cm\tpF\t\tcm/h\n";
  const double h_Sat = 0;
  tmp << h_Sat << "\t" << "\t" << hydraulic->Theta (h_Sat) 
      << "\t" << hydraulic->K (h_Sat) << "\n";
  const double pF_Zero = 0;
  const double h_Zero = pF2h (pF_Zero);
  tmp << h_Zero << "\t" << pF_Zero << "\t" << hydraulic->Theta (h_Zero) 
      << "\t" << hydraulic->K (h_Zero) << "\n";
  const double pF_One = 1;
  const double h_One = pF2h (pF_One);
  tmp << h_One << "\t" << pF_One << "\t" << hydraulic->Theta (h_One) 
      << "\t" << hydraulic->K (h_One) << "\n";
  const double pF_FC = 2.0;
  const double h_FC = pF2h (pF_FC);
  tmp << h_FC << "\t" << pF_FC << "\t" << hydraulic->Theta (h_FC) 
      << "\t" << hydraulic->K (h_FC) << "\n";
  const double pF_Three = 3;
  const double h_Three = pF2h (pF_Three);
  tmp << h_Three << "\t" << pF_Three << "\t" << hydraulic->Theta (h_Three) 
      << "\t" << hydraulic->K (h_Three) << "\n";
  const double pF_WP = 4.2;
  const double h_WP = pF2h (pF_WP);
  tmp << h_WP << "\t" << pF_WP << "\t" << hydraulic->Theta (h_WP) 
      << "\t" << hydraulic->K (h_WP);
  msg.debug (tmp.str ());
}

//.........这里部分代码省略.........
            tmp << "J_secondary[" << edge << "] = " << J_secondary[edge]
                << ", in_sign = " << in_sign << ", J_above = " << J_above;
            msg.bug (tmp.str ());
          }
      }
    J_sum /= geometry ().surface_area (); // [g/cm^2 S/h]
    daisy_approximate (-J_above, J_sum);
  }

  // We set a fixed concentration below lower boundary, if specified.
  std::map<size_t, double> C_border;

  const double C_below = chemical.C_below ();
  if (C_below >= 0.0)
    {
      const std::vector<size_t>& edge_below 
        = geometry ().cell_edges (Geometry::cell_below);
      const size_t edge_below_size = edge_below.size ();

      for (size_t i = 0; i < edge_below_size; i++)
        {
          const size_t edge = edge_below[i];
          C_border[edge] = C_below;
        }
    }

  // Tertiary transport.
  tertiary->solute (geometry (), soil_water, J_tertiary, dt, chemical, msg);

  // Fully adsorbed.
  if (chemical.adsorption ().full ())
    {
      static const symbol solid_name ("immobile transport");
      Treelog::Open nest (msg, solid_name);
      if (!iszero (J_above))
        {
          std::ostringstream tmp;
          tmp << "J_above = " << J_above << ", expected 0 for full sorbtion";
          msg.error (tmp.str ());
        }

      // Secondary "transport".
      std::vector<double> J2 (edge_size, 0.0); // Flux delivered by flow.
      std::vector<double> Mn (cell_size); // New content.
      for (size_t c = 0; c < cell_size; c++)
        {
          Mn[c] = chemical.M_secondary (c) + chemical.S_secondary (c) * dt;
          if (Mn[c] < 0.0)
            {
              S_extra[c] = Mn[c] / dt;
              Mn[c] = 0.0;
            }
          else
            S_extra[c] = 0.0;
        }
      chemical.set_secondary (soil, soil_water, Mn, J2);

      // Primary "transport".
      primary_transport (geometry (), soil, soil_water,
                         *matrix_solid, sink_sorbed, 0, J_primary, C_border,
                         chemical, S_extra, dt, scope, msg);
      return;
    }

  // Secondary transport activated.
  secondary_transport (geometry (), soil, soil_water, J_secondary, C_border,
                       chemical, S_extra, dt, scope, msg);

  // Solute primary transport.
  for (size_t transport_iteration = 0; 
       transport_iteration < 2; 
       transport_iteration++)
    for (size_t i = 0; i < matrix_solute.size (); i++)
      {
        solute_attempt (i);
        static const symbol solute_name ("solute");
        Treelog::Open nest (msg, solute_name, i, matrix_solute[i]->objid);
        try
          {
            primary_transport (geometry (), soil, soil_water, 
                               *matrix_solute[i], sink_sorbed, 
                               transport_iteration,
                               J_primary, C_border,
                               chemical, S_extra, dt, scope, msg);
            if (i > 0)
              msg.debug ("Succeeded");
            return;
          }
        catch (const char* error)
          {
            msg.debug (std::string ("Solute problem: ") + error);
          }
        catch (const std::string& error)
          {
            msg.debug(std::string ("Solute trouble: ") + error);
          }
        solute_failure (i);
      }
  throw "Matrix solute transport failed";
}

void
MovementSolute::primary_transport (const Geometry& geo, const Soil& soil,
                                   const SoilWater& soil_water,
                                   const Transport& transport,
                                   const bool sink_sorbed,
                                   const size_t transport_iteration,
                                   const std::map<size_t, double>& J_forced,
                                   const std::map<size_t, double>& C_border,
                                   Chemical& solute, 
                                   const std::vector<double>& S_extra,
                                   const double dt,
                                   const Scope& scope, Treelog& msg)
{ 
  
  // Edges.
  const size_t edge_size = geo.edge_size ();

  std::vector<double> q (edge_size); // Water flux [cm].
  std::vector<double> J (edge_size); // Flux delivered by flow.

  for (size_t e = 0; e < edge_size; e++)
    {
      q[e] = soil_water.q_primary (e);
      daisy_assert (std::isfinite (q[e]));
      J[e] = 0.0;
    }

  // Cells.
  const size_t cell_size = geo.cell_size ();

  std::vector<double> Theta_old (cell_size); // Water content at start...
  std::vector<double> Theta_new (cell_size); // ...and end of timestep.
  std::vector<double> C (cell_size); // Concentration given to flow.
  std::vector<double> A (cell_size); // Sorbed mass not given to flow.
  std::vector<double> S (cell_size); // Source given to flow.

  for (size_t c = 0; c < cell_size; c++)
    {
      Theta_old[c] = soil_water.Theta_primary_old (c);
      daisy_assert (Theta_old[c] > 0.0);
      Theta_new[c] = soil_water.Theta_primary (c);
      daisy_assert (Theta_new[c] > 0.0);
      C[c] = solute.C_primary (c);
      daisy_assert (C[c] >= 0.0);
      const double M = solute.M_primary (c);
      daisy_assert (M >= 0.0);
      A[c] = M - C[c] * Theta_old[c];
      daisy_assert (std::isfinite (A[c]));
      if (A[c] < 0.0)
        {
          daisy_approximate (M,  C[c] * Theta_old[c]);
          A[c] = 0.0;
        }
      daisy_assert (A[c] >= 0.0);
      S[c] = solute.S_primary (c) + S_extra[c];
      if (sink_sorbed && S[c] < 0.0)
        {
          A[c] += S[c] * dt;
          S[c] = 0.0;
          if (A[c] < 0.0)
            {
              S[c] = A[c] / dt;
              A[c] = 0.0;
            }
        }

      daisy_assert (std::isfinite (S[c]));
    }
  
  // Flow.
  transport.flow (geo, soil, Theta_old, Theta_new, q, solute.objid, 
                  S, J_forced, C_border, C, J, 
                  solute.diffusion_coefficient (), 
                  dt, msg);

  // Check fluxes.
  for (size_t e = 0; e < edge_size; e++)
    daisy_assert (std::isfinite (J[e]));


  // Update with new content.
  std::vector<double> M (cell_size);
  for (size_t c = 0; c < cell_size; c++)
    {
      daisy_assert (std::isfinite (C[c]));
      M[c] = A[c] + C[c] * Theta_new[c];

      if (M[c] < 0.0)
        {
          std::ostringstream tmp;
          tmp << "M[" << c << "] = " << M[c] 
              << " @ " << geo.cell_name (c)
              << ", C = " << C[c]
              << ", A = " << A[c]
              << ", M_new = " << M[c]
              << ", M_old = " << solute.M_primary (c) << ", dt " << dt
              << ", S = " << S[c] 
              << ", S_extra = " << S_extra[c];
          solute.debug_cell (tmp, c);
          tmp << ", Theta_old " << Theta_old[c]
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  const double estar = FAO::SaturationVapourPressure (Tl); // [Pa]
  daisy_assert (gbw_ms >= 0.0);
  gbw = Resistance::ms2molly (Tl, Ptot, gbw_ms);
  daisy_assert (gbw >= 0.0);
  // CAI in each interval.
  const double dCAI = PAR_LAI / No;
  
  for (int i = 0; i < No; i++)
    {
      const double height = PAR_height[i+1];
      daisy_assert (height < PAR_height[i]);

      // Leaf Area index for a given leaf layer
      const double LA = prevLA - LAIvsH (height);
      daisy_assert (LA >= 0.0);

      prevLA = LAIvsH (height);
      accCAI += LA;

      if (LA * fraction [i] > 0)
	{  
	  // PAR in mol/m2/s = PAR in W/m2 * 0.0000046
	  const double dPAR = (PAR[i] - PAR[i+1])/dCAI * 0.0000046; //W/m2->mol/m²leaf/s

          if (dPAR < 0)
            {
              std::stringstream tmp;
              tmp << "Negative dPAR (" << dPAR
                  << " [mol/m^2 leaf/h])" << " PAR[" << i << "] = " << PAR[i]
                  << " PAR[" << i+1 << "] = " << PAR[i+1] 
                  << " dCAI = " << dCAI
                  << " LA = " << LA
                  << " fraction[" << i << "] = " << fraction [i];
              msg.debug (tmp.str ());
              continue;
            }

	  // log variable
	  PAR_ += dPAR * dCAI * 3600.0; //mol/m²area/h/fraction

	  // Photosynthetic rubisco capacity 
	  const double vmax25 = crop_Vm_total[i]*fraction[i];//[mol/m²leaf/s/fracti.]
	  daisy_assert (vmax25 >= 0.0);

	  // leaf respiration
	  const double rd = respiration_rate(vmax25, Tl);
	  daisy_assert (rd >= 0.0);

	  //solving photosynthesis and stomatacondctance model for each layer
          double& pn = pn_vector[i];
	  double ci  = 0.5 * CO2_atm;//first guess for ci, [Pa]
          double& hs = hs_vector[i];
          hs = 0.5;              // first guess of hs []
          double& cs = cs_vector[i];
          //first gues for stomatal cond,[mol/s/m²leaf]
	  double gsw = Stomatacon->minimum () * 2.0;
	  // double gsw = 2 * b; // old value
	  const int maxiter = 150;
	  int iter = 0;
	  double lastci;
          double lasths;
          double lastgs;
	  do
	    {
	      lastci = ci; //Stomata CO2 pressure 
              lasths = hs;

void 
Movement1D::tick_water (const Geometry1D& geo,
                        const Soil& soil, const SoilHeat& soil_heat, 
                        Surface& surface, Groundwater& groundwater,
                        const std::vector<double>& S,
                        std::vector<double>& h_old,
                        const std::vector<double>& Theta_old,
                        const std::vector<double>& h_ice,
                        std::vector<double>& h,
                        std::vector<double>& Theta,
                        std::vector<double>& q,
                        std::vector<double>& q_p,
                        const double dt, 
                        Treelog& msg)
{
  const size_t top_edge = 0U;
  const size_t bottom_edge = geo.edge_size () - 1U;

  // Limit for groundwater table.
  size_t last  = soil.size () - 1;

  // Limit for ridging.
  const size_t first = (surface.top_type (geo, 0U) == Surface::soil)
    ? surface.last_cell (geo, 0U) 
    : 0U;
  // Calculate matrix flow next.

  for (size_t m = 0; m < matrix_water.size (); m++)
    {
      water_attempt (m);
      Treelog::Open nest (msg, matrix_water[m]->name);
      try
        {
          matrix_water[m]->tick (msg, geo, soil, soil_heat,
                                 first, surface, top_edge, 
                                 last, groundwater, bottom_edge,
                                 S, h_old, Theta_old, h_ice, h, Theta, 0U, q, 
                                 dt);

          for (size_t i = last + 2; i <= soil.size (); i++)
            {
              q[i] = q[i-1];
              q_p[i] = q_p[i-1];
            }

          // Update surface and groundwater reservoirs.
          surface.accept_top (q[0] * dt, geo, 0U, dt, msg);
          surface.update_pond_average (geo);
          const double q_down = q[soil.size ()] + q_p[soil.size ()];
          groundwater.accept_bottom (q_down * dt, geo, soil.size ());
          if (m > 0)
            msg.debug ("Reserve model succeeded");
          return;
        }
      catch (const char* error)
        {
          msg.debug (std::string ("UZ problem: ") + error);
        }
      catch (const std::string& error)
        {
          msg.debug (std::string ("UZ trouble: ") + error);
        }
      
      water_failure (m);
    }
  throw "Water matrix transport failed"; 
}