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

void
Fetch::initialize (std::vector<Select*>& select, Treelog& msg)
{ 
  Treelog::Open nest (msg, this->tag);
  bool found = false;
  
  for (size_t j = 0; j != select.size (); j++)
    {
      Treelog::Open nest (msg, select[j]->tag ());

      if (this->tag == select[j]->tag ())
        {
          if (found)
            msg.warning ("Duplicate tag ignored");
          else
            {	
              select[j]->add_dest (this);
              this->select_dimension = select[j]->dimension ();
              this->type = ((select[j]->handle != Select::Handle::current)
                            && !select[j]->accumulate)
                ? Fetch::Flux 
                : Fetch::NewContent;
              found = true;
            }
        }
    }
  if (!found)
    msg.warning ("No tag found");
}

void 
MovementSolute::element (const Soil& soil, const SoilWater& soil_water,
                         DOE& element, 
                         const double diffusion_coefficient, double dt, 
                         Treelog& msg)
{
  for (size_t i = 0; i < matrix_solute.size (); i++)
    {
      Treelog::Open nest (msg, "element", i, matrix_solute[i]->library_id ());
      try
        {
          matrix_solute[i]->element (geometry (), soil, soil_water, element, 
                                     diffusion_coefficient, dt, 
                                     msg);
          if (i > 0)
            msg.message ("Succeeded");
          return;
        }
      catch (const char* error)
        {
          msg.warning (std::string ("DOM problem: ") + error);
        }
      catch (const std::string& error)
        {
          msg.warning (std::string ("DOM trouble: ") + error);
        }
    }
  throw "Matrix element transport failed";
}

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

bool
GnuplotProfile::plot (std::ostream& out, Treelog& msg)
{ 
  if (xplus.size () < 1 || zplus.size () < 1)
    {
      msg.warning ("Nothing to plot");
      return false;
    }

  // Header.
  plot_header (out);
  out << "\
set pm3d map\n\
set pm3d corners2color c4\n\
unset colorbox\n";
  // Same size axes.
  out << "\
set size ratio -1\n";

  // Legend.
  if (legend != "auto")
    out << "set key " << legend_table[legend] << "\n";

  // Extra.
  for (size_t i = 0; i < extra.size (); i++)
    out << extra[i].name () << "\n";

  // Plot.
  out << "splot '-' using 2:1:3 title \"\"\n";

  // Data.
  daisy_assert (value.size () == xplus.size () * zplus.size ());

  // Cell corners only.
  daisy_assert (value.size () > 0);
  out << "0 0 " << value[0] << "\n";
  daisy_assert (value.size () >= zplus.size ());
  for (size_t iz = 0; iz < zplus.size (); iz++)
    out << zplus[iz] << " 0 " << value[iz * xplus.size ()] << "\n";
  for (size_t ix = 0; ix < xplus.size (); ix++)
    {
      out << "\n0 " << xplus[ix] << " " << value[ix] << "\n";
      for (size_t iz = 0; iz < zplus.size (); iz++)
        {
          const size_t c = ix + iz * xplus.size ();
          daisy_assert (c < value.size ());
          out << zplus[iz] << " " << xplus[ix] << " " << value[c] << "\n";
        }
    }

  out << "e\n";

  // The end.
  if (interactive ())
    out << "pause mouse\n";

  return true;
}

void
ABAProdUptake::production (const Units& units,
                           const Geometry& geo, const SoilWater& soil_water,
                           const std::vector<double>& S /* [cm^3/cm^3/h] */,
                           const std::vector<double>& /* l [cm/cm^3] */,
                           std::vector<double>& ABA    /* [g/cm^3/h] */,
                           Treelog& msg) const
{
  // Check input.
  const size_t cell_size = geo.cell_size ();
  daisy_assert (ABA.size () == cell_size);
  daisy_assert (S.size () == cell_size);
  
  // For all cells.
  for (size_t c = 0; c < cell_size; c++)
    {
      const double h = soil_water.h (c);
      // Set up 'h' in scope.
      scope.set (h_name, soil_water.h (c));

      // Find soil value.
      double value = 0.0;
      if (!expr->tick_value (units, value, ABA_unit, scope, msg))
	msg.error ("No ABA production value found");
      if (!std::isfinite (value) || value < 0.0)
        {
          if (h > -14000) // We are not concerned with ABA below wilting point.
            {
              std::ostringstream tmp;
              tmp << "ABA in cell " << c << " with h = " << h
                  << " was " << value << " [" << ABA_unit << "], using 0 instead";
              msg.warning (tmp.str ());
            }
          value = 0.0;
        }
      daisy_assert (std::isfinite (S[c]));

      // Find ABA uptake.
      ABA[c] = value * S[c];
      // [g/cm^3 S/h] = [g/cm^3 W] * [cm^3 W/cm^3 S/h]
      daisy_assert (std::isfinite (ABA[c]));
    }
}

