// biopore_drain.C --- Static vertical biopores that ends in drain pipes.
// 
// Copyright 2008 Per Abrahamsen and KU.
//
// This file is part of Daisy.
// 
// Daisy is free software; you can redistribute it and/or modify
// the Free Software Foundation; either version 2.1 of the License, or
// (at your option) any later version.
// 
// Daisy is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser Public License for more details.
// 
// You should have received a copy of the GNU Lesser Public License
// along with Daisy; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#define BUILD_DLL

#include "biopore.h"
#include "block_model.h"
#include "vcheck.h"
#include "librarian.h"
#include "submodeler.h"
#include "check.h"
#include "geometry.h"
#include "soil.h"
#include "anystate.h"
#include "chemical.h"
#include "groundwater.h"
#include "treelog.h"
#include "frame.h"
#include "assertion.h"
#include "mathlib.h"
#include <sstream>

// The 'drain' model.

struct BioporeDrain : public Biopore
{
  // Parameters.
  /* const */ double pipe_position;   // [cm]

  // State
  Anystate get_state () const
  { return Anystate::none (); }
  void set_state (const Anystate&)
  { }

  // Simulation.
  double total_water () const
  { return 0.0; }
  void get_solute (IM&) const
  { }
  double air_bottom (size_t) const    // Lowest point with air [cm]
  { return pipe_position; }
  
  double infiltration_capacity (const Geometry& geo, size_t e, 
                                const double dt) const
  { 
    // Any macropores that reach the surface.
    if (height_start < 0.0)
      return 0.0;

    return max_infiltration_rate (geo, e) * dt; 
  }

  double matrix_biopore_drain (size_t c, const Geometry& geo, 
                               bool active, 
                               const double h_barrier, double pressure_limit,
                               double K_xx, double h) const;
  void forward_sink (const Geometry& geo,    
                     const std::vector<bool>& active,
                     const std::vector<double>& K, 
                     const std::vector<double>& K_crack, 
                     const double h_barrier,
                     const double pressure_limit,
                     const std::vector<double>& h, 
                     std::vector<double>& S3) const;
  void tick_source (const Geometry&, const std::vector<bool>&,
                    const std::vector<double>&)
  { }
  void update_matrix_sink (const Geometry& geo,    
                           const std::vector<bool>& active,
                           const std::vector<double>& K,
                           const std::vector<double>& K_crack,
                           const double h_barrier,
                           const double pressure_limit,
                           const std::vector<double>& h,
                           const double dt);
  void update_water ()
  { }
  void add_to_sink (std::vector<double>&,
                    std::vector<double>& S_drain) const
  {
    const size_t cell_size = S.size ();
    daisy_assert (S_drain.size () == cell_size);
    for (size_t c = 0; c < cell_size; c++)
      S_drain[c] += S[c];
  }
  void add_solute (symbol, size_t, const double)
  { }
  void matrix_solute (const Geometry&, const double, 
                      const Chemical&, std::vector<double>&,
                      Treelog&)
  { /* Handled by S_drain. */ }
  void output (Log& log) const
  { output_base (log); }

  // Create and Destroy.
  bool initialize (const Units& units, const Geometry& geo, const Scope& scope,
                   const Groundwater& groundwater, Treelog& msg)
  {
    bool ok = initialize_base (units, geo, scope, msg); 
    if (pipe_position > 0)
      {
        if (groundwater.is_pipe ())
          // Pipe height not specified here, use value from column.
          pipe_position = groundwater.pipe_height ();
        else
          {
            msg.error ("Unknown pipe position");
            ok = false;
          }
      }
    return ok;
  }
  bool check (const Geometry& geo, Treelog& msg) const
  { return check_base (geo, msg); }
  BioporeDrain (const BlockModel& al);
};

double 
BioporeDrain::matrix_biopore_drain (size_t c, const Geometry& geo, 
                                    bool active, 
                                    const double h_barrier, 
                                    const double pressure_limit,
                                    double K_xx, double h) const
{
  const double M_c = density_cell[c];
  if (!std::isnormal (M_c))
    // No biopores here.
    return 0.0;
  const double r_c = diameter / 2.0;
  const double h_3 = air_bottom (c) - geo.cell_z (c);

#if 0
  //--------------------------------------
  std::ostringstream tmp;
  tmp << "c = " << c 
      << " airbottom = " << air_bottom (c) 
      << " cell_z = " <<  geo.cell_z (c)
      << " h = " << h
      << " h_3 = " << h_3 
      << " h + p_end = " << h+pressure_limit; 
  Assertion::message (tmp.str ());  
  //--------------------------------------
#endif


  double S;
  if (active && h>h_3 + h_barrier)
    {
      // The largest pressure gradient between the domains are
      // pressure_limit, above that we claim air will disrupt the suction.
      const double h_3_suck = std::max (h_3, h + pressure_limit);

      S = matrix_to_biopore (K_xx, M_c, r_c, h, h_3_suck)
        * geo.fraction_in_z_interval (c, height_start, height_end);
    }
  else 
    S = 0.0;
  return S;
}

void 
BioporeDrain::forward_sink (const Geometry& geo,    
                            const std::vector<bool>& active,
                            const std::vector<double>& K, 
                            const std::vector<double>& /* K_crack */, 
                            const double h_barrier,
                            const double pressure_limit,
                            const std::vector<double>& h, 
                            std::vector<double>& S3) const
{
  const size_t cell_size = geo.cell_size ();
  for (size_t c = 0; c < cell_size; c++)
    S3[c] += matrix_biopore_drain (c, geo, active[c], h_barrier, 
                                   pressure_limit, K[c], h[c]);
}

void
BioporeDrain::update_matrix_sink (const Geometry& geo,    
                                  const std::vector<bool>& active,
                                  const std::vector<double>& K, 
                                  const std::vector<double>& K_crack, 
                                  const double h_barrier,
                                  const double pressure_limit,
                                  const std::vector<double>& h, 
                                  const double /* dt */)
{
  std::fill (S.begin (), S.end (), 0.0);
  forward_sink (geo, active, K, K_crack, h_barrier, pressure_limit, h, S);
}

BioporeDrain::BioporeDrain (const BlockModel& al)
  : Biopore (al),
    pipe_position (al.number ("pipe_position", 42.42e42))
{ }

static struct BioporeDrainSyntax : DeclareModel
{
  Model* make (const BlockModel& al) const
  { return new BioporeDrain (al); }

  BioporeDrainSyntax ()
    : DeclareModel (Biopore::component, "drain", "\
Biopores that ends in the drain pipes.")
  { }
  void load_frame (Frame& frame) const
  { 
    frame.declare ("pipe_position", "cm", Check::negative (), Attribute::Const,
                "Height pipes are placed in the soil (a negative number).\n\
By default, use the height specified for pipes in the column.");
  }
} BioporeDrain_syntax;

// biopore_drain.C ends here.
