From fe40ec1293cfdedba7eed20f8ae56e19deebff65 Mon Sep 17 00:00:00 2001 From: Sylvain Plessis Date: Mon, 9 Mar 2015 13:06:50 -0500 Subject: [PATCH 1/4] EG lib correction for small molar fractions --- .../molecular_binary_diffusion_utils.h | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h b/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h index aaf0ee03..64125c80 100644 --- a/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h +++ b/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h @@ -191,9 +191,28 @@ namespace Antioch VectorStateType molar_fractions = zero_clone(mass_fractions); mixture.X(mixture.M(mass_fractions),mass_fractions,molar_fractions); +// EGlib traces management, see doc: http://blanche.polytechnique.fr/www.eglib/manual.ps +// page 5 + typename value_type::type eps(1e-16); + typename value_type::type mol_frac_sum = zero_clone(mass_fractions[0]); + for(unsigned int s = 0; s < molar_fractions.size(); s++) + { + mol_frac_sum += molar_fractions[s]; + } + mol_frac_sum /= (typename rebind< unsigned int, typename value_type::type >::type)(molar_fractions.size()); + + typename value_type::type M_tr = zero_clone(mass_fractions[0]); + for(unsigned int s = 0; s < molar_fractions.size(); s++) + { + molar_fractions[s] += eps * (mol_frac_sum - molar_fractions[s]); // add perturbation + M_tr += molar_fractions[s] * mixture.M(s); + } + + +// mass_fraction = molar_fraction * molar_mass / perturbed_molar_mass for(unsigned int s = 0; s < ds.size(); s++) { - ds[s] = constant_clone(mass_fractions[s],1) - mass_fractions[s]; + ds[s] = constant_clone(mass_fractions[s],1) - mixture.M(s) / M_tr * molar_fractions[s]; typename value_type::type denom = zero_clone(mass_fractions[0]); for(unsigned int j = 0; j < ds.size(); j++) { From faff8937cab19c4aa33fda8862d6be5bc018b6f3 Mon Sep 17 00:00:00 2001 From: Sylvain Plessis Date: Tue, 10 Mar 2015 12:17:34 -0500 Subject: [PATCH 2/4] verbosity added mostly, a very small change to simplify (StateType is defined) --- .../molecular_binary_diffusion_utils.h | 25 +++++++++++++------ 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h b/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h index 64125c80..9dda3ebf 100644 --- a/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h +++ b/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h @@ -26,10 +26,14 @@ #ifndef ANTIOCH_MOLECULAR_BINARY_DIFFUSION_UTILS_H #define ANTIOCH_MOLECULAR_BINARY_DIFFUSION_UTILS_H +// Antioch #include "antioch/transport_species.h" #include "antioch/molecular_binary_diffusion_utils_decl.h" #include "antioch/molecular_binary_diffusion_building.h" +// C++ +#include + namespace Antioch { // getting tag @@ -188,20 +192,26 @@ namespace Antioch antioch_assert_equal_to(ds.size(),mixture.n_species()); antioch_assert_equal_to(ds.size(),mass_fractions.size()); + // convenient + typedef typename value_type::type StateType; + VectorStateType molar_fractions = zero_clone(mass_fractions); mixture.X(mixture.M(mass_fractions),mass_fractions,molar_fractions); // EGlib traces management, see doc: http://blanche.polytechnique.fr/www.eglib/manual.ps // page 5 - typename value_type::type eps(1e-16); - typename value_type::type mol_frac_sum = zero_clone(mass_fractions[0]); +// EGlib uses eps = 1e-16 + typename raw_type::type eps(std::numeric_limits::epsilon() * 10); + StateType mol_frac_sum = zero_clone(mass_fractions[0]); for(unsigned int s = 0; s < molar_fractions.size(); s++) { mol_frac_sum += molar_fractions[s]; } - mol_frac_sum /= (typename rebind< unsigned int, typename value_type::type >::type)(molar_fractions.size()); + mol_frac_sum /= (typename rebind< unsigned int, StateType >::type)(molar_fractions.size()); - typename value_type::type M_tr = zero_clone(mass_fractions[0]); +// (i) evaluate perturbed mole fractions +// (ii) evaluate the perturbed mean molar weight ... + StateType M_tr = zero_clone(mass_fractions[0]); for(unsigned int s = 0; s < molar_fractions.size(); s++) { molar_fractions[s] += eps * (mol_frac_sum - molar_fractions[s]); // add perturbation @@ -209,11 +219,12 @@ namespace Antioch } -// mass_fraction = molar_fraction * molar_mass / perturbed_molar_mass +// (ii) .. evaluate perturbed mass_fraction [= molar_fraction * molar_mass / perturbed_molar_mass] +// (iii) use perturbed values to evaluate the transport properties [ (1 - y_s) / sum_i x_i / D_{is} ] for(unsigned int s = 0; s < ds.size(); s++) { - ds[s] = constant_clone(mass_fractions[s],1) - mixture.M(s) / M_tr * molar_fractions[s]; - typename value_type::type denom = zero_clone(mass_fractions[0]); + ds[s] = StateType(1) - mixture.M(s) / M_tr * molar_fractions[s]; + StateType denom = zero_clone(mass_fractions[0]); for(unsigned int j = 0; j < ds.size(); j++) { if(j == s)continue; From d8d10d789e51ea7f5f83e2ccacf79fcb70fe535c Mon Sep 17 00:00:00 2001 From: Sylvain Plessis Date: Wed, 11 Mar 2015 19:12:04 -0500 Subject: [PATCH 3/4] typo fix --- .../include/antioch/molecular_binary_diffusion_utils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h b/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h index 9dda3ebf..3be53ac4 100644 --- a/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h +++ b/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h @@ -201,7 +201,7 @@ namespace Antioch // EGlib traces management, see doc: http://blanche.polytechnique.fr/www.eglib/manual.ps // page 5 // EGlib uses eps = 1e-16 - typename raw_type::type eps(std::numeric_limits::epsilon() * 10); + typename raw_value_type::type eps(std::numeric_limits::epsilon() * 10); StateType mol_frac_sum = zero_clone(mass_fractions[0]); for(unsigned int s = 0; s < molar_fractions.size(); s++) { From c111960be0033c2e206bc4e1398873fea2aaf04e Mon Sep 17 00:00:00 2001 From: Sylvain Plessis Date: Wed, 25 Mar 2015 11:48:06 -0500 Subject: [PATCH 4/4] StateType can be vectorized and might require rules to be constant-copied --- .../include/antioch/molecular_binary_diffusion_utils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h b/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h index 3be53ac4..3c088b6d 100644 --- a/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h +++ b/src/diffusion/include/antioch/molecular_binary_diffusion_utils.h @@ -223,7 +223,7 @@ namespace Antioch // (iii) use perturbed values to evaluate the transport properties [ (1 - y_s) / sum_i x_i / D_{is} ] for(unsigned int s = 0; s < ds.size(); s++) { - ds[s] = StateType(1) - mixture.M(s) / M_tr * molar_fractions[s]; + ds[s] = constant_clone(mass_fractions[0],1) - mixture.M(s) / M_tr * molar_fractions[s]; StateType denom = zero_clone(mass_fractions[0]); for(unsigned int j = 0; j < ds.size(); j++) {