Is it possible to use dynamic core option without having an earth solid core?

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.