void
Rootdens_G_P::set_density  (const Geometry& geo,
                            double SoilDepth, double CropDepth,
                            const double /* CropWidth [cm] */,
                            double WRoot, double /* DS */,
                            std::vector<double>& Density, Treelog& msg)
{
    const double Depth = std::min (SoilDepth, CropDepth);
    // Dimensional conversion.
    static const double m_per_cm = 0.01;

    const double MinLengthPrArea = (DensRtTip * 1.2) * CropDepth;
    const double LengthPrArea
        = std::max (m_per_cm * SpRtLength * WRoot, MinLengthPrArea); // [cm/cm^2]

    // We find a * depth first.
    const double ad = density_distribution_parameter (LengthPrArea /
                      (CropDepth * DensRtTip));
    // We find L0 from a d.
    //
    // L0 * exp (-a d) = DensRtTip
    // => L0 = DensRtTip / exp (-a d)
    L0 = DensRtTip * exp (ad);	//  1 / exp (-x) = exp (x)
    a = ad / CropDepth;

    if (Depth < CropDepth)
    {
        double Lz = L0 * exp (-a * Depth);
        a = density_distribution_parameter (LengthPrArea / (Depth * Lz)) / Depth;
    }

    // Check minimum density
    double extra = 0.0;
    if (MinDens > 0.0 && WRoot > 0.0)
    {
        daisy_assert (L0 > 0.0);
        daisy_assert (a > 0.0);
        const double too_low = -log (MinDens / L0) / a; // [cm]

        if (too_low < Depth)
        {
            // We don't have MinDens all the way down.
            const double NewLengthPrArea
                =  LengthPrArea - MinDens * Depth; // [cm/cm^2]
#if 1
            Treelog::Open nest (msg, "RootDens G+P");
            std::ostringstream tmp;
            tmp << "too_low = " << too_low
                << ", NewLengthPrArea = " << NewLengthPrArea
                << "MinLengthPrArea = " << MinLengthPrArea;
            msg.warning (tmp.str ());
#endif
            if (too_low > 0.0 && NewLengthPrArea > too_low * DensRtTip * 1.2)
            {
                // There is enough to have MinDens all the way, spend
                // the rest using the standard model until the point
                // where the standard model would give too little..
                a = density_distribution_parameter (NewLengthPrArea
                                                    / (too_low * DensRtTip));
                L0 = DensRtTip * exp (a);
                a /= too_low;
                extra = MinDens;
            }
            else
            {
                // There is too little, use uniform density all the way.
                L0 = 0.0;
                extra = LengthPrArea / Depth;
            }
        }
    }

    const size_t size = geo.cell_size ();
    daisy_assert (Density.size () == size);
#if 1
    for (size_t i = 0; i < size; i++)
    {
        const double f = geo.fraction_in_z_interval (i, 0.0, -Depth);
        const double d = -geo.cell_z (i);
        if (f > 0.01)
            Density[i] = f * (extra + L0 * exp (- a * d));
        else
            Density[i] = 0.0;
    }
#else // 0
    unsigned int i = 0;
    for (; i == 0 || -geo.zplus (i-1) < Depth; i++)
    {
        daisy_assert (i < geo.size ());
        Density[i] = extra + L0 * exp (a * geo.cell_z (i));
    }

    for (; i < geo.size (); i++)
        Density[i] = 0.0;
#endif // 0
}

