Merge pull request #297 from VirtualPlanetaryLaboratory/EqtideThermint

Fixed Eqtide-Thermint coupling so the CPL model can work with Thermint.
VirtualPlanetaryLaboratory · Jul 30, 2024 · cf99e6a · cf99e6a
2 parents 4c18ee4 + 0aaf16b
commit cf99e6a
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