<commit_before>// biopore.C --- A single class of biopores.
// 
// Copyright 2008 <NAME> and KU.
//
// This file is part of Daisy.
// 
// Daisy is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser Public License as published by
// 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 "frame.h"
#include "librarian.h"
#include "scope_multi.h"
#include "scope_id.h"
#include "units.h"
#include "check.h"
#include "geometry.h"
#include "log.h"
#include "treelog.h"
#include "assertion.h"
#include "mathlib.h"
#include <sstream>

// biopore component.

const char *const Biopore::component = "biopore";

symbol
Biopore::library_id () const
{
  static const symbol id (component);
  return id;
}

symbol
Biopore::x_symbol ()
{
  static const symbol x ("x");
  return x;
}

double
Biopore::max_infiltration_rate (const Geometry& geo, size_t e) const // [cm/h]
{
  const size_t cell = geo.edge_other (e, Geometry::cell_above);
  daisy_assert (cell < geo.cell_size ());

  // Based on Poiseuille equation.
  // Q = pi r^4 rho g L / ( 8 L mu)
  //
  //        Q [cm^3/h]    Infiltration rate
  //  delta P [cm]        Pressure drop.
  //        L [cm]        Length op cylinder.
  //       mu [cm^3/cm/h] Dynamic viscosity.
  //        r [cm]        Biopore radius.
  //        M [cm^-2]     Macropore density.
  //     

  const double rho = 1.0;   // [g/cm^3]
  const double g = 9.81;      // [m/s^2]
  const double g1 = g * 100.0 * 3600.0 * 3600.0;        // [cm/h^2]
  const double mu = (1.0020 / 100.0) * 3600.0;          // [g/cm/h] At 20 dg C. 
  const double M = density (cell);                      // [cm^-2]
  const double r = diameter * 0.5;                      // [cm]
  const double r4 = r * r * r * r;                      // [cm^4]
  const double Q = M_PI * r4 * rho * g1 / (8 * mu);     // [cm^3/h]
  const double q = Q * M;                               // [cm/h]
  
  return q;
}

void 
Biopore::infiltrate (const Geometry& geo, size_t e,
                     const double amount, const double dt)
{
  daisy_assert (std::isfinite (amount));
  const double edge_area = geo.edge_area (e);
  const double total_area = geo.surface_area ();
  const double edge_flux = amount / dt;
  const double total_flux = edge_flux * edge_area / total_area;
  infiltration += total_flux;
}

void 
Biopore::solute_infiltrate (const symbol chem, 
                            const Geometry& geo, const size_t e,
                            const double amount /* [g] */, 
                            const double dt)
{
  daisy_assert (std::isfinite (amount));
  const double total_area = geo.surface_area ();
  solute_infiltration.add_value (chem, solute_infiltration.unit (), 
                                 amount / total_area / dt);
}

double 
Biopore::top_density (const size_t c) const
{
  if (height_start < 0.0)
    return 0.0;
  
  return density (c);
}

void 
Biopore::clear ()
{ 
  infiltration = 0.0; 
  solute_infiltration.clear ();
}

double 
Biopore::matrix_to_biopore (double K_xx, double M_c, double r_c, 
                            double h, double h_3)
{
  const double S = - 4*M_PI*M_c*K_xx*(h-h_3) / 
    (log(M_PI*M_c*r_c*r_c));
  daisy_assert (std::isfinite (S));
  
  if (S < 0.0)
    {
      std::ostringstream tmp;
      tmp << "Bug: S = " << S << " (should be positive), M_c = " << M_c
          << ", r_c = " << r_c << ", K_xx = " << K_xx << ", h = " << h
          << ", h_3 = " << h_3 << ", pi M_c r_c^2 = " << M_PI*M_c*r_c*r_c;
      Assertion::error (tmp.str ());
      return 0.0;
    }
  return S;
}

double 
Biopore::biopore_to_primary (const double K_matrix,
                             const double K_wall_rel, //  K_wall_rel []
                             const double M_c, const double r_c,
                             const double h, const double h_3)
{
  const double a1 = 1.1;  // Relative biopore wall size: r_i / r_c [] 
  const double a2 = K_wall_rel; // Relative biopore wall conductivity []
  const double A1 
    = a2 * (0.5 * log (M_PI * M_c * sqr (r_c)) - std::log (a1))
    / std::log (a1);
  const double B1 = A1 / (1.0 + A1);
  const double K = B1 * K_matrix;
  const double r_matrix = std::pow (M_PI * M_c, -0.5);
  const double S = M_c * 2.0 * M_PI * K * (h_3 - h)
    / std::log (r_c / r_matrix);
  if (S < 0.0)
    {
      std::ostringstream tmp;
      tmp << "K_matrix = " << K_matrix << " [cm/h]\n"
          << "K_wall_rel = " << K_wall_rel << " []\n"
          << "M_c = " << M_c << " [cm^-2]\n"
          << "r_c = " << r_c << " [cm]\n"
          << "h = " << h << " [cm]\n"
          << "h_3 = " << h_3 << " [cm]\n"
          << "a1 = " << a1 << " []\n"
          << "a2 = " << a2 << " []\n"
          << "A1 = " << A1 << " []\n"
          << "B1 = " << B1 << " []\n"
          << "K = " << K << " [cm/h]\n"
          << "r_matrix = " << r_matrix << " [cm]\n"
          << "S = " << S << " [] (expected >= 0)";
      daisy_bug (tmp.str ());
    }
  return S;
}
<commit_msg>make<commit_after>// biopore.C --- A single class of biopores.
// 
// Copyright 2008 <NAME> and KU.
//
// This file is part of Daisy.
// 
// Daisy is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser Public License as published by
// 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 "frame.h"
#include "librarian.h"
#include "scope_multi.h"
#include "scope_id.h"
#include "units.h"
#include "check.h"
#include "geometry.h"
#include "log.h"
#include "treelog.h"
#include "assertion.h"
#include "mathlib.h"
#include <sstream>