//.........这里部分代码省略.........
#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;  

	  // Q_Cw_h is shorthand for Qmatrix * Cw * h
	  const ublas::vector<double> Q_Cw_h = prod (Q_Cw, h);

	  //Initialize b-vector
	  ublas::vector<double> b (cell_size);  
	  //b = sumvec + (1.0 / ddt) * (Qmatrix * Cw * h + Qmatrix *(Wxx-Wyy));
	  b = sumvec + (1.0 / ddt) * (Q_Cw_h
				      + prod (Qmat, Theta_previous-Theta));

	  // Force active drains to zero h.
          drain (geo, drain_cell, drain_water_level,
		 h, Theta_previous, Theta, S_vol,
#ifdef TEST_OM_DEN_ER_BRUGT
                 S_macro,
#endif
                 dq, ddt, drain_cell_on, A, b, debug, msg);  
          
          try {
            solver->solve (A, b, h); // Solve Ah=b with regard to h.
          } catch (const char *const error) {
              std::ostringstream tmp;
              tmp << "Could not solve equation system: " << error;
              msg.warning (tmp.str ());
              // Try smaller timestep.
              iterations_used = max_loop_iter + 100;
              break;
          }

	  for (int c=0; c < cell_size; c++) // update Theta 
	    Theta (c) = soil.Theta (c, h (c), h_ice (c)); 

	  if (debug > 1)
	    {
	      std::ostringstream tmp;
	      tmp << "Time left = " << time_left << ", ddt = " << ddt 
		  << ", iteration = " << iterations_used << "\n";
	      tmp << "B = " << B << "\n";
	      tmp << "h = " << h << "\n";
	      tmp << "Theta = " << Theta;
	      msg.message (tmp.str ());
	    }
          
          for (int c=0; c < cell_size; c++)
            {
              if (h (c) < min_pressure_potential || h (c) > max_pressure_potential)
                {
                  std::ostringstream tmp;
                  tmp << "Pressure potential out of realistic range, h[" 
                      << c << "] = " << h (c);
                  msg.debug (tmp.str ());
                  iterations_used = max_loop_iter + 100;
                  break;
                } 
            }
        }

