Does the option 'dOutputTime' influence the calculated data in the module magmoc?

[zjdzzz123]

When I tried to modify the parameters in Example MagmOc_Earth to use on other planets, I found that changing option dOutputTime in vpl.in would cause the output data to change. This change does not seem to be caused by computational precision. See Figures 1-4 below for details, which correspond to the cases where dOutputTime is 1e5, 1e4, 1e3, and 1e2 respectively. Please forgive me for using screenshots to display the pictures.

Figure 1 dOutputTime=1e5

Figure 2 dOutputTime=1e4

Figure 3 dOutputTime=1e3

Figure 4 dOutputTime=1e2

It seems that option dOutputTime is not only used to control the output step size, but also to control the calculation step size. In my opinion, data calculation should not be related to the output frequency of the data, so modifying option d will only affect the accuracy of the image, and will not produce such abnormal oscillations. So I raised this issue, hoping someone can point out whether it is a problem with my code or something else causing this.

Here is the code I used:
In file vpl.in:

sSystemName   System1  # change to System2, System3, ... to adapt python code
iVerbose      5                            
bOverwrite    1                             

saBodyFiles   star.in b.in

sUnitMass      solar                        
sUnitLength    aU                           
sUnitTime      YEARS                       
sUnitAngle     d                            

bDoLog         1                            
iDigits            6                            
dMinValue     1e-10                        

bDoForward    1                             
bVarDt             1                             
dEta                 1.0
dStopTime       1e8
dOutputTime   1e2           # I only changed this option from 1e5 to 1e2

In file star.in:

# The host star
sName	        star                # Body's name
saModules	stellar, eqtide       # Modules to apply, exact spelling required
# saModules	stellar       # Modules to apply, exact spelling required

# Output
saOutputOrder Age -Luminosity -LXUVStellar -Radius Temperature

# Physical Parameters
dMass         0.333                       # Mass, solar masses
dAge          5e6

# STELLAR Parameters
sStellarModel             baraffe                     # Stellar evolution model: `baraffe` or `none`
sMagBrakingModel          reiners
# dSatXUVFrac   1.e-3.03                       # Saturation level of the XUV luminosity
# dSatXUVTime   -3.146			  # End of XUV-Saturation (- in Gyr, default .1Gyr)
# dXUVBeta      1.172

# EQTIDE Parameters
dTidalQ         100          # Tidal phase lag
dK2             0.5          # Love number of degree 2
dMaxLockDiff    0.01
sTideModel      p2           # Tidal model, p2=CPL, t8=CTL
saTidePerts     b        # Body name(s) of tidal perturbers

In file b.in:

# Planet a parameters
sName         b                     # Body's name
saModules     magmoc atmesc eqtide                      # Modules to apply, exact spelling required

# Physical Properties
dMass         -2.82                       # Mass, negative -> Earth masses
sPlanetRadiusModel      lopez12      # Mass-radius model

# MAGMOC Properties
dMassFracFeOIni		0.0788		  # Initial mass fraction of FeO in the mantle (default: 0.0788)
dWaterMassAtm		-500		  # Initial water mass in the system (neg: terrestrial oceans, default: -1)
dSurfTemp		3000		  # Initial surface temperature (default: 4000K)
dManMeltDensity		4000		  # Mantle melt density (default: 4000kg/m^3)
bHaltMantleSolidified 	1		  # Halt when mantle solidified? (default: 0 = no, 1 = yes)
bHaltMantleMeltFracLow  1     		  # Halt when melt fraction smaller than 0.4 at surface? (default: 0 = no, 1 = yes)
bHaltAtmDesiSurfCool    1                 # Halt when atmosphere desiccated? (default: 0 = no, 1 = yes)
bHaltEnterHabZone	0		  # Halt when planet enters habitable zone? (default: 0 = no, 1 = yes)
sRadioHeatModel		schaefer          # Radiogenic heating model (default: none = RadHeat; schaefer = BSE composition)
sMagmOcAtmModel 	grey		  # Atmospheric net flux model (default: grey; petit only for GJ1132b-H2O)
bOptManQuasiSol		1