// biopore component.

const char *const Biopore::component = "biopore";

symbol
Biopore::library_id () const
{
  static const symbol id (component);
  return id;
}

symbol
Biopore::x_symbol ()
{
  static const symbol x ("x");
  return x;
}

double
Biopore::max_infiltration_rate (const Geometry& geo, size_t e) const // [cm/h]
{
  const size_t cell = geo.edge_other (e, Geometry::cell_above);
  daisy_assert (cell < geo.cell_size ());

  // Based on Poiseuille equation.
  // Q = pi r^4 rho g L / ( 8 L mu)
  //
  //        Q [cm^3/h]    Infiltration rate
  //  delta P [cm]        Pressure drop.
  //        L [cm]        Length op cylinder.
  //       mu [cm^3/cm/h] Dynamic viscosity.
  //        r [cm]        Biopore radius.
  //        M [cm^-2]     Macropore density.
  //     

  const double rho = 1.0;   // [g/cm^3]
  const double g = 9.81;      // [m/s^2]
  const double g1 = g * 100.0 * 3600.0 * 3600.0;        // [cm/h^2]
  const double mu = (1.0020 / 100.0) * 3600.0;          // [g/cm/h] At 20 dg C. 
  const double M = density (cell);                      // [cm^-2]
  const double r = diameter * 0.5;                      // [cm]
  const double r4 = r * r * r * r;                      // [cm^4]
  const double Q = M_PI * r4 * rho * g1 / (8 * mu);     // [cm^3/h]
  const double q = Q * M;                               // [cm/h]
  
  return q;
}

void 
Biopore::infiltrate (const Geometry& geo, size_t e,
                     const double amount, const double dt)
{
  daisy_assert (std::isfinite (amount));
  const double edge_area = geo.edge_area (e);
  const double total_area = geo.surface_area ();
  const double edge_flux = amount / dt;
  const double total_flux = edge_flux * edge_area / total_area;
  infiltration += total_flux;
}

void 
Biopore::solute_infiltrate (const symbol chem, 
                            const Geometry& geo, const size_t e,
                            const double amount /* [g] */, 
                            const double dt)
{
  daisy_assert (std::isfinite (amount));
  const double total_area = geo.surface_area ();
  solute_infiltration.add_value (chem, solute_infiltration.unit (), 
                                 amount / total_area / dt);
}

double 
Biopore::top_density (const size_t c) const
{
  if (height_start < 0.0)
    return 0.0;
  
  return density (c);
}

void 
Biopore::clear ()
{ 
  infiltration = 0.0; 
  solute_infiltration.clear ();
}

double 
Biopore::matrix_to_biopore (double K_xx, double M_c, double r_c, 
                            double h, double h_3)
{
  const double S = - 4*M_PI*M_c*K_xx*(h-h_3) / 
    (log(M_PI*M_c*r_c*r_c));
  daisy_assert (std::isfinite (S));
  
  if (S < 0.0)
    {
      std::ostringstream tmp;
      tmp << "Bug: S = " << S << " (should be positive), M_c = " << M_c
          << ", r_c = " << r_c << ", K_xx = " << K_xx << ", h = " << h
          << ", h_3 = " << h_3 << ", pi M_c r_c^2 = " << M_PI*M_c*r_c*r_c;
      Assertion::error (tmp.str ());
      return 0.0;
    }
  return S;
}

double 
Biopore::biopore_to_primary (const double K_matrix,
                             const double K_wall_rel, //  K_wall_rel []
                             const double M_c, const double r_c,
                             const double h, const double h_3)
{
  const double r_matrix = std::pow (M_PI * M_c, -0.5);
  const double a1 = 1.1;  // Relative biopore wall size: r_i / r_c [] 
  const double a2 = K_wall_rel; // Relative biopore wall conductivity []
  const double A1 
    = a2 * (std::log (r_matrix / r_c) - std::log (a1))
    / std::log (a1);
  const double B1 = A1 / (1.0 + A1);
  const double K = B1 * K_matrix;
  const double S = -M_c * 2.0 * M_PI * K * (h_3 - h)
    / std::log (r_c / r_matrix);
  if (S < 0.0)
    {
      std::ostringstream tmp;
      tmp << "K_matrix = " << K_matrix << " [cm/h]\n"
          << "K_wall_rel = " << K_wall_rel << " []\n"
          << "M_c = " << M_c << " [cm^-2]\n"
          << "r_c = " << r_c << " [cm]\n"
          << "h = " << h << " [cm]\n"
          << "h_3 = " << h_3 << " [cm]\n"
          << "a1 = " << a1 << " []\n"
          << "a2 = " << a2 << " []\n"
          << "A1 = " << A1 << " []\n"
          << "B1 = " << B1 << " []\n"
          << "K = " << K << " [cm/h]\n"
          << "r_matrix = " << r_matrix << " [cm]\n"
          << "S = " << S << " [] (expected >= 0)";
      daisy_bug (tmp.str ());
    }
  return S;
}

{'fraction_non_alphanumeric': 0.15408209291835814, 'fraction_numerical': 0.02237257555254849, 'mean_word_length': 1.9841184387617765, 'pattern_counts': {'":': 0, '<': 227, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5054882', 'n_tokens_mistral': 4259, 'n_tokens_neox': 4032, 'n_words': 1159}