Merge pull request #6224 from gassmoeller/fix_dynamic_core

Fix dynamic core boundary temperature plugin
geodynamics · Feb 7, 2025 · 4701e9c · 4701e9c
2 parents 2ac3b92 + 1665d04
commit 4701e9c
Show file tree

Hide file tree

Showing 8 changed files with 198 additions and 40 deletions.
diff --git a/doc/modules/changes/20250206_francyrad b/doc/modules/changes/20250206_francyrad
@@ -0,0 +1,5 @@
+Fixed: The boundary temperature plugin 'dynamic core' used to
+crash when the inner core was completely molten or completely
+solid. This is fixed now.
+<br>
+(Francesco Radica, Rene Gassmoeller, 2025/02/06)
diff --git a/source/boundary_temperature/dynamic_core.cc b/source/boundary_temperature/dynamic_core.cc
@@ -356,6 +356,8 @@ namespace aspect
       return (r0+r1)/2.;
     }
 
+
+
     template <int dim>
     bool
     DynamicCore<dim>::solve_time_step(double &X, double &T, double &R) const
@@ -386,6 +388,8 @@ namespace aspect
       double dT1 = get_dT(R_1);
       double dT2 = get_dT(R_2);
 
+      // If the temperature difference at the core-mantle boundary and at the
+      // inner-outer core boundary have the same sign, we have a fully molten or fully solid core.
       if (dT0 >= 0. && dT2 >= 0.)
         {
           // Fully molten core
@@ -399,48 +403,64 @@ namespace aspect
           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;
-          }
+        {
+          // Use bisection method to find R_1 such that dT1 = 0
+          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;
+                }
+
+              // Update R_1 and recalculate 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);
 
+      // Check the signs of dT at the boundaries to classify the solution
       if (dT0<0. && dT2>0.)
         {
-          // Normal solution
+          // Core partially molten, freezing from the inside, normal solution
           return true;
         }
       else if (dT0>0. && dT2<0.)
         {
-          // Snowing core solution
+          // Core partially molten, snowing core solution
           return false;
         }
+      else if (dT0 >= 0. && dT2 >= 0.)
+        {
+          // Core fully molten, normal solution
+          return true;
+        }
+      else if (dT0 <= 0. && dT2 <= 0.)
+        {
+          // Core fully solid, normal solution
+          return true;
+        }
       else
         {
           // No solution found.
@@ -681,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;
     }
 
 
@@ -705,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))

diff --git a/tests/dynamic_core_fully_molten.prm b/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
diff --git a/tests/dynamic_core_fully_molten/screen-output b/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
+
+
+
diff --git a/tests/dynamic_core_fully_molten/statistics b/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 
diff --git a/tests/dynamic_core_fully_solid.prm b/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
diff --git a/tests/dynamic_core_fully_solid/screen-output b/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
+
+
+
diff --git a/tests/dynamic_core_fully_solid/statistics b/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