From 0f91a7fa84fea11e6d5b20c8c77203d9ec2b0c76 Mon Sep 17 00:00:00 2001 From: Sylvain Plessis Date: Wed, 2 Apr 2014 22:48:09 -0500 Subject: [PATCH] Derivatives added, need tests now: - eddy diff eval - mol diff eval - diff eval - atm mixture, atmospheric scale height - atm kin eval - physics eval binary is finally useless for the calculations, but it can't hurt to have it tested too --- src/core/include/planet/atmospheric_mixture.h | 31 ++++++ .../include/planet/binary_diffusion.h | 14 +-- .../include/planet/diffusion_evaluator.h | 102 ++++++++++++++++- .../include/planet/eddy_diffusion_evaluator.h | 13 +++ .../planet/molecular_diffusion_evaluator.h | 104 +++++++++++++++++- .../include/planet/planet_physics_evaluator.h | 52 ++++++++- .../include/planet/atmospheric_kinetics.h | 35 +++++- 7 files changed, 335 insertions(+), 16 deletions(-) diff --git a/src/core/include/planet/atmospheric_mixture.h b/src/core/include/planet/atmospheric_mixture.h index b41c40d..e59fbbc 100644 --- a/src/core/include/planet/atmospheric_mixture.h +++ b/src/core/include/planet/atmospheric_mixture.h @@ -189,6 +189,10 @@ namespace Planet template const CoeffType atmospheric_scale_height(const VectorStateType &molar_densities,const StateType &z) const; + //! + template + 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; @@ -323,6 +327,33 @@ namespace Planet return (this->H(Mm,_temperature.neutral_temperature(z),z)); } + template + template + inline + void AtmosphericMixture::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 inline const CoeffType AtmosphericMixture::total_bottom_density() const diff --git a/src/diffusion/include/planet/binary_diffusion.h b/src/diffusion/include/planet/binary_diffusion.h index 5469049..f538c73 100644 --- a/src/diffusion/include/planet/binary_diffusion.h +++ b/src/diffusion/include/planet/binary_diffusion.h @@ -98,12 +98,12 @@ class BinaryDiffusion{ //! template 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 - 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; //! @@ -274,12 +274,12 @@ void BinaryDiffusion::set_binary_diffusion(const Antioch::Species &mo template template inline -void BinaryDiffusion::binary_coefficient_and_derivatives(const StateType &T, const StateType &P, - StateType &Dij, StateType &Dij_dT, StateType &Dij_dP) const +void BinaryDiffusion::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; } diff --git a/src/diffusion/include/planet/diffusion_evaluator.h b/src/diffusion/include/planet/diffusion_evaluator.h index c406afd..ed354a0 100644 --- a/src/diffusion/include/planet/diffusion_evaluator.h +++ b/src/diffusion/include/planet/diffusion_evaluator.h @@ -69,6 +69,15 @@ namespace Planet{ const VectorStateType &dmolar_concentrations_dz, const StateType &z, VectorStateType &omegas) const; + //! + template + 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 @@ -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 * [ @@ -150,6 +159,97 @@ namespace Planet{ return; } + template + template + inline + void DiffusionEvaluator::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; + + } + } diff --git a/src/diffusion/include/planet/eddy_diffusion_evaluator.h b/src/diffusion/include/planet/eddy_diffusion_evaluator.h index e9e44e5..32cd0c6 100644 --- a/src/diffusion/include/planet/eddy_diffusion_evaluator.h +++ b/src/diffusion/include/planet/eddy_diffusion_evaluator.h @@ -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 + 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 + 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 &mix, const CoeffType &K0 = -1.); diff --git a/src/diffusion/include/planet/molecular_diffusion_evaluator.h b/src/diffusion/include/planet/molecular_diffusion_evaluator.h index c527370..459edbd 100644 --- a/src/diffusion/include/planet/molecular_diffusion_evaluator.h +++ b/src/diffusion/include/planet/molecular_diffusion_evaluator.h @@ -102,6 +102,14 @@ namespace Planet void set_medium_species(const std::vector &medium_species); + template + void Dtilde_and_derivs_dn(const VectorStateType &molar_concentrations, const StateType &T, const StateType &nTot, + VectorStateType &Dtilde, MatrixStateType &dD_dns) const; + + template + 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 @@ -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]; @@ -199,6 +208,99 @@ namespace Planet return; } + template + template + inline + void MolecularDiffusionEvaluator::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 + template + inline + void MolecularDiffusionEvaluator::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() * 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 diff --git a/src/grins_interface/include/planet/planet_physics_evaluator.h b/src/grins_interface/include/planet/planet_physics_evaluator.h index a7855a0..b30f1e5 100644 --- a/src/grins_interface/include/planet/planet_physics_evaluator.h +++ b/src/grins_interface/include/planet/planet_physics_evaluator.h @@ -45,6 +45,15 @@ namespace Planet libMesh::Real chemical_term(unsigned int s) const; + //! domega_s_dn_i = A_TERM + B_TERM * d(dn_s_dz)_dn_i + libMesh::Real ddiffusion_term_s_d_n_i_A_TERM(unsigned int s, unsigned int i) const; + + //! domega_s_dn_i = A_TERM + B_TERM * d(dn_s_dz)_dn_i + libMesh::Real ddiffusion_term_s_d_n_i_B_TERM(unsigned int s, unsigned int i) const; + + //! domega_dot_s_dn_i + libMesh::Real dchemical_term_dn_i(unsigned int s, unsigned int i) const; + const EddyDiffusionEvaluator& eddy_diffusion() const; const MolecularDiffusionEvaluator& molecular_diffusion() const; @@ -87,6 +96,10 @@ namespace Planet VectorCoeffType _omegas; VectorCoeffType _omegas_dots; + MatrixCoeffType _domegas_dn_A_TERM; + MatrixCoeffType _domegas_dn_B_TERM; + MatrixCoeffType _domegas_dots_dn; + MatrixCoeffType _cache_composition; VectorCoeffType _cache_altitudes; std::map _cache; @@ -113,8 +126,18 @@ namespace Planet _diffusion(_molecular_diffusion,_eddy_diffusion,_composition,helper.temperature()), _scaling_factor(helper.scaling_factor()) { - _omegas.resize(_kinetics.neutral_kinetics().reaction_set().n_species()); - _omegas_dots.resize(_kinetics.neutral_kinetics().reaction_set().n_species()); + _omegas.resize(_kinetics.neutral_kinetics().n_species()); + _omegas_dots.resize(_kinetics.neutral_kinetics().n_species()); + + _domegas_dots_dn.resize(_kinetics.neutral_kinetics().n_species()); + _domegas_dn_A_TERM.resize(_kinetics.neutral_kinetics().n_species()); + _domegas_dn_B_TERM.resize(_kinetics.neutral_kinetics().n_species()); + for(unsigned int s = 0; s < _kinetics.neutral_kinetics().n_species(); s++) + { + _domegas_dots_dn[s].resize(_kinetics.neutral_kinetics().n_species()); + _domegas_dn_A_TERM[s].resize(_kinetics.neutral_kinetics().n_species()); + _domegas_dn_B_TERM[s].resize(_kinetics.neutral_kinetics().n_species()); + } _photon.set_photon_flux_at_top(helper.lambda_hv(), helper.phy1AU(), Constants::Saturn::d_Sun()); @@ -147,9 +170,8 @@ namespace Planet molar[i] = molar_concentrations[i] * _scaling_factor; dmolar[i] = dmolar_concentrations_dz[i] * _scaling_factor; } - _diffusion.diffusion(molar,dmolar,z,_omegas); - _kinetics.chemical_rate(molar,this->get_cache(z),z,_omegas_dots); - + _diffusion.diffusion_and_derivs(molar,dmolar,z,_omegas,_domegas_dn_A_TERM,_domegas_dn_B_TERM); + _kinetics.chemical_rate_and_derivs(molar,this->get_cache(z),z,_omegas_dots,_domegas_dots_dn); this->update_cache(molar,z); @@ -228,6 +250,26 @@ namespace Planet return _omegas_dots[s] / _scaling_factor; } + // this is dw_s / dn_i + template + libMesh::Real PlanetPhysicsEvaluator::ddiffusion_term_s_d_n_i_A_TERM(unsigned int s, unsigned int i) const + { + return _domegas_dn_A_TERM[s][i] / _scaling_factor; + } + + // this is factor such that factor * d(dn_s / dz) / dn_i + template + libMesh::Real PlanetPhysicsEvaluator::ddiffusion_term_s_d_n_i_B_TERM(unsigned int s, unsigned int i) const + { + return _domegas_dn_B_TERM[s][i] / _scaling_factor; + } + + template + libMesh::Real PlanetPhysicsEvaluator::dchemical_term_dn_i(unsigned int s, unsigned int i) const + { + return _domegas_dots_dn[s][i] / _scaling_factor; + } + template const EddyDiffusionEvaluator & PlanetPhysicsEvaluator::eddy_diffusion() const { diff --git a/src/kinetics/include/planet/atmospheric_kinetics.h b/src/kinetics/include/planet/atmospheric_kinetics.h index 1c6de32..88b5826 100644 --- a/src/kinetics/include/planet/atmospheric_kinetics.h +++ b/src/kinetics/include/planet/atmospheric_kinetics.h @@ -77,6 +77,10 @@ namespace Planet void chemical_rate(const VectorStateType &molar_concentrations, const VectorStateType &sum_concentrations, const StateType &z, VectorStateType &kin_rates) const; + template + void chemical_rate_and_derivs(const VectorStateType &molar_concentrations, const VectorStateType &sum_concentrations, + const StateType &z, VectorStateType &kin_rates, MatrixStateType &dkin_rates_dn) const; + //! Newton solver for the ionic system template void add_ionic_contribution(const VectorStateType &molar_concentrations, const StateType &z, VectorStateType &kin_rates) const; @@ -139,8 +143,8 @@ namespace Planet const StateType &z, VectorStateType &kin_rates) const { - kin_rates.resize(_composition.neutral_composition().n_species(),0.L); - VectorCoeffType dummy; + antioch_assert_equal_to(kin_rates.size(),_composition.neutral_composition().n_species()); + VectorStateType dummy; dummy.resize(_composition.neutral_composition().n_species(),0.L); //everything is irreversible _photon.update_photon_flux(molar_concentrations, sum_concentrations, z); _neutral_reactions.compute_mole_sources(_temperature.neutral_temperature(z), @@ -151,6 +155,33 @@ namespace Planet return; } + template + template + inline + void AtmosphericKinetics::chemical_rate_and_derivs(const VectorStateType &molar_concentrations, const VectorStateType &sum_concentrations, + const StateType &z, VectorStateType &kin_rates, MatrixStateType &dkin_rates_dn) const + { + antioch_assert_equal_to(kin_rates.size(),_composition.neutral_composition().n_species()); + antioch_assert_equal_to(dkin_rates_dn.size(),_composition.neutral_composition().n_species()); +#ifdef NDEBUG +#else + for(unsigned int s = 0; s < _composition.neutral_composition().n_species(); s++) + { + antioch_assert_equal_to(dkin_rates_dn[s].size(),_composition.neutral_composition().n_species()); + } +#endif + VectorStateType dummy; + VectorStateType ddummy_dT; + VectorStateType dkin_dT; + dummy.resize(_composition.neutral_composition().n_species(),0.L); //everything is irreversible + dkin_dT.resize(_composition.neutral_composition().n_species(),0.L); //no temp + _photon.update_photon_flux(molar_concentrations, sum_concentrations, z); + _neutral_reactions.compute_mole_sources_and_derivs(_temperature.neutral_temperature(z),molar_concentrations,dummy,ddummy_dT, + kin_rates,dkin_dT,dkin_rates_dn); + + this->add_ionic_contribution(molar_concentrations,z,kin_rates); + } + template template inline