Fix bug for fully solid case

source/boundary_temperature/dynamic_core.cc

@@ -701,21 +701,29 @@ namespace aspect
     DynamicCore<dim>::
     get_heat_solution(const double Tc, const double r, const double X, double &Eh) const
     {
-      double It = numbers::signaling_nan<double>();
-      if (D>L)
+      if (r==Rc)
         {
-          const double B = std::sqrt(1./(1./Utilities::fixed_power<2>(L)-1./Utilities::fixed_power<2>(D)));
-          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));
-          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<3>(B)/std::sqrt(numbers::PI)/4*std::erf(r/B));
+          // No energy change rate if the inner core is fully frozen.
+          Eh = 0.;
         }
       else
         {
-          const double B = std::sqrt(1./(Utilities::fixed_power<-2>(D)-Utilities::fixed_power<-2>(L)));
-          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);
-          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);
+          double It = numbers::signaling_nan<double>();
+          if (D>L)
+            {
+              const double B = std::sqrt(1./(1./Utilities::fixed_power<2>(L)-1./Utilities::fixed_power<2>(D)));
+              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));
+              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<3>(B)/std::sqrt(numbers::PI)/4*std::erf(r/B));
+            }
+          else
+            {
+              const double B = std::sqrt(1./(Utilities::fixed_power<-2>(D)-Utilities::fixed_power<-2>(L)));
+              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);
+              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);
+            }
+          const double Cc = 4*numbers::PI*Utilities::fixed_power<2>(r)*get_rho(r)*X/(Mc-get_mass(r));
+          Eh = Rh*(It-(Mc-get_mass(r))/get_T(Tc,r))*Cc;
         }
-      const double Cc = 4*numbers::PI*Utilities::fixed_power<2>(r)*get_rho(r)*X/(Mc-get_mass(r));
-      Eh = Rh*(It-(Mc-get_mass(r))/get_T(Tc,r))*Cc;
     }
 
 
