Fixed Eqtide-Thermint coupling so the CPL model can work with Thermint.

Added clean exit if user tries to select CTL model with Thermint. Added unit tests for the CPL model with Radheat and Thermint.
VirtualPlanetaryLaboratory · Jul 23, 2024 · 0aaf16b · 0aaf16b
1 parent 4c18ee4
commit 0aaf16b
Show file tree

Hide file tree

Showing 7 changed files with 605 additions and 12 deletions.
diff --git a/examples/TidalEarth/README.rst b/examples/TidalEarth/README.rst
@@ -15,8 +15,8 @@ from its host star.
 ===================   ============
 
 This example shows the thermal, magnetic, and orbital evolution of an
-Earth-like planet using **ThermInt**, **RadHeat**, and **EqTide**. Here we use the "orbit-only" tidal model, which is highly idealized.
-The interior model is 1-D and many free parameters have been
+Earth-like planet using **ThermInt**, **RadHeat**, and **EqTide**. Here we use the "orbit-only" (DB15) tidal model, 
+which is highly idealized. The interior model is 1-D and many free parameters have been
 tuned to match `the modern Earth <../EarthInterior>`_.  This example reproduces the results of
 `Driscoll & Barnes (2015) <https://ui.adsabs.harvard.edu/abs/2015AsBio..15..739D/abstract>`_.  The generated
 figures should be compared to their Figures 4 and 5.

diff --git a/src/body.c b/src/body.c
@@ -771,15 +771,11 @@ double fdHflowSecMan(BODY *body, EVOLVE *evolve, int iBody) {
   double dHflowSecMan = 0;
 
   if (body[iBody].bThermint) {
-    dHflowSecMan += fdPowerThermint(body, iBody);
-    // dHflowSecMan += fdHflowSecManThermint(body,iBody);
+    dHflowSecMan += fdPowerThermint(body, iBody); 
   }
   if (body[iBody].bEqtide) {
-    // dHflowSecMan -= fdTidePower(body,iBody,evolve->iEqtideModel); // formerly
-    // dTidalPowerMan
-    dHflowSecMan -= body[iBody].dTidalPowMan;
-  }
-  // Should add RadHeat here
+    dHflowSecMan -= fdTidePower(body,iBody,evolve->iEqtideModel);
+  }  // XXX Should add RadHeat here
   return dHflowSecMan;
 }
 

diff --git a/src/module.c b/src/module.c
@@ -981,8 +981,6 @@ void VerifyModuleMultiRadheatThermint(BODY *body, UPDATE *update,
                                       OPTIONS *options, int iBody,
                                       int *iModuleProps, int *iModuleForce) {
 
-  /* This will need modification if material can move between layers */
-
   if (body[iBody].bThermint) {
     if (!body[iBody].bRadheat) {
       if (control->Io.iVerbose > VERBINPUT) {
@@ -1007,7 +1005,11 @@ void VerifyModuleMultiEqtideThermint(BODY *body, UPDATE *update,
                                      int iBody, int *iModuleProps,
                                      int *iModuleForce) {
 
-  // Initialize mantle tidal power
+  if (control->Evolve.iEqtideModel == CTL) {
+    fprintf(stderr,"ERROR: The CTL EqTide model cannot be coupled to ThermInt.\n");
+    exit(EXIT_INPUT);
+  }
+
   body[iBody].dTidalPowMan = 0;
 
   if (body[iBody].bEqtide && body[iBody].bThermint) {

diff --git a/tests/EqtideRadheatThermint/CPLRadheatThermint/c.in b/tests/EqtideRadheatThermint/CPLRadheatThermint/c.in
@@ -0,0 +1,66 @@
+# Planet c's parameters
+sName		c			# Body's name
+saModules 	radheat thermint eqtide
+
+# Physical Properties
+dMass       -2.8857142857142857
+dRadius		-1.367  			# Radius, negative -> Earth radii
+# Orbital Properties
+dEcc        0.08571428571428572
+dSemi       -0.03984 		# Semi-major axis, negative -> AU
+
+# mantle parameters that we will change
+#dActViscMan    500000
+
+
+# EQTIDE Parameters
+bForceEqSpin  1
+dObliquity    0
+dRadGyra     0.5
+dTidalQ      10     #
+dK2          0.299   # Love number of degree 2 (Yoder 1995)
+saTidePerts  sun    # Names of perturbing bodies (only central body may have >1)
+bFixOrbit    0
+#bTideLock   1
+
+# RADHEAT Parameters
+# *Num* are in numbers of atoms, negative -> Earth vals
+### 40K
+d40KPowerMan      -1
+d40KPowerCore     -1
+d40KPowerCrust    -1
+### 232Th
+d232ThPowerMan	  -1
+d232ThPowerCore	  -1
+d232ThPowerCrust  -1
+### 235U
+d235UPowerMan     -1
+d235UPowerCore	  -1
+d235UPowerCrust	  -1
+### 238U
+d238UPowerMan	    -1
+d238UPowerCore	  -1
+d238UPowerCrust	  -1
+
+### THERMINT inputs.
+dTMan               3500
+dTCore              5500
+dEruptEff           0.2
+dViscJumpMan        2.0
+dTrefLind           5600
+dImK2ManOrbModel    2
+dShModRef           1e5 
+dStiffness          1.71e13 
+dActViscMan         309285.71428571426
+
+saOutputOrder -Time -TMan -TUMan -TLMan -TCMB -TCore $
+    -HflowUMan -HflowMeltMan -RadPowerMan -RadPowerCore -RadPowerCrust $
+    -HflowCMB -HflowSecMan -HFlowSurf $
+    -TDotMan -TDotCore -TJumpLMan -TJumpUMan -RIC -RayleighMan -ViscUMan -ViscLMan $
+    -MeltMassFluxMan FMeltUMan $
+    -MagMom -CoreBuoyTherm -CoreBuoyCompo -CoreBuoyTotal -MagPauseRad $
+    -BLUMan -BLLMan $
+     K2 ImK2 SemiMajorAxis Eccentricity ShmodUMan MeanMotion ViscUManArr $
+     ChiOC ChiIC MassChiOC MassChiIC MassOC MassIC $
+    -RadPowerTotal SurfEnFluxEqtide -PowerEqtide TideLock TidalQ -SurfEnFluxTotal 
+
diff --git a/tests/EqtideRadheatThermint/CPLRadheatThermint/sun.in b/tests/EqtideRadheatThermint/CPLRadheatThermint/sun.in
@@ -0,0 +1,11 @@
+# host star parameters
+sName        sun
+dMass        0.118
+dRadius      0.00072354
+dLuminosity  5.50944e30
+sStellarModel none
+saModules    stellar eqtide
+saTidePerts  c
+
+dTidalQ       1e6
+dK2           1.5