[Test] Fix solute parameters const correctness and test functions

dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/eff_hydromechanic_dispersion/isotropic.hh

@@ -82,7 +82,7 @@ public:
                 for (int j = 0; j < Traits::dim; ++j) {
                     D[i][j] = 0.;
                     if (Dune::FloatCmp::gt(v_abs,0.)) {
-                        D[i][j] += difference*water_flux.value[i]*water_flux.value[j]/v_abs; // !!!
+                        D[i][j] += difference*water_flux.value[i]*water_flux.value[j]/(v_abs*water_content.value); // !!!
                         if (i==j)
                             D[i][j] += lambda_t.value * v_abs/water_content.value;
                     }
@@ -117,7 +117,7 @@ public:
     std::unique_ptr<EffectiveHydromechanicDispersion<Traits>> clone () const override
     {
         using ThisType = IsotropicEffectiveHydromechanicDispersion<Traits>;
-        return std::make_unique<ThisType>(*this);
+        return std::make_unique<ThisType>(this->get_name(),_long_disp->clone(),_trans_disp->clone());
     }
 };
 

dune/dorie/model/transport/parameterization/hydrodynamic_dispersion/superposition.hh

@@ -84,7 +84,7 @@ public:
     std::unique_ptr<HydrodynamicDispersion<Traits>> clone () const override
     {
         using ThisType = HydrodynamicDispersionSuperposition<Traits>;
-        return std::make_unique<ThisType>(*this);
+        return std::make_unique<ThisType>(this->get_name(),_eff_diff->clone(),_eff_hyrm_disp->clone());
     }
 };
 

dune/dorie/model/transport/parameterization/interface.hh

@@ -113,7 +113,7 @@ public:
     virtual std::unique_ptr<Transport<Traits>> clone () const
     {
         using ThisType = Transport<Traits>;
-        return std::make_unique<ThisType>(*this);
+        return std::make_unique<ThisType>(this->get_name(),_peclet->clone(),_hyrd_disp->clone());
     }
 private:
     const std::string _name;

dune/dorie/model/transport/parameterization/peclet.hh

@@ -84,6 +84,11 @@ public:
         };
     }
 
+    std::unique_ptr<Peclet<Traits>> clone () const
+    {
+        using ThisType = Peclet<Traits>;
+        return std::make_unique<ThisType>(*this);
+    }
 
 private:
     const std::string _name;

dune/dorie/model/transport/solute_parameters.hh

@@ -121,7 +121,7 @@ struct TransportParameterizationFactory
 
       hd_dips = std::make_shared<Parameterization::HydrodynamicDispersionSuperposition<Traits>>(name,eff_diff,eff_hm_disp);
     } else {
-            DUNE_THROW(NotImplemented,"Invalid hydrodyn_disp");
+            DUNE_THROW(NotImplemented,"Invalid hydrodyn_disp" );
     }
     return std::make_shared<Dorie::Parameterization::Transport<Traits>>(name,peclet,hd_dips);
   }

dune/dorie/test/test-param-transport.cc

@@ -21,14 +21,14 @@
 
 #include "test-parameterization.hh"
 
-/// Return parameter values based on the medium name
-/** \param medium Name of the medium
+/// Return parameter values based on the material index
+/** \param index index assigned to the material
  *  \return Map of parameters and associated values
- *  \throw runtime_error Medium name is not known
+ *  \throw runtime_error index is not known
  */