@@ -725,12 +733,12 @@ namespace aspect
     DynamicCore<dim>::
     get_gravity_heating(const double Tc, const double r, const double X, double &Qg, double &Eg) const
     {
-      const double Cc = 4*numbers::PI*Utilities::fixed_power<2>(r)*get_rho(r)*X/(Mc-get_mass(r));
-      const double C_2 = 3./16.*Utilities::fixed_power<2>(L) - 0.5*Utilities::fixed_power<2>(Rc)*(1.-3./10.*Utilities::fixed_power<2>(Rc/L));
       if (r==Rc)
         Qg = 0.;
       else
         {
+          const double Cc = 4*numbers::PI*Utilities::fixed_power<2>(r)*get_rho(r)*X/(Mc-get_mass(r));
+          const double C_2 = 3./16.*Utilities::fixed_power<2>(L) - 0.5*Utilities::fixed_power<2>(Rc)*(1.-3./10.*Utilities::fixed_power<2>(Rc/L));
           Qg = (8./3.*Utilities::fixed_power<2>(numbers::PI*Rho_cen)*constants::big_g*(
                   ((3./20.*Utilities::fixed_power<5>(Rc)-Utilities::fixed_power<2>(L)*Utilities::fixed_power<3>(Rc)/8.-C_2*Utilities::fixed_power<2>(L)*Rc)*std::exp(-Utilities::fixed_power<2>(Rc/L))
                    +C_2/2.*Utilities::fixed_power<3>(L)*std::sqrt(numbers::PI)*std::erf(Rc/L))

tests/dynamic_core_fully_molten.prm

@@ -0,0 +1,12 @@
+# A simple setup for Earth convection model in a 2d shell
+# where the core mantle boundary (CMB) temperature dynamically evolves through time.
+# 
+# This is a variation of dynamic_core.prm with a fully molten core.
+
+include $ASPECT_SOURCE_DIR/tests/dynamic_core.prm
+
+subsection Boundary temperature model
+  subsection Dynamic core
+    set Inner temperature                 = 6000
+  end
+end

tests/dynamic_core_fully_molten/screen-output

@@ -0,0 +1,38 @@
+
+Number of active cells: 768 (on 4 levels)
+Number of degrees of freedom: 10,656 (6,528+864+3,264)
+
+   Dynamic core initialized as:
+     Tc(K)          Ri(km)         Xi             dT/dt(K/year)  dR/dt(km/year) dX/dt(1/year)  
+     6000           0              0.042          0              0              0              
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving Stokes system (GMG)... 16+0 iterations.
+
+   Postprocessing:
+     CMB heat flux out of the core 1.43 TW,
+
+*** Timestep 1:  t=598995 years, dt=598995 years
+   Dynamic core data updated.
+     Tc(K)          Ri(km)         Xi             dT/dt(K/year)  dR/dt(km/year) dX/dt(1/year)  
+     6000.07        0              0.042          1.12778e-07    0              0              
+   Solving temperature system... 14 iterations.
+   Solving Stokes system (GMG)... 16+0 iterations.
+
+   Postprocessing:
+     CMB heat flux out of the core 1.43 TW,
+
+*** Timestep 2:  t=1e+06 years, dt=401005 years
+   Dynamic core data updated.
+     Tc(K)          Ri(km)         Xi             dT/dt(K/year)  dR/dt(km/year) dX/dt(1/year)  
+     6000.11        0              0.042          1.14788e-07    0              0              
+   Solving temperature system... 17 iterations.
+   Solving Stokes system (GMG)... 15+0 iterations.
+
+   Postprocessing:
+     CMB heat flux out of the core 1.31 TW,
+
+Termination requested by criterion: end time
+
+
+

tests/dynamic_core_fully_molten/statistics

@@ -0,0 +1,18 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Iterations for temperature solver
+# 8: Iterations for Stokes solver
+# 9: Velocity iterations in Stokes preconditioner
+# 10: Schur complement iterations in Stokes preconditioner
+# 11: CMB heat flux out of the core (TW)
+# 12: CMB Temperature (K)
+# 13: Inner core radius (km)
+# 14: Light element concentration (%)
+# 15: Excess entropy (W/K)
+0 0.000000000000e+00 0.000000000000e+00 768 7392 3264  0 15 17 17 1.429e+00 6000.00 0.00 4.2000 -2.864e+07 
+1 5.989948582733e+05 5.989948582733e+05 768 7392 3264 14 15 17 17 1.429e+00 6000.07 0.00 4.2000 -1.799e+08 
+2 1.000000000000e+06 4.010051417267e+05 768 7392 3264 17 14 16 16 1.305e+00 6000.11 0.00 4.2000 -1.827e+08 

tests/dynamic_core_fully_solid.prm

@@ -0,0 +1,12 @@
+# A simple setup for Earth convection model in a 2d shell
+# where the core mantle boundary (CMB) temperature dynamically evolves through time.
+# 
+# This is a variation of dynamic_core.prm with a fully frozen core.
+
+include $ASPECT_SOURCE_DIR/tests/dynamic_core.prm
+
+subsection Boundary temperature model
+  subsection Dynamic core
+    set Inner temperature                 = 1000
+  end
+end

tests/dynamic_core_fully_solid/screen-output

@@ -0,0 +1,28 @@
+
+Number of active cells: 768 (on 4 levels)
+Number of degrees of freedom: 10,656 (6,528+864+3,264)
+
+   Dynamic core initialized as:
+     Tc(K)          Ri(km)         Xi             dT/dt(K/year)  dR/dt(km/year) dX/dt(1/year)  
+     1000           3481           0.084          0              0              0              
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving Stokes system (GMG)... 16+0 iterations.
+
+   Postprocessing:
+     CMB heat flux out of the core 0.175 TW,
+
+*** Timestep 1:  t=1e+06 years, dt=1e+06 years
+   Dynamic core data updated.
+     Tc(K)          Ri(km)         Xi             dT/dt(K/year)  dR/dt(km/year) dX/dt(1/year)  
+     1000.13        3481           0.084          1.33889e-07    0              0              
+   Solving temperature system... 9 iterations.
+   Solving Stokes system (GMG)... 16+0 iterations.
+
+   Postprocessing:
+     CMB heat flux out of the core 0.175 TW,
+
+Termination requested by criterion: end time
+
+
+

tests/dynamic_core_fully_solid/statistics

@@ -0,0 +1,17 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Iterations for temperature solver
+# 8: Iterations for Stokes solver
+# 9: Velocity iterations in Stokes preconditioner
+# 10: Schur complement iterations in Stokes preconditioner
+# 11: CMB heat flux out of the core (TW)
+# 12: CMB Temperature (K)
+# 13: Inner core radius (km)
+# 14: Light element concentration (%)
+# 15: Excess entropy (W/K)
+0 0.000000000000e+00 0.000000000000e+00 768 7392 3264 0 15 17 17 1.755e-01 1000.00 3481.00 8.4000  8.412e+08 
+1 1.000000000000e+06 1.000000000000e+06 768 7392 3264 9 15 17 17 1.755e-01 1000.13 3481.00 8.4000 -2.367e+08