Skip to content

Commit

Permalink
Merge pull request Planet-INC#3 from SylvainPlessis/planet_grins
Browse files Browse the repository at this point in the history 
Derivatives added, need tests now:
  • Loading branch information
@pbauman
pbauman committed Apr 3, 2014
2 parents c82327a + 0f91a7f commit f3f93cd
Show file tree
Hide file tree
Showing 7 changed files with 335 additions and 16 deletions.
31 changes: 31 additions & 0 deletions src/core/include/planet/atmospheric_mixture.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -189,6 +189,10 @@ namespace Planet
template<typename StateType, typename VectorStateType>
const CoeffType atmospheric_scale_height(const VectorStateType &molar_densities,const StateType &z) const;

//!
template<typename StateType, typename VectorStateType>
void datmospheric_scale_height_dn_i(const VectorStateType &molar_densities, const StateType &z, StateType & Ha, VectorStateType & dHa_dn_i) const;

//!
const CoeffType total_bottom_density() const;

Expand Down Expand Up @@ -323,6 +327,33 @@ namespace Planet
return (this->H(Mm,_temperature.neutral_temperature(z),z));
}

template<typename CoeffType, typename VectorCoeffType, typename MatrixCoeffType>
template<typename StateType, typename VectorStateType>
inline
void AtmosphericMixture<CoeffType,VectorCoeffType,MatrixCoeffType>::datmospheric_scale_height_dn_i(const VectorStateType &molar_densities,
const StateType &z,
StateType & Ha, VectorStateType & dHa_dn_i) const
{
antioch_assert_equal_to(molar_densities.size(),_neutral_composition.n_species());

CoeffType Mm;
Antioch::set_zero(Mm);
CoeffType nTot;
Antioch::set_zero(nTot);
for(unsigned int s = 0; s < _neutral_composition.n_species(); s++)
{
Mm += molar_densities[s] * _neutral_composition.M(s);
nTot += molar_densities[s];
}
Mm /= nTot;

Ha = this->H(Mm,_temperature.neutral_temperature(z),z);
for(unsigned int s = 0; s < _neutral_composition.n_species(); s++)
{
dHa_dn_i[s] = (_neutral_composition.M(s) - Ha) / nTot;
}
}