-std::multimap<std::string, double> get_parameter_map (const std::string medium)
+std::multimap<std::string, double> get_parameter_map (const int index)
 {
-    if (medium == "material_0") {
+    if (index == 0) {
         return std::multimap<std::string, double> {
             {"char_length", 1.5E-11},
             {"mol_diff", 2.E-9},
@@ -38,7 +38,7 @@ std::multimap<std::string, double> get_parameter_map (const std::string medium)
             {"hydrodynamic_disp_yx", 0}
         };
     }
-    else if (medium == "material_1") {
+    else if (index == 1) {
         return std::multimap<std::string, double> {
             {"char_length", 1.5E-11},
             {"mol_diff", 2.E-9},
@@ -46,7 +46,7 @@ std::multimap<std::string, double> get_parameter_map (const std::string medium)
             {"alpha", 1.17}
         };
     }
-    else if (medium == "material_2") {
+    else if (index == 2) {
         return std::multimap<std::string, double> {
             {"char_length", 1.5E-11},
             {"mol_diff", 2.E-9},
@@ -57,7 +57,7 @@ std::multimap<std::string, double> get_parameter_map (const std::string medium)
             {"eff_hydromechanic_disp_yx", 0}
         };
     }
-    else if (medium == "material_3") {
+    else if (index == 3) {
         return std::multimap<std::string, double> {
             {"char_length", 1.5E-11},
             {"mol_diff", 2.E-9},
@@ -66,7 +66,7 @@ std::multimap<std::string, double> get_parameter_map (const std::string medium)
             {"lambda_l", 0.025}
         };
     }
-    else if (medium == "material_4") {
+    else if (index == 4) {
         return std::multimap<std::string, double> {
             {"char_length", 1.5E-11},
             {"mol_diff", 2.E-9},
@@ -78,10 +78,21 @@ std::multimap<std::string, double> get_parameter_map (const std::string medium)
         };
     }
     else {
-        throw std::runtime_error("Unknown medium name " + medium);
+        throw std::runtime_error("Unknown index " + std::to_string(index));
     }
 }
 
+template<typename ParamMap, typename ParamGetter>
+bool compare_parameter_maps(const ParamMap& param_map, const ParamGetter& param_getter)
+{
+    bool pass = true;
+    for (const auto& [index,param] : param_map)
+    {
+        pass &= compare_maps(param->parameters(),param_getter(index));
+    }
+    return pass;
+}
+
 /// Test the parameters and parameterization functions
 /** The standard parameters for Sand and Silt are hardcoded here.
  *  We check the actual values and then the return values of the
@@ -93,35 +104,135 @@ template<class SoluteParameters, class Grid>
 void test_new_parameters (const SoluteParameters& sparam,
                                 std::shared_ptr<Grid> grid)
 {
-    // define parameter sets to check against
-    const auto param_mat_0 = get_parameter_map("material_0");
-    const auto param_mat_1 = get_parameter_map("material_1");
-    const auto param_mat_2 = get_parameter_map("material_2");
-    const auto param_mat_3 = get_parameter_map("material_3");
-    const auto param_mat_4 = get_parameter_map("material_4");
-
     // iterate over grid cells and verify parameterization
     auto level_gv = grid->leafGridView();
     for (auto&& cell: elements(level_gv))
     {
+        constexpr int dim = Grid::dimension;
+        using Scalar = Dune::FieldVector<double,1>;
+        using Vector = Dune::FieldVector<double,dim>;
+        using Tensor = Dune::FieldMatrix<double,dim,dim>;
+
         // bind parameterization to cell and fetch parameters
         sparam.bind(cell);
         auto par = sparam.cache();
+        auto name = par->get_name();
         auto parameters = par->parameters();
 
-        // compare values to defined parameter sets
-        assert(compare_maps(param_mat_0, parameters)
-            or compare_maps(param_mat_1, parameters)
-            or compare_maps(param_mat_2, parameters)
-            or compare_maps(param_mat_3, parameters)
-            or compare_maps(param_mat_4, parameters));
+        Scalar   wc(0.5);
+        Vector   flux(0.);
+        flux[dim-1] = -1e-8;
+        flux[0] = 1e-8;
+        auto fluxn = flux.two_norm();
+        // Compute peclet and dispersion with some parameters
+        Tensor disp = sparam.hydrodynamic_dispersion_f()(flux,wc);
+        Scalar pe = sparam.peclet_f()(flux,wc);
+
         namespace FloatCmp = Dune::FloatCmp;
 
-        DUNE_THROW(Dune::NotImplemented,
-            "Test is not finished");
+        // compare values to defined parameter sets
+        if (name == "material_0") { // const hd disp
+            assert(compare_maps(get_parameter_map(0), parameters));
+            
+            for (int i = 0; i < dim; ++i)
+                for (int j = 0; j < dim; ++j)
+                    if (i==j){
+                        assert(FloatCmp::eq(disp[i][j],2.E-8));
+                    }
+                    else
+                        assert(FloatCmp::eq(disp[i][j],0.));
+            auto char_len = parameters.find("char_length")->second;
+            auto mol_diff = parameters.find("mol_diff")->second;
+            assert(FloatCmp::eq(pe,char_len*fluxn/(mol_diff*wc)));
+        } 
+        else if (name == "material_1") // power_law hd disp
+        {
+            assert(compare_maps(get_parameter_map(1), parameters));
+
+            auto alpha = parameters.find("alpha")->second;
+            auto gamma = parameters.find("gamma")->second;
+
+            auto _disp = gamma*std::pow(pe,alpha);
+            
+            for (int i = 0; i < dim; ++i)
+                for (int j = 0; j < dim; ++j)
+                    if (i==j)
+                        assert(FloatCmp::eq(disp[i][j],_disp));
+                    else
+                        assert(FloatCmp::eq(disp[i][j],0.));
+
+            auto char_len = parameters.find("char_length")->second;
+            auto mol_diff = parameters.find("mol_diff")->second;
+            assert(FloatCmp::eq(pe,char_len*fluxn/(mol_diff*wc)));
+        }
+        else if (name == "material_2") // const diff + const hm disp
+        {
+            assert(compare_maps(get_parameter_map(2), parameters));
+            for (int i = 0; i < dim; ++i)
+                for (int j = 0; j < dim; ++j)
+                    if (i==j)
+                        assert(FloatCmp::eq(disp[i][j],2*2.E-8));
+                    else
+                        assert(FloatCmp::eq(disp[i][j],0.));
+
+            auto char_len = parameters.find("char_length")->second;
+            auto mol_diff = parameters.find("mol_diff")->second;
+            assert(FloatCmp::eq(pe,char_len*fluxn/(mol_diff*wc)));
+        } 
+        else if (name == "material_3") // MC1 + iso
+        {
+            assert(compare_maps(get_parameter_map(3), parameters));
 
-        // check evaluation of functions
-        // TODO: ¿¿ HOW ??
+            // MC1
+            auto mol_diff = parameters.find("mol_diff")->second;
+            auto porosity = parameters.find("phi")->second;
+            auto diff = mol_diff*std::pow(wc,7./3.)/std::pow(porosity,2./3.);
+
+            // iso
+            auto lt = parameters.find("lambda_t")->second;
+            auto ll = parameters.find("lambda_l")->second;
+            Tensor _disp(0.);
+            for (int i = 0; i < dim; ++i) {
+                for (int j = 0; j < dim; ++j) {
+                    _disp[i][j] += (ll-lt)*flux[i]*flux[j]/(fluxn*wc);
+                    if (i==j)
+                        _disp[i][j] += lt*fluxn/wc + diff;
+                }
+            }
+
+            for (int i = 0; i < dim; ++i)
+                for (int j = 0; j < dim; ++j)
+                        assert(FloatCmp::eq(disp[i][j],_disp[i][j]));
+
+            auto char_len = parameters.find("char_length")->second;
+            
+            assert(FloatCmp::eq(pe,char_len*fluxn/(mol_diff*wc)));
+        }
+        else if (name == "material_4") // MC2 + const hm disp
+        {
+            assert(compare_maps(get_parameter_map(4), parameters));
+
+            // MC2
+            auto mol_diff = parameters.find("mol_diff")->second;
+            auto porosity = parameters.find("phi")->second;
+            auto diff = mol_diff*wc/std::pow(porosity,2./3.);
+
+            // const hm disp
+            for (int i = 0; i < dim; ++i)
+                for (int j = 0; j < dim; ++j)
+                    if (i==j)
+                        assert(FloatCmp::eq(disp[i][j],double(diff)+2.E-8));
+                    else
+                        assert(FloatCmp::eq(disp[i][j],0.));
+
+            auto char_len = parameters.find("char_length")->second;
+            
+            assert(FloatCmp::eq(pe,char_len*fluxn/(mol_diff*wc)));
+        } 
+        else 
+        {
+            DUNE_THROW(Dune::IOError,"Transport parameterization has an error");
+        }
     }
 }
 
@@ -143,36 +254,29 @@ void test_parameter_manipulation (SoluteParameters& sparam,
     // manipulate via returned map (this is the preferred way)
     auto par = sparam.cache();
     auto parameters = par->parameters();
-    {
+
+    // modify all keys
+    for (auto [key,value] : parameters)
+        value = 0.;
+
+    { // modify specific key
         auto [begin,end] = parameters.equal_range("mol_diff");
         for (auto p = begin; p != end; ++p)
             p->second = 2.5E-9;
-        assert(par->_peclet->_mol_diff.value == 2.5E-9);
     }
-    {
+    { // modify specific key
         auto [begin,end] = parameters.equal_range("char_length");
         for (auto p = begin; p != end; ++p)
             p->second = 1E-11;
         assert(par->_peclet->_char_length.value == 1E-11);
     }
 
-    DUNE_THROW(Dune::NotImplemented,
-            "Test is not finished");
-
-    // manipulate directly via casting (this is not recommended)
-    // using Sup = Dune::Dorie::Parameterization::HydrodynamicDispersionSuperposition<Traits<2>>;
-    // auto& par_sup = dynamic_cast<Sup&>(*par);
-    // auto& par_eff_diff = par_sup._eff_diff;
-    // using MQ2 = Dune::Dorie::Parameterization::MillingtonQuirk2<Traits<2>>;
-    // auto& par_mq2 = dynamic_cast<MQ2&>(*par_eff_diff);
-    // assert(par_mq2._mol_diff.value = 2.5E-9);
-
 }
 
 template<class SoluteParameters, class Grid>
 bool compare_parameters (const SoluteParameters& sparam1,
                          const SoluteParameters& sparam2,
-                               std::shared_ptr<Grid> grid)
+                         const std::shared_ptr<Grid>& grid)
 {
     // iterate over grid cells and verify parameterization
     auto level_gv = grid->leafGridView();
@@ -215,9 +319,7 @@ int main (int argc, char** argv)
         const Dune::Dorie::SoluteParameters<Traits<2>> const_param(inifile,
                                                                    grid,
                                                                    index_map);
-
         Dune::Dorie::SoluteParameters<Traits<2>> param(const_param);
-
         // perform the actual tests
         test_new_parameters(const_param, grid);
         test_new_parameters(param, grid);

dune/dorie/test/test-param-transport.yml

@@ -12,12 +12,12 @@ volumes:
   material_1:
     index: 1
     transport: 
-      type: sthor
+      type: power_law
       parameters:
         char_length: 1.5E-11
         mol_diff: 2.E-9
-        gamma_st: 0.8
-        alpha_st: 1.17
+        gamma: 0.8
+        alpha: 1.17
 
   material_2:
     index: 2