Skip to content

Commit ed7c465

Browse files
bangerthbangerth
committed
Restructure code and check inputs.
1 parent b7d12ee commit ed7c465

File tree

2 files changed

+108
-40
lines changed

2 files changed

+108
-40
lines changed

‎include/aspect/boundary_temperature/dynamic_core.h‎

Lines changed: 14 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -467,40 +467,46 @@ namespace aspect
467467
* (regarding the core cooling rated Tc/dt) for a given core-mantle boundary (CMB)
468468
* temperature @p Tc
469469
*/
470-
void get_specific_heating(const double Tc, double &Qs, double &Es) const;
470+
std::pair<double,double>
471+
get_specific_heating(const double Tc) const;
471472

472473
/**
473474
* Calculate energy (@p Qr) and entropy (@p Er) change rate factor (regarding the
474475
* radioactive heating rate H) for a given CMB temperature @p Tc
475476
*/
476-
void get_radio_heating(const double Tc, double &Qr, double &Er) const;
477+
std::pair<double,double>
478+
get_radio_heating(const double Tc) const;
477479

478480
/**
479481
* Calculate energy (@p Qg) and entropy (@p Eg) change rate factor
480482
* (regarding the inner core growth rate dR/dt) for a given
481483
* @p Tc (CMB temperature), @p r (inner core radius), and @p X
482484
* (light element concentration)
483485
*/
484-
void get_gravity_heating(const double Tc, const double r, const double X, double &Qg, double &Eg) const;
486+
std::pair<double,double>
487+
get_gravity_heating(const double Tc, const double r, const double X) const;
485488

486489
/**
487-
* Calculate energy (@p Qk) and entropy (@p Ek) change rate factor
490+
* Calculate entropy (@p Ek) and energy (@p Qk) change rate factor
488491
* (regarding the core cooling rate Tc/dt) for a given @p Tc (CMB temperature)
489492
*/
490-
void get_adiabatic_heating(const double Tc, double &Ek, double &Qk) const;
493+
std::pair<double,double>
494+
get_adiabatic_heating(const double Tc) const;
491495

492496
/**
493-
* Calculate energy (@p Ql) and entropy (@p El) change rate factor
497+
* Calculate entropy (@p El) and energy (@p Ql) change rate factor
494498
* (regarding the inner core growth rate dR/dt) for a given @p Tc (CMB temperature)
495499
* and @p r (inner core radius)
496500
*/
497-
void get_latent_heating(const double Tc, const double r, double &El, double &Ql) const;
501+
std::pair<double,double>
502+
get_latent_heating(const double Tc, const double r) const;
498503

499504
/**
500505
* Calculate entropy of heat of solution @p Eh for a given @p Tc (CMB temperature),
501506
* @p r (inner core radius), and @p X (light element concentration)
502507
*/
503-
void get_heat_solution(const double Tc, const double r, const double X, double &Eh) const;
508+
double
509+
get_heat_solution(const double Tc, const double r, const double X) const;
504510

505511
/**
506512
* return radiogenic heating rate at the current time

‎source/boundary_temperature/dynamic_core.cc‎

Lines changed: 94 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -512,12 +512,12 @@ namespace aspect
512512
void
513513
DynamicCore<dim>::update_core_data()
514514
{
515-
get_specific_heating(core_data.Ti,core_data.Qs,core_data.Es);
516-
get_radio_heating(core_data.Ti,core_data.Qr,core_data.Er);
517-
get_gravity_heating(core_data.Ti,core_data.Ri,core_data.Xi,core_data.Qg,core_data.Eg);
518-
get_adiabatic_heating(core_data.Ti,core_data.Ek,core_data.Qk);
519-
get_latent_heating(core_data.Ti,core_data.Ri,core_data.El,core_data.Ql);
520-
get_heat_solution(core_data.Ti,core_data.Ri,core_data.Xi,core_data.Eh);
515+
std::tie(core_data.Qs,core_data.Es) = get_specific_heating(core_data.Ti);
516+
std::tie(core_data.Qr,core_data.Er) = get_radio_heating(core_data.Ti);
517+
std::tie(core_data.Qg,core_data.Eg) = get_gravity_heating(core_data.Ti,core_data.Ri,core_data.Xi);
518+
std::tie(core_data.Ek,core_data.Qk) = get_adiabatic_heating(core_data.Ti);
519+
std::tie(core_data.El,core_data.Ql) = get_latent_heating(core_data.Ti,core_data.Ri);
520+
core_data.Eh = get_heat_solution(core_data.Ti,core_data.Ri,core_data.Xi);
521521
}
522522

523523

@@ -660,51 +660,83 @@ namespace aspect
660660

661661

662662
template <int dim>
663-
void
663+
std::pair<double,double>
664664
DynamicCore<dim>::
665-
get_specific_heating(const double Tc, double &Qs,double &Es) const
665+
get_specific_heating(const double Tc) const
666666
{
667+
// The object is initialized with invalid values. Make sure that
668+
// by the time we get here, everything we need has been set to
669+
// reasonable values.
670+
Assert (numbers::is_finite(Tc), ExcInternalError());
671+
Assert (numbers::is_finite(L), ExcInternalError());
672+
Assert (numbers::is_finite(D), ExcInternalError());
673+
Assert (numbers::is_finite(Rho_cen), ExcInternalError());
674+
Assert (numbers::is_finite(Rc), ExcInternalError());
675+
667676
const double A = std::sqrt(1./(Utilities::fixed_power<-2>(L)+Utilities::fixed_power<-2>(D)));
668677
const double Is = 4.*numbers::PI*get_T(Tc,0.)*Rho_cen*(-Utilities::fixed_power<2>(A)*Rc/2.*std::exp(-Utilities::fixed_power<2>(Rc/A))+Utilities::fixed_power<3>(A)*std::sqrt(numbers::PI)/4.*std::erf(Rc/A));
669678

670-
Qs = -Cp/Tc*Is;
671-
Es = Cp/Tc*(Mc-Is/Tc);
679+
const double Qs = -Cp/Tc*Is;
680+
const double Es = Cp/Tc*(Mc-Is/Tc);
681+
682+
return { Qs, Es };
672683
}
673684

674685

675686

676687
template <int dim>
677-
void
688+
std::pair<double,double>
678689
DynamicCore<dim>::
679-
get_radio_heating(const double Tc, double &Qr, double &Er) const
690+
get_radio_heating(const double Tc) const
680691
{
692+
// The object is initialized with invalid values. Make sure that
693+
// by the time we get here, everything we need has been set to
694+
// reasonable values.
695+
Assert (numbers::is_finite(Tc), ExcInternalError());
696+
Assert (numbers::is_finite(L), ExcInternalError());
697+
Assert (numbers::is_finite(D), ExcInternalError());
698+
Assert (numbers::is_finite(Rho_cen), ExcInternalError());
699+
Assert (numbers::is_finite(Rc), ExcInternalError());
700+
681701
double It = numbers::signaling_nan<double>();
682702
if (D>L)
683703
{
684-
const double B = std::sqrt(1/(1/Utilities::fixed_power<2>(L)-1/Utilities::fixed_power<2>(D)));
685-
It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(-Utilities::fixed_power<2>(B)*Rc/2*std::exp(-Utilities::fixed_power<2>(Rc/B))+Utilities::fixed_power<3>(B)/std::sqrt(numbers::PI)/4*std::erf(Rc/B));
704+
const double B = std::sqrt(1./(1./Utilities::fixed_power<2>(L)-1./Utilities::fixed_power<2>(D)));
705+
It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(-Utilities::fixed_power<2>(B)*Rc/2*std::exp(-Utilities::fixed_power<2>(Rc/B))
706+
+ Utilities::fixed_power<3>(B)/std::sqrt(numbers::PI)/4*std::erf(Rc/B));
686707
}
687708
else
688709
{
689710
const double B = std::sqrt(1/(Utilities::fixed_power<-2>(D)-Utilities::fixed_power<-2>(L)));
690-
It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(Utilities::fixed_power<2>(B)*Rc/2*std::exp(Utilities::fixed_power<2>(Rc/B))-Utilities::fixed_power<2>(B)*fun_Sn(B,Rc,100)/2);
711+
It = 4*numbers::PI*Rho_cen/get_T(Tc,0)*(Utilities::fixed_power<2>(B)*Rc/2*std::exp(Utilities::fixed_power<2>(Rc/B))
712+
- Utilities::fixed_power<2>(B)*fun_Sn(B,Rc,100)/2);
691713
}
692714

693-
Qr = Mc*core_data.H;
694-
Er = (Mc/Tc-It)*core_data.H;
715+
const double Qr = Mc*core_data.H;
716+
const double Er = (Mc/Tc-It)*core_data.H;
717+
return { Qr, Er };
695718
}
696719

697720

698721

699722
template <int dim>
700-
void
723+
double
701724
DynamicCore<dim>::
702-
get_heat_solution(const double Tc, const double r, const double X, double &Eh) const
725+
get_heat_solution(const double Tc, const double r, const double X) const
703726
{
727+
// The object is initialized with invalid values. Make sure that
728+
// by the time we get here, everything we need has been set to
729+
// reasonable values.
730+
Assert (numbers::is_finite(Tc), ExcInternalError());
731+
Assert (numbers::is_finite(r), ExcInternalError());
732+
Assert (numbers::is_finite(X), ExcInternalError());
733+
Assert (numbers::is_finite(L), ExcInternalError());
734+
Assert (numbers::is_finite(Rc), ExcInternalError());
735+
704736
if (r==Rc)
705737
{
706738
// No energy change rate if the inner core is fully frozen.
707-
Eh = 0.;
739+
return 0.;
708740
}
709741
else
710742
{
@@ -722,17 +754,29 @@ namespace aspect
722754
It -= 4*numbers::PI*Rho_cen/get_T(Tc,0)*(Utilities::fixed_power<2>(B)*r/2*std::exp(Utilities::fixed_power<2>(r/B))-Utilities::fixed_power<2>(B)*fun_Sn(B,r,100)/2);
723755
}
724756
const double Cc = 4*numbers::PI*Utilities::fixed_power<2>(r)*get_rho(r)*X/(Mc-get_mass(r));
725-
Eh = Rh*(It-(Mc-get_mass(r))/get_T(Tc,r))*Cc;
757+
return Rh*(It-(Mc-get_mass(r))/get_T(Tc,r))*Cc;
726758
}
727759
}
728760

729761

730762

731763
template <int dim>
732-
void
764+
std::pair<double,double>
733765
DynamicCore<dim>::
734-
get_gravity_heating(const double Tc, const double r, const double X, double &Qg, double &Eg) const
766+
get_gravity_heating(const double Tc, const double r, const double X) const
735767
{
768+
// The object is initialized with invalid values. Make sure that
769+
// by the time we get here, everything we need has been set to
770+
// reasonable values.
771+
Assert (numbers::is_finite(Tc), ExcInternalError());
772+
Assert (numbers::is_finite(r), ExcInternalError());
773+
Assert (numbers::is_finite(X), ExcInternalError());
774+
Assert (numbers::is_finite(Rc), ExcInternalError());
775+
Assert (numbers::is_finite(r), ExcInternalError());
776+
Assert (numbers::is_finite(L), ExcInternalError());
777+
Assert (numbers::is_finite(Beta_c), ExcInternalError());
778+
779+
double Qg;
736780
if (r==Rc)
737781
Qg = 0.;
738782
else
@@ -747,29 +791,47 @@ namespace aspect
747791
-(Mc-get_mass(r))*get_gravity_potential(r))*Beta_c*Cc;
748792
}
749793

750-
Eg = Qg/Tc;
794+
const double Eg = Qg/Tc;
795+
796+
return { Qg, Eg };
751797
}
752798

753799

754800

755801
template <int dim>
756-
void
802+
std::pair<double,double>
757803
DynamicCore<dim>::
758-
get_adiabatic_heating(const double Tc, double &Ek, double &Qk) const
804+
get_adiabatic_heating(const double Tc) const
759805
{
760-
Ek = 16*numbers::PI*k_c*Utilities::fixed_power<5>(Rc)/5/Utilities::fixed_power<4>(D);
761-
Qk = 8*numbers::PI*Utilities::fixed_power<3>(Rc)*k_c*Tc/Utilities::fixed_power<2>(D);
806+
// The object is initialized with invalid values. Make sure that
807+
// by the time we get here, everything we need has been set to
808+
// reasonable values.
809+
Assert (numbers::is_finite(Tc), ExcInternalError());
810+
Assert (numbers::is_finite(Rc), ExcInternalError());
811+
Assert (numbers::is_finite(D), ExcInternalError());
812+
813+
const double Ek = 16*numbers::PI*k_c*Utilities::fixed_power<5>(Rc)/5/Utilities::fixed_power<4>(D);
814+
const double Qk = 8*numbers::PI*Utilities::fixed_power<3>(Rc)*k_c*Tc/Utilities::fixed_power<2>(D);
815+
return { Ek, Qk };
762816
}
763817

764818

765819

766820
template <int dim>
767-
void
821+
std::pair<double,double>
768822
DynamicCore<dim>::
769-
get_latent_heating(const double Tc, const double r, double &El, double &Ql) const
823+
get_latent_heating(const double Tc, const double r) const
770824
{
771-
Ql = 4.*numbers::PI*Utilities::fixed_power<2>(r)*Lh*get_rho(r);
772-
El = Ql*(get_T(Tc,r)-Tc)/(Tc*get_T(Tc,r));
825+
// The object is initialized with invalid values. Make sure that
826+
// by the time we get here, everything we need has been set to
827+
// reasonable values.
828+
Assert (numbers::is_finite(Tc), ExcInternalError());
829+
Assert (numbers::is_finite(r), ExcInternalError());
830+
Assert (numbers::is_finite(Lh), ExcInternalError());
831+
832+
const double Ql = 4.*numbers::PI*Utilities::fixed_power<2>(r)*Lh*get_rho(r);
833+
const double El = Ql*(get_T(Tc,r)-Tc)/(Tc*get_T(Tc,r));
834+
return { El, Ql };
773835
}
774836

775837

0 commit comments

Comments
 (0)