void 
MovementSolute::secondary_flow (const Geometry& geo, 
                                const std::vector<double>& Theta_old,
                                const std::vector<double>& Theta_new,
                                const std::vector<double>& q,
                                const symbol name,
                                const std::vector<double>& S, 
                                const std::map<size_t, double>& J_forced,
                                const std::map<size_t, double>& C_border,
                                std::vector<double>& M, 
                                std::vector<double>& J, 
                                const double dt,
                                Treelog& msg)
{
  const size_t cell_size = geo.cell_size ();
  const size_t edge_size = geo.edge_size ();
  
  // Full timstep left.
  daisy_assert (dt > 0.0);
  double time_left = dt;

  // Initial water content.
  std::vector<double> Theta (cell_size);
  for (size_t c = 0; c < cell_size; c++)
    Theta[c] = Theta_old[c];

  // Small timesteps.
  for (;;)
    {
      // Are we done yet?
      const double min_timestep_factor = 1e-19;
      if (time_left < 0.1 * min_timestep_factor * dt)
        break;

      // Find new timestep.
      double ddt = time_left;
  
      // Limit timestep based on water flux.
      for (size_t e = 0; e < edge_size; e++)
        {
          const int cell = (q[e] > 0.0 ? geo.edge_from (e) : geo.edge_to (e));
          if (geo.cell_is_internal (cell) 
              && Theta[cell] > 1e-6 && M[cell] > 0.0)
            {
              const double loss_rate = std::fabs (q[e]) * geo.edge_area (e);
              const double content = Theta[cell] * geo.cell_volume (cell); 
              const double time_to_empty = content / loss_rate;
              if (time_to_empty < min_timestep_factor * dt)
                {
                  msg.warning ("Too fast water movement in secondary domain");
                  ddt = min_timestep_factor * dt;
                  break;
                }
              
              // Go down in timestep while it takes less than two to empty cell.
              while (time_to_empty < 2.0 * ddt)
                ddt *= 0.5;
            }
        }

      // Cell source.  Must be before transport to avoid negative values.
      for (size_t c = 0; c < cell_size; c++)
        M[c] += S[c] * ddt;

      // Find fluxes using new values (more stable).
      std::vector<double> dJ (edge_size, -42.42e42);
      for (size_t e = 0; e < edge_size; e++)
        {
          std::map<size_t, double>::const_iterator i = J_forced.find (e);
          if (i != J_forced.end ())
            // Forced flux.
            {
              dJ[e] = (*i).second;
              daisy_assert (std::isfinite (dJ[e]));
              continue;
            }

          const int edge_from = geo.edge_from (e);
          const int edge_to = geo.edge_to (e);
          const bool in_flux = q[e] > 0.0;
          int flux_from = in_flux ? edge_from : edge_to;
          double C_flux_from = -42.42e42;

          if (geo.cell_is_internal (flux_from))
            // Internal cell, use its concentration.
            {
              if (Theta[flux_from] > 1e-6 && M[flux_from] > 0.0)
                // Positive content in positive water.
                C_flux_from = M[flux_from] / Theta[flux_from];
              else
                // You can't cut the hair of a bald guy.
                C_flux_from = 0.0;
            }
          else
            {
              i = C_border.find (e);
              if (i != C_border.end ())
                // Specified by C_border.
                C_flux_from = (*i).second;
              else
//.........这里部分代码省略.........

 void doIt (Daisy&, const Scope&, Treelog& out)
 { 
   out.warning (message.name ());
 }

//.........这里部分代码省略.........
	      h[i] = std::max (h_old[i], 0.0);
              daisy_assert (std::isfinite (h[i]));
	    }
	  else
	    {
	      q[i+1] = -K_new;
	      Theta[i] = Theta_next;
	      h[i] = soil.h (i, Theta[i]);
              daisy_assert (std::isfinite (h[i]));
	    }
	}
      daisy_assert (std::isfinite (h[i]));
      daisy_assert (std::isfinite (Theta[i]));
      daisy_assert (std::isfinite (q[i+1]));
      daisy_assert (Theta[i] <= Theta_sat);
      daisy_assert (Theta[i] > Theta_res);
    }

  // Lower border.
  q_down = q[last + 1];

  if (bottom.bottom_type () == Groundwater::forced_flux)
    // Ensure forced bottom.
    {
      double extra_water = (bottom.q_bottom (bottom_edge) - q_down) * dt;

      for (int i = last; true; i--)
	{
	  q[i+1] += extra_water / dt;
	  if (i < static_cast<int> (first))
            {
              if (extra_water > 0.0 && overflow_warn)
                {
                  msg.warning ("Soil profile saturated, water flow to surface");
                  overflow_warn = false;
                }
              break;
            }
	  const double dz = geo.dz (i);
          const double h_min = pF2h (10.0);
          const double Theta_min = soil.Theta (i, h_min, h_ice[i]);
	  const double Theta_sat = soil.Theta (i, 0.0, h_ice[i]);
	  Theta[i] += extra_water / dz;
          if (Theta[i] <= Theta_min)
            {
              extra_water = (Theta[i] - Theta_min) * dz;
	      Theta[i] = Theta_min;
	      h[i] = h_min;
            }
	  else if (Theta[i] <= Theta_sat)
	    {
	      extra_water = 0;
	      h[i] = soil.h (i, Theta[i]);
              break;
	    }
	  else
	    {
	      extra_water = (Theta[i] - Theta_sat) * dz;
	      Theta[i] = Theta_sat;
	      h[i] = 0.0;
	    }
	}
      q_up = q[first];
      q_down = q[last + 1];
    }

//.........这里部分代码省略.........
                      pn, ci, dPAR /*[mol/m²leaf/s]*/, 
                      gsw, gbw, Tl, vmax25, rd, msg);//[mol/m²leaf/s/fraction]

              // Vapour pressure at leaf surface. [Pa]
              const double es = (gsw * estar + gbw * ec) / (gsw + gbw);
              // Vapour defecit at leaf surface. [Pa]
              const double Ds = bound (0.0, estar - es, estar);
              // Relative humidity at leaf surface. []
              hs = es / estar;
              const double hs_use = bound (0.0, hs, 1.0);

              // Boundary layer resistance. [s*m2 leaf/mol]
              daisy_assert (gbw >0.0);
              const double rbw = 1./gbw;   //[s*m2 leaf/mol]
  
              // leaf surface CO2 [Pa]

              // We really should use CO2_canopy instead of CO2_atm
              // below.  Adding the resitence from canopy point to
              // atmostphere is not a good workaround, as it will
              // ignore sources such as the soil and stored CO2 from
              // night respiration.
              cs = CO2_atm - (1.4 * pn * Ptot * rbw); //[Pa] 
              daisy_assert (cs > 0.0);

              //stomatal conductance
              gsw = Stomatacon->stomata_con (ABA /*g/cm^3*/,
                                             h_x /* MPa */, 
                                             hs_use /*[]*/,
                                             pn /*[mol/m²leaf/s]*/, 
                                             Ptot /*[Pa]*/, 
                                             cs /*[Pa]*/, Gamma /*[Pa]*/, 
                                             Ds/*[Pa]*/,
                                             msg); //[mol/m²leaf/s] 

	      iter++;
	      if(iter > maxiter)
		{
		  std::ostringstream tmp;
		  tmp << "total iterations in assimilation model exceed "
                      << maxiter;
		  msg.warning (tmp.str ());
		  break;
		}
	    }
	  // while (std::fabs (lastci-ci)> 0.01);
          while (std::fabs (lastci-ci)> 0.01
                 || std::fabs (lasths-hs)> 0.01
                 || std::fabs (lastgs-gs)> 0.01);

	  // Leaf brutto photosynthesis [gCO2/m2/h] 
	  /*const*/ double pn_ = (pn+rd) * molWeightCO2 * 3600.0;//mol CO2/m²leaf/s->g CO2/m²leaf/h
	  const double rd_ = (rd) * molWeightCO2 * 3600.0;   //mol CO2/m²/s->g CO2/m²/h
	  const double Vm_ = V_m(vmax25, Tl); //[mol/m² leaf/s/fraction]
	  const double Jm_ = J_m(vmax25, Tl); //[mol/m² leaf/s/fraction]

	  if (pn_ < 0.0)
            {
              std::stringstream tmp;
              tmp << "Negative brutto photosynthesis (" << pn_ 
                  << " [g CO2/m²leaf/h])" << " pn " << pn << " rd " << rd
                  << " CO2_atm " << CO2_atm << "  O2_atm " <<  O2_atm 
                  << "  Ptot " <<  Ptot << "  pn " <<  pn << "  ci " <<  ci 
                  << "  dPAR " <<  dPAR << "  gsw " <<  gsw 
                  << "  gbw " <<  gbw << "  Tl " <<  Tl 
                  << "  vmax25 " <<  vmax25 << "  rd " <<  rd; 
              msg.error (tmp.str ());
              pn_ = 0.0;
            }
	  Ass_ += LA * pn_; // [g/m²area/h] 
	  Res += LA * rd_;  // [g/m²area/h] 
	  daisy_assert (Ass_ >= 0.0);

	  //log variables:
	  Ass_vector[i]+= pn_* (molWeightCH2O / molWeightCO2) * LA;//[g CH2O/m²area/h]
	  Nleaf_vector[i]+= rubisco_Ndist[i] * LA * fraction[i]; //[mol N/m²area]OK
	  gs_vector[i]+= gsw /* * LA * fraction[i] */;     //[mol/m² area/s]
	  ci_vector[i]+= ci /* * fraction[i] */;  //[Pa] OK
	  Vm_vector[i]+= Vm_ * 1000.0 * LA * fraction[i]; //[mmol/m² area/s]OK
	  Jm_vector[i]+= Jm_ * 1000.0 * LA * fraction[i]; //[mmol/m² area/s]OK
	  LAI_vector[i] += LA * fraction[i];//OK
	  
	  ci_middel += ci * fraction[i]/(No + 0.0);// [Pa]   OK
	  gs += LA * gsw * fraction[i]; 
	  Ass += LA * pn_ * (molWeightCH2O / molWeightCO2);//[g CH2O/m2 area/h] OK
	  LAI += LA * fraction[i];//OK
	  Vmax += 1000.0 * LA * fraction[i] * Vm_;   //[mmol/m² area/s]
	  jm += 1000.0 * LA * fraction[i] * Jm_;     //[mmol/m² area/s]
	  leafPhotN += rubisco_Ndist[i] * LA *fraction[i]; //[mol N/m²area]; 
	  fraction_total += fraction[i]/(No + 0.0);
	}
    }
  daisy_assert (approximate (accCAI, canopy.CAI));
  daisy_assert (Ass_ >= 0.0);

  // Omregning af gs(mol/(m2s)) til gs_ms (m/s) foretages ved 
  // gs_ms = gs * (R * T)/P:
  gs_ms = Resistance::molly2ms (Tl, Ptot, gs);
  return (molWeightCH2O / molWeightCO2)* Ass_;    // Assimilate [g CH2O/m2/h]
}