# ATMESC Properties
dXFrac                    1.0             # X-Ray/XUV absorption radius (fraction of planet radius)
dSurfWaterMass            -1.0            # Initial surface water (Earth oceans)
dEnvelopeMass             0               # Initial envelope mass (Earth masses)
bHaltSurfaceDesiccated    0               # Halt when dry?
bHaltEnvelopeGone         0               # Halt when evaporated?
dMinSurfWaterMass         -1.e-5          # Planet is desiccated when water content drops below this (Earth oceans)
sWaterLossModel           lbexact
# sPlanetRadiusModel        none
bInstantO2Sink            0
sAtmXAbsEffH2OModel       none
dAtmXAbsEffH2O 		  0.3

# EQTIDE Parameters
# dTidalQ         100		   # Tidal phase lag
dTidalQEnv      100
# dK2             0.5	           # Love number of degree 2
dK2Env          0.5
dMaxLockDiff    0.01
saTidePerts     star          # Body name(s) of tidal perturbers
bFixOrbit 	1		   # Keep dEcc constant?
bForceEqSpin	1		   # Planet tidally locked, no rotation?

# Orbital Properties
dSemi         -0.01713              # Semi-major axis, negative -> AU
dEcc          0.0                # Eccentricity
# dOrbPeriod    -1.510826

# Output
saOutputOrder Time -PotTemp -SurfTemp -SolidRadius -WaterMassMOAtm -WaterMassSol -OxygenMassMOAtm -OxygenMassSol -PressWaterAtm -PressOxygenAtm -HydrogenMassSpace -OxygenMassSpace FracFe2O3Man NetFluxAtmo MeltFraction $
WaterFracMelt

I have write a python script for my own debug, the python code I used to draw the picture:

OutPutString = 'saOutputOrder Time -PotTemp -SurfTemp -SolidRadius -WaterMassMOAtm -WaterMassSol -OxygenMassMOAtm -OxygenMassSol -PressWaterAtm -PressOxygenAtm -HydrogenMassSpace -OxygenMassSpace FracFe2O3Man NetFluxAtmo MeltFraction WaterFracMelt'
OutPutString_Star = 'saOutputOrder Age -Luminosity -LXUVStellar -Radius Temperature'

OutPutString = OutPutString.replace('saOutputOrder', '')
OutPutString_Star = OutPutString_Star.replace('saOutputOrder', '')

OutPutString = OutPutString.replace('$', '').split()
OutPutString_Star = OutPutString_Star.replace('$', '').split()

OutPutList = [item.lstrip('-') for item in OutPutString]
OutPutList_Star = [item.lstrip('-') for item in OutPutString_Star]
# print(OutPutList)

import matplotlib.pyplot as plt
import numpy as np

filename = input("the number of model is:")

file_path = f'System{ filename }.b.forward'
file_path_star = f'System{ filename }.star.forward'

data = np.loadtxt(file_path)
data_star = np.loadtxt(file_path_star)


time = data[:, 0]
time_star = data_star[:, 0]


columns = OutPutList
columns_star = OutPutList_Star

fig, axes = plt.subplots(nrows=5, ncols=4, figsize=(18, 15))
axes = axes.flatten()

# Caculation
n_time = len(time)
TO = 1.39e21  # mass of 1 Terr. Ocean [kg]
REARTH = 6.3781e6  # Earth radius [m]
MEARTH = 5.972186e24  # Earth mass [kg]
BIGG = 6.67428e-11  # Gravitational constant [m**3/kg/s**2]
r_p = 1.305 * REARTH  # Planet radius [m]
r_core = 3.4e6 / r_p
m_p = 2.82 * MEARTH  # Planet mass [kg]
g = (BIGG * m_p) / (r_p ** 2)  # Surface gravity [m/s**2]

for i in range(n_time):
    data[i, 8] = data[i, 8] * 1e5 * 4 * np.pi * r_p ** 2 / g
    data[i, 9] = data[i, 9] * 1e5 * 4 * np.pi * r_p ** 2 / g

columns[8] = 'WaterAtm'
columns[9] = 'OxygenAtm'

for i in range(1, data.shape[1]):
    axes[i-1].plot(time, data[:, i])
    axes[i-1].set_title(columns[i])
    axes[i-1].set_xlabel('Time')
    axes[i-1].set_ylabel('Value')
    axes[i-1].grid(True)

for i in range(1, data_star.shape[1]):
    axes[i-1+data.shape[1]].plot(time_star, data_star[:, i])
    axes[i-1+data.shape[1]].set_title(columns_star[i])
    axes[i-1+data.shape[1]].set_xlabel('Time')
    axes[i-1+data.shape[1]].set_ylabel('Value')
    axes[i-1+data.shape[1]].grid(True)

plt.tight_layout()
plt.show()