template<typename CoeffType, typename VectorCoeffType, typename MatrixCoeffType>
inline
const CoeffType AtmosphericMixture<CoeffType,VectorCoeffType,MatrixCoeffType>::total_bottom_density() const
Expand Down
14 changes: 7 additions & 7 deletions src/diffusion/include/planet/binary_diffusion.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,12 +98,12 @@ class BinaryDiffusion{
//!
template<typename StateType>
ANTIOCH_AUTO(StateType)
binary_coefficient_deriv_P(const StateType &T, const StateType &P) const
ANTIOCH_AUTOFUNC(StateType,this->binary_coefficient(T,P) / P)
binary_coefficient_deriv_n(const StateType &T, const StateType &P, const StateType &nTot) const
ANTIOCH_AUTOFUNC(StateType, - this->binary_coefficient(T,P) / nTot)

//!
template<typename StateType>
void binary_coefficient_and_derivatives(const StateType &T, const StateType &P,
void binary_coefficient_and_derivatives(const StateType &T, const StateType &P, const StateType &nTot,
StateType &Dij, StateType &Dij_dT, StateType &Dij_dP) const;

//!
Expand Down Expand Up @@ -274,12 +274,12 @@ void BinaryDiffusion<CoeffType>::set_binary_diffusion(const Antioch::Species &mo
template<typename CoeffType>
template<typename StateType>
inline
void BinaryDiffusion<CoeffType>::binary_coefficient_and_derivatives(const StateType &T, const StateType &P,
StateType &Dij, StateType &Dij_dT, StateType &Dij_dP) const
void BinaryDiffusion<CoeffType>::binary_coefficient_and_derivatives(const StateType &T, const StateType &P, const StateType &nTot,
StateType &Dij, StateType &Dij_dT, StateType &Dij_dns) const
{
Dij = this->binary_coefficient(T,P);
Dij_dT = this->binary_coefficient_deriv_T(T,P);
Dij_dP = this->binary_coefficient_deriv_P(T,P);
Dij_dT = this->binary_coefficient_deriv_T(T,P);
Dij_dns = this->binary_coefficient_deriv_P(T,P,nTot);
return;
}

Expand Down
102 changes: 101 additions & 1 deletion src/diffusion/include/planet/diffusion_evaluator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,15 @@ namespace Planet{
const VectorStateType &dmolar_concentrations_dz,
const StateType &z, VectorStateType &omegas) const;

//!
template<typename StateType, typename VectorStateType, typename MatrixStateType>
void diffusion_and_derivs(const VectorStateType &molar_concentrations,
const VectorStateType &dmolar_concentrations_dz,
const StateType &z,
VectorStateType &omegas,
MatrixStateType &domegas_dn_i_A_TERM,
MatrixStateType &domegas_dn_i_B_TERM) const;

};

template <typename CoeffType, typename VectorCoeffType, typename MatrixCoeffType>
Expand Down Expand Up @@ -129,7 +138,7 @@ namespace Planet{
// eddy diff
CoeffType eddy_K = _eddy_diffusion.K(nTot);

// in km-2.s-1
// in cm-3.km.s-1
for(unsigned int s = 0; s < _mixture.neutral_composition().n_species(); s++)
{
omegas[s] = Antioch::constant_clone(T,1e-10) * (//omega = - ns * Dtilde * [
Expand All @@ -150,6 +159,97 @@ namespace Planet{
return;
}

template <typename CoeffType, typename VectorCoeffType, typename MatrixCoeffType>
template <typename StateType, typename VectorStateType, typename MatrixStateType>
inline
void DiffusionEvaluator<CoeffType, VectorCoeffType,MatrixCoeffType>::diffusion_and_derivs(const VectorStateType &molar_concentrations,
const VectorStateType &dmolar_concentrations_dz,
const StateType &z,
VectorStateType &omegas,
MatrixStateType &domegas_dn_i_A_TERM,
MatrixStateType &domegas_dn_i_B_TERM) const
{

antioch_assert_equal_to(_mixture.neutral_composition().n_species(),molar_concentrations.size());
antioch_assert_equal_to(_mixture.neutral_composition().n_species(),dmolar_concentrations_dz.size());
//params
StateType nTot;
Antioch::set_zero(nTot);
for(unsigned int s = 0; s < molar_concentrations.size();)
{
nTot += molar_concentrations[s];
}
StateType T = _temperature.neutral_temperature(z);
StateType dT_dz = _temperature.dneutral_temperature_dz(z);
StateType dT_dz_T = dT_dz / T;

//eddy
CoeffType dK_dn = _eddy_diffusion.K_deriv_ns(nTot);
CoeffType K = _eddy_diffusion.K(nTot);

//molecular
VectorStateType Dtilde;
Dtilde.resize(_mixture.neutral_composition().n_species(),0.);
MatrixCoeffType dDtilde_dn;
dDtilde_dn.resize(_mixture.neutral_composition().n_species());
for(unsigned int s = 0; s < _mixture.neutral_composition().n_species(); s++)
{
dDtilde_dn[s].resize(_mixture.neutral_composition().n_species());
}

_molecular_diffusion.Dtilde_and_derivs_dn(molar_concentrations,T,nTot,Dtilde,dDtilde_dn);

//scale heights
CoeffType Ha;
VectorCoeffType dHa_dn_i;
VectorCoeffType Hs;
Hs.resize(_mixture.neutral_composition().n_species(),0.);
dHa_dn_i.resize(_mixture.neutral_composition().n_species(),0.);

_mixture.datmospheric_scale_height_dn_i(molar_concentrations,z,Ha,dHa_dn_i);
_mixture.scale_heights(z,Hs);

for(unsigned int s = 0; s < _mixture.neutral_composition().n_species(); s++)
{
omegas[s] = //omega = - ns * Dtilde * [
- Dtilde[s] *
(
dmolar_concentrations_dz[s] // 1/ns * dns_dz
+ molar_concentrations[s]/Hs[s] // + 1/Hs
+ molar_concentrations[s] * dT_dz_T // + 1/T * dT_dz * (
* (Antioch::constant_clone(T,1.) + ((nTot - molar_concentrations[s])/nTot) * _mixture.thermal_coefficient()[s]) //1 + (1 - xs)*alphas ) ]
)
- K * // - ns * K * (
(
dmolar_concentrations_dz[s] // 1/ns * dns_dz
+ molar_concentrations[s]/Ha // + 1/Ha
+ molar_concentrations[s] * dT_dz_T //+1/T * dT_dz )
);
domegas_dn_i_B_TERM[s].clear();
domegas_dn_i_B_TERM[s].resize(_mixture.neutral_composition().n_species(),- (Dtilde[s] + K) * Antioch::constant_clone(T,1e-10)); //to cm-3.km.s-1
for(unsigned int i = 0; i < _mixture.neutral_composition().n_species(); i++)
{
domegas_dn_i_A_TERM[s][i] = - dmolar_concentrations_dz[s] * (dDtilde_dn[s][i] + dK_dn)
- omegas[s] / molar_concentrations[s]
- molar_concentrations[s] * (
dDtilde_dn[s][i] * ( Antioch::constant_clone(Hs[s],1.)/Hs[s]
+ dT_dz_T * (Antioch::constant_clone(T,1.) + ((nTot - molar_concentrations[s])/nTot) * _mixture.thermal_coefficient()[s])
)
+ dK_dn * (Antioch::constant_clone(Ha,1) / Ha + dT_dz / T)
+ Dtilde[s] * dT_dz_T * _mixture.thermal_coefficient()[s] /nTot * (Antioch::constant_clone(nTot,1) - molar_concentrations[s] / nTot)
- K * dHa_dn_i[i] / (Ha * Ha)
);
// in cm-3.km.s-1
domegas_dn_i_A_TERM[s][i] *= Antioch::constant_clone(T,1e-10);
}
// in cm-3.km.s-1
omegas[s] *= Antioch::constant_clone(T,1e-10);
}

return;

}


}

Expand Down
13 changes: 13 additions & 0 deletions src/diffusion/include/planet/eddy_diffusion_evaluator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,19 @@ namespace Planet{
K(const StateType &ntot) const
ANTIOCH_AUTOFUNC(StateType,_K0 * Antioch::ant_sqrt(Antioch::ant_abs(_mixture.total_bottom_density()/ntot)))

//! \return deddy_dT coefficient in cm2.s-1.K-1
template<typename StateType>
ANTIOCH_AUTO(StateType)
K_deriv_T(const StateType &T) const
ANTIOCH_AUTOFUNC(StateType,Antioch::zero_clone(T))

//! \return deddy_dn coefficient in cm2.s-1.K-1.(cm-3)-1
template<typename StateType>
ANTIOCH_AUTO(StateType)
K_deriv_ns(const StateType &ntot) const
ANTIOCH_AUTOFUNC(StateType,- this->K(ntot) / (ntot * Antioch::constant_clone(ntot,2)))


//!
EddyDiffusionEvaluator(const AtmosphericMixture<CoeffType,VectorCoeffType,MatrixCoeffType> &mix,
const CoeffType &K0 = -1.);
Expand Down
104 changes: 103 additions & 1 deletion src/diffusion/include/planet/molecular_diffusion_evaluator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,14 @@ namespace Planet

void set_medium_species(const std::vector<std::string> &medium_species);

template<typename StateType, typename VectorStateType, typename MatrixStateType>
void Dtilde_and_derivs_dn(const VectorStateType &molar_concentrations, const StateType &T, const StateType &nTot,
VectorStateType &Dtilde, MatrixStateType &dD_dns) const;

template<typename StateType, typename VectorStateType>
void Dtilde_s_dn_i(unsigned int s, const VectorStateType &molar_concentrations, const StateType &T,const StateType &nTot,
StateType & Dtilde_s, VectorStateType &dDtilde_s_dn_i) const;

};

template<typename CoeffType, typename VectorCoeffType, typename MatrixCoeffType>
Expand Down Expand Up @@ -161,7 +169,8 @@ namespace Planet
antioch_assert_equal_to(molar_concentrations.size(),_mixture.neutral_composition().n_species());

Dtilde.resize(molar_concentrations.size(),0.L);
CoeffType nTot(0.L);
CoeffType nTot;
Antioch::set_zero(nTot);
for(unsigned int s = 0; s < molar_concentrations.size(); s++)
{
nTot += molar_concentrations[s];
Expand Down Expand Up @@ -199,6 +208,99 @@ namespace Planet
return;
}

template<typename CoeffType, typename VectorCoeffType, typename MatrixCoeffType>
template<typename StateType, typename VectorStateType, typename MatrixStateType>
inline
void MolecularDiffusionEvaluator<CoeffType, VectorCoeffType,MatrixCoeffType>::Dtilde_and_derivs_dn(const VectorStateType &molar_concentrations, const StateType &T, const StateType &nTot,
VectorStateType &Dtilde, MatrixStateType &dD_dns) const
{
antioch_assert_equal_to(molar_concentrations.size(),_mixture.neutral_composition().n_species());
antioch_assert_equal_to(Dtilde.size(),_mixture.neutral_composition().n_species());
antioch_assert_equal_to(dD_ns.size(),_mixture.neutral_composition().n_species());
#ifdef NDEBUG
#else
for(unsigned int s = 0; s < dD_dns.size(); s++)
{
antioch_assert_equal_to(dD_dns[s].size(),_mixture.neutral_composition().n_species());
}
#endif

for(unsigned int s = 0; s < dD_dns.size(); s++)
{
this->Dtilde_s_dn_i(s,molar_concentrations,T,nTot,Dtilde[s],dD_dns[s]);
}
}

template<typename CoeffType, typename VectorCoeffType, typename MatrixCoeffType>
template<typename StateType, typename VectorStateType>
inline
void MolecularDiffusionEvaluator<CoeffType, VectorCoeffType,MatrixCoeffType>::Dtilde_s_dn_i(unsigned int s,
const VectorStateType &molar_concentrations,
const StateType &T,const StateType &nTot,
StateType & Dtilde_s, VectorStateType &dDtilde_s_dn) const
{

antioch_assert_equal_to(molar_concentrations.size(),_mixture.neutral_composition().n_species());
antioch_assert_equal_to(dDtilde_s_dn_i.size(),_mixture.neutral_composition().n_species());

//D_s = (nT - ns) / (sum_{medium} n_{medium}/D_{medium,s})
CoeffType Ds = (nTot - molar_concentrations[s]); //cm-3
//Ds denominator : sum_{j_m} n_{j_m}/D_{s,j_m}
CoeffType n_D;
Antioch::set_zero(n_D);
CoeffType p = nTot * Antioch::constant_clone(nTot,1e6) //cm-3 -> m-3
* Constants::Universal::kb<CoeffType>() * T;

for(unsigned int m = 0; m < _n_medium; m++)
{
if(_i_medium[m] == s)continue;
n_D += molar_concentrations[_i_medium[m]] / this->binary_coefficient(m,s,T,p);
}
Ds /= n_D;

CoeffType meanM;
Antioch::set_zero(meanM);
CoeffType ntot_s = nTot - molar_concentrations[s]; //ntot - ns
for(unsigned int j = 0; j < _mixture.neutral_composition().n_species(); j++)
{
if(j == s)continue;
meanM += _mixture.neutral_composition().M(j) * molar_concentrations[j] / ntot_s; //x_i without s: ni/(ntot - ns)
}

CoeffType sum = (Antioch::constant_clone(nTot,1) - molar_concentrations[s]/nTot *
(Antioch::constant_clone(nTot,1) - _mixture.neutral_composition().M(s)/meanM));

Dtilde_s = Ds / sum;

for(unsigned int i = 0; i < molar_concentrations.size(); i++)
{

CoeffType first = (i == s)?Antioch::constant_clone(nTot,1.):Antioch::zero_clone(nTot);

CoeffType second;
Antioch::set_zero(second);
for(unsigned int m = 0; m < _n_medium; m++)
{
if(_i_medium[m] == s)continue;
if(_i_medium[m] == i)second = Antioch::constant_clone(nTot,1)/ this->binary_coefficient(m,s,T,p);
}

CoeffType dDs_dn_i = Ds * ( first - (second + n_D/nTot)/n_D);

CoeffType third = _mixture.neutral_composition().M(i) * molar_concentrations[i];
if(i == s)Antioch::set_zero(third);
CoeffType dmeanM_dn_i = (third - meanM * first) / (nTot - molar_concentrations[s]);

CoeffType dsum_dn_i = molar_concentrations[s]/(nTot * nTot) + _mixture.neutral_composition().M(s)/(meanM * meanM) * dmeanM_dn_i;
if(i == s)dsum_dn_i -= Antioch::constant_clone(nTot,1) / nTot;


dDtilde_s_dn[i] = Dtilde_s * (dDs_dn_i/Ds - dsum_dn_i/sum);
}

return;
}

}

#endif
Loading

0 comments on commit f3f93cd

Please sign in to comment.