Dear users
I wanted to use the “dynamic core” option to cool my CMB. Like some other authors did in the literature, i wanted to start using a wide range of initial temperatures (6000 K, 5000K, 4000 K) and see how the temperature evolves with time in my thermochemical model. However, if I exceed around 4200 K, ASPECT says that it didn’t find a solution for the inner core.
The point is that other simulations done with StagYY (literature) had this purpose: see how the temperature of the CMB evolved independently if the inner core is present or not (it will eventually form in the future).
Is there someway to use this option without having the necessity of having an inner core? I didn’t change any option from the standard. Also, the “test” file provided by ASPECT start from 4100K, while in the manual there is written that the initial temperature is 6000 K.
subsection Dynamic core
set Inner temperature = 4100
set Outer temperature = 300
set Core density = 12.5e3
set Gravity acceleration = 9.81 # For the core
set CMB pressure = 0.14e12
set Initial light composition = 0.042 # wt% for light elements
set Max iteration = 30
set Core heat capacity = 840
set K0 = 4.111e11
set Alpha = 1.35e-5 # Thermal expansivity for the core
set Core conductivity = 50
subsection Geotherm parameters
set Tm0= 1695
set Tm1= 10.9e-12
set Tm2= -8.0e-24
set Theta= 0.11
set Composition dependency= false
end
subsection Radioactive heat source
# 100ppm K in core
set Number of radioactive heating elements = 1
set Heating rates = 2.92e-5
set Half life times = 1.25 #Unit: Gyr
set Initial concentrations = 0.1442962
end
Thank you in advance
Francesco
@Francyrad Like you, I see no particularly good reason why the core initial temperature should matter. Could you post the exact input file you are using, plus the error you are getting?
As for the difference between the manual and the corresponding test file, would you mind opening a github issue that points at the two files that should be, but are not, in agreement?
Best
W.
@bangerth thank you for the reply. This is the .prm file that i use:
# A simple setup for Earth convection model in a 2d shell
# with the core mantle boundary (CMB) temperature dynamically evolves through time.
# The 'dynamic core' boundary temperature plugin and 'core statistics' postprocessor
# are used to solving and tracing the core evolution.
set Dimension = 2
set Use years in output instead of seconds = true
set End time = 4.5e9
set Output directory = dynamic_core_1
#set Resume computation = true
subsection Material model
set Model name = simple compressible
subsection Simple compressible model
set Thermal expansion coefficient = 4e-5
set Viscosity = 1e22
end
end
subsection Geometry model
set Model name = spherical shell
subsection Spherical shell
set Inner radius = 3481000
set Outer radius = 6336000
end
end
subsection Boundary velocity model
set Zero velocity boundary indicators = inner
set Tangential velocity boundary indicators = outer
end
subsection Boundary temperature model
set Fixed temperature boundary indicators = inner, outer
set List of model names = dynamic core
subsection Dynamic core
set Inner temperature = 4100 # or 6000, like the default value. Above this value the simulation will fail.
set Outer temperature = 300
set Core density = 12.5e3
set Gravity acceleration = 9.81 # For the core
set CMB pressure = 0.14e12
set Initial light composition = 0.042 # wt% for light elements
set Max iteration = 30
set Core heat capacity = 840
set K0 = 4.111e11
set Alpha = 1.35e-5 # Thermal expansivity for the core
set Core conductivity = 50
subsection Geotherm parameters
set Tm0= 1695
set Tm1= 10.9e-12
set Tm2= -8.0e-24
set Theta= 0.11
set Composition dependency= false
end
subsection Radioactive heat source
# 100ppm K in core
set Number of radioactive heating elements = 1
set Heating rates = 2.92e-5
set Half life times = 1.25 #Unit: Gyr
set Initial concentrations = 0.1442962
end
end
end
subsection Heating model
set List of model names = adiabatic heating
end
subsection Initial temperature model
set Model name = spherical hexagonal perturbation
end
subsection Gravity model
set Model name = ascii data
end
subsection Mesh refinement
set Initial global refinement = 5
set Initial adaptive refinement = 0 #4
set Time steps between mesh refinement = 0 #15
end
subsection Postprocess
set List of postprocessors = core statistics, visualization, velocity statistics, temperature statistics, heat flux statistics, depth average
subsection Dynamic core statistics
set Excess entropy only = true
end
subsection Visualization
set Output format = vtu
set Time between graphical output = 1e6
set Number of grouped files = 1
end
subsection Depth average
set Time between graphical output = 1.5e6
set Output format = vtu
end
end
subsection Checkpointing
set Steps between checkpoint = 50
end
subsection Termination criteria
set Termination criteria = user request
end
The error that i get is:
An error occurred in line <494> of file </home/francesco-radica/Documenti/aspect/source/boundary_temperature/dynamic_core.cc> in function
bool aspect::BoundaryTemperature::DynamicCore<dim>::solve_time_step(double&, double&, double&) [with int dim = 2]
The violated condition was:
false
Additional information:
[Dynamic core] No inner core radius solution found!
Stacktrace:
-----------
#0 aspect: aspect::BoundaryTemperature::DynamicCore<2>::solve_time_step(double&, double&, double&)
#1 aspect: aspect::BoundaryTemperature::DynamicCore<2>::update()
#2 aspect: aspect::BoundaryTemperature::Manager<2>::update()
#3 aspect: aspect::Simulator<2>::compute_current_constraints()
#4 aspect: aspect::Simulator<2>::start_timestep()
#5 aspect: aspect::Simulator<2>::run()
#6 aspect: void run_simulator<2>(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, bool, bool, bool)
#7 aspect: main
--------------------------------------------------------
Aborting!
----------------------------------------------------
@Francyrad I don’t know that plugin, and I’m not sure anyone else does well enough to really help you. But you might want to look into the code to see what can be done. The function that triggers the error solves a nonlinear problem by way of the bisection method. What the error reports is that for the values you use, it cannot seem to find a solution to the nonlinear problem. I suspect that the output you should have gotten just above the error you show might illuminate why that is so.
Best
W.
Are the roots of the problem just these lines?:
template <int dim>
bool
DynamicCore<dim>::solve_time_step(double &X, double &T, double &R)
{
// Well solving the change in core-mantle boundary temperature T, inner core radius R, and
// light component (e.g. S, O, Si) composition X, the following relations has to be respected:
// 1. At the inner core boundary the adiabatic temperature should be equal to solidus temperature
// 2. The following energy production rate should be balanced in core:
// Heat flux at core-mantle boundary Q
// Specific heat Qs*dT/dt
// Radioactive heating Qr
// Gravitational contribution Qg*dR/dt
// Latent heat Ql*dR/dt
// So that Q+Qs*dT/dt+Qr+Qg*dR/dt*Ql*dR/dt=0
// 3. The light component composition X depends on inner core radius (See function get_X() ),
// and core solidus may dependent on X as well
// This becomes a small nonlinear problem. Directly iterate through the above three system doesn't
// converge well. Alternatively we solve the inner core radius by bisection method.
int steps=1;
double R_0,R_1,R_2;
// dT is the temperature difference between adiabatic and solidus at
// inner-outer core boundary. If dT=0 then we found our solution.
double dT0,dT1,dT2;
R_0 = 0.;
R_1 = core_data.Ri;
R_2 = Rc;
dT0 = get_dT(R_0);
dT1 = get_dT(R_1);
dT2 = get_dT(R_2);
if (dT0 >= 0. && dT2 >= 0.)
{
// Fully molten core
R_1 = R_0;
dT1 = 0;
}
else if (dT2 <= 0. && dT0 <= 0. )
{
//Completely solid core
R_1 = R_2;
dT1 = 0;
}
else while (!(dT1==0 || steps>max_steps))
{
// If solution is out of the interval, then something is wrong.
if (dT0*dT2>0)
{
this->get_pcout()<<"Step: "<<steps<<std::endl
<<" R=["<<R_0/1e3<<","<<R_2/1e3<<"]"<<"(km)"
<<" dT0="<<dT0<<", dT2="<<dT2<<std::endl
<<"Q_CMB="<<core_data.Q<<std::endl
<<"Warning: Solution for inner core radius can not be found! Mid-point is used."<<std::endl;
AssertThrow(dT0*dT2<=0,ExcMessage("No single solution for inner core!"));
}
else if (dT0*dT1<0.)
{
R_2 = R_1;
dT2 = dT1;
}
else if (dT2*dT1<0.)
{
R_0 = R_1;
dT0 = dT1;
}
R_1 = (R_0 + R_2) / 2.;
dT1 = get_dT(R_1);
++steps;
}
// Calculate new R,T,X
R = R_1;
T = get_Tc(R);
X = get_X(R);
if (dT0<0. && dT2>0.)
{
// Normal solution
return true;
}
else if (dT0>0. && dT2<0.)
{
// Snowing core solution
return false;
}
else
{
// No solution found.
this->get_pcout()<<"[Dynamic core] Step: "<<steps<<std::endl
<<" R=["<<R_0/1e3<<","<<R_2/1e3<<"]"<<"(km)"
<<" dT0="<<dT0<<", dT2="<<dT2<<std::endl
<<"Q_CMB="<<core_data.Q<<std::endl;
AssertThrow(false, ExcMessage("[Dynamic core] No inner core radius solution found!"));
}
return false;
}
That’s the file/function where the error message is produced. I don’t know what algorithm it uses, but you could presumably go through the function (it’s not that long) and try to understand why it does not want to solve the problem you have. Understanding is the first step towards fixing.
Best
W.
Hello there! I ran into the same problem as well when I was modeling the lunar core evolution. Although I didn’t figure out the reason, I did find an another problem.In the test the Tm1= 10.9e-12 and Tm2= -8.0e-24 shouldn’t be Tm1= 10.9 and Tm2= -8.0? Since their units are TPa−1 and TPa−2, and also the default values are 10.8 and -8.0. However, even after changing them to 10.8 and -8.0 it reported the same error again .
For what I understood the simulation won’t run if you don’t have a solid core, the point is that i don’t care if there is or not. This thing is solvable just editing the code allowing to run the sim if the inner core is present or not, and if it will born to truck the growth like some other studies did.
I solved the problem, now the code works perfectly. I please ask to accept the new PR
The previous code would halt if no solution was found for the inner core, preventing the simulation from proceeding. This limitation made it impossible to observe the Earth’s mantle evolution during cooling when the inner core had not yet formed.
The updated code addresses this issue. Now, if the inner core has not formed, the simulation continues, allowing the core-mantle boundary (CMB) to cool. If the inner core eventually forms during the simulation, its radius will be recorded in the statistics.
![image|690x382](upload://pA8DAinedKNwYdacDMY722QYu8s.png
I updated the code using english, in PR is not in italian
best
Francyrad
Hi @Francyrad, I tried to reproduce your error message while reviewing your PR, but for me the model runs fine for at least 50 timesteps. Can you post the complete log.txt file that ASPECT writes during the run that fails? I will still finish the review for now, but I have some questions that I can only check if I can reproduce the problem. Thanks