Commit 9366c3ca authored by Laura Ketzer's avatar Laura Ketzer
Browse files

updated mass_evolution_function. test if problems are fixed.

parent 989ccadc
......@@ -78,7 +78,12 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, initial_
R_core = planet_object.core_radius # stop when this radius is reached!
dt = initial_step_size
i = 1 # counter
# NOTE: minimum and maximum step size are HARDCODED for now (see further down in code for details)
min_step_size = 1e-2
max_step_size = 10.
i = 1 # counter to track how many tracked RK iterations have been performed
while t <= t_final:
#print(i, ' - ', t, "- dt = ", dt)
################
......@@ -89,9 +94,6 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, initial_
# for first step use the parameters initialized above (beta=beta0, K=K0, M=M0)
###############################################################################################################
Mdot1 = mass_loss_rate_forward_LO14(t, epsilon, K_on, beta_on, planet_object, f_env, track_dict)
k1 = (dt*Myr_to_sec * Mdot1)/M_earth # mass lost in one timestep in earth masses
M_05k1 = M + 0.5 * k1
# when the initial time step is too large OR the planet mass becomes very close to the core mass, it can happen
# that one of the RK substeps leads to such a large mass lost that the new planet mass is smaller than the core mass.
......@@ -107,113 +109,152 @@ def mass_planet_RK4_forward_LO14(epsilon, K_on, beta_on, planet_object, initial_
# Laura: the second case should already be taken care? or maybe make step size go smaller and smaller again until
# this situation is reached with the minimum step size
print("1:", M_05k1)
M_env_05k1 = M_05k1 - M_core
f_env_05k1 = (M_env_05k1/M_05k1) * 100 # new envelope mass fraction
print("f_env1:", f_env_05k1)
Mdot2 = mass_loss_rate_forward_LO14(t+0.5*dt, epsilon, K_on, beta_on, planet_object, f_env_05k1, track_dict)
k2 = (dt*Myr_to_sec * Mdot2)/M_earth
M_05k2 = M + 0.5*k2
print("2:", M_05k2)
M_env_05k2 = M_05k2 - M_core
f_env_05k2 = (M_env_05k2/M_05k2)*100 # new mass fraction
print("f_env2:", f_env_05k2)
Mdot3 = mass_loss_rate_forward_LO14(t+0.5*dt, epsilon, K_on, beta_on, planet_object, f_env_05k2, track_dict)
k3 = (dt*Myr_to_sec * Mdot3)/M_earth
M_k3 = M + k3
print("3:", M_k3)
M_env_k3 = M_k3 - M_core
f_env_k3 = (M_env_k3/M_k3)*100 # new mass fraction
print("f_env3:", f_env_k3)
Mdot4 = mass_loss_rate_forward_LO14(t+dt, epsilon, K_on, beta_on, planet_object, f_env_k3, track_dict)
k4 = (dt*Myr_to_sec * Mdot4)/M_earth
###########################################################################################
M_lost = (k1 + 2*k2 + 2*k3 + k4)/6. # after time-step dt
print(M_lost)
###########################################################################################
# this part is new compared to the one used in the PAPER (there we used a fixed step size!)
###########################################################################################
# if radius change is too drastic, decrease step size
M_new = M + M_lost
M_env_new = M_new - M_core
f_env_new = (M_env_new/M_new)*100 # in %
print("first ", M_new, f_env_new)
R_new = plmo14.calculate_planet_radius(M_core, f_env_new, t, flux_at_planet_earth(planet_object.Lbol,
planet_object.distance),
planet_object.metallicity)
#print("time = ", t)
#print(R, R_new, f_env, f_env_new)
#print(abs((R-R_new)/R)*100)
# adjust step size, if radius change is too drastic or too little
if (abs((R-R_new)/R)*100 >= 1.) & (t < track_dict["t_curr"]) & (dt > 1e-2):#10**(-1)): # larger than 1%, make step size smaller
dt = dt/10.
#print('make step size smaller: ', dt)
# don't write anything to file, do iteration again with new step size
elif (abs((R-R_new)/R)*100 < (0.01)) & (t < track_dict["t_curr"]) & (dt < 10.): # smaller than 0.01%, make step size bigger, but not bigger than 10 Myr
dt = dt*10.
#print('make step size bigger: ', dt)
elif (abs((R-R_new)/R)*100 >= 1.) & (t >= track_dict["t_curr"]) & (dt > 1e-2): # larger than 1%, make step size smaller
dt = dt/10.
elif (abs((R-R_new)/R)*100 < (0.01)) & (t >= track_dict["t_curr"]) & (dt < 10.): # smaller than 0.01%, make step size bigger, but not bigger than 10 Myr
dt = dt*10
else: # if radius change is ok, check if there is still envelope left
if ((M + M_lost) - M_core) >= 0:
#print("envelope left")
# then planet still has some atmosphere left -> continue
M = M + M_lost # new planet mass (M_lost is negative)
M_arr.append(M)
M_env = M - M_core # new envelope mass
#print("(t_arr[-1]", t_arr[-1])
t = t_arr[-1] + dt #t_start + i*dt # t_i_plus_1 - update time value
# if t > track_dict["t_curr"]:
# print(dt, "new time: ", t)
t_arr.append(t) # new time t
Lx_arr.append(Lx_evo(t=t, track_dict=track_dict)) # Lx at new time t
while (envelope_left = True):
Mdot1 = mass_loss_rate_forward_LO14(t, epsilon, K_on, beta_on, planet_object, f_env, track_dict)
k1 = (dt*Myr_to_sec * Mdot1)/M_earth # mass lost in one timestep in earth masses
M_05k1 = M + 0.5 * k1
M_env_05k1 = M_05k1 - M_core
f_env_05k1 = (M_env_05k1/M_05k1) * 100 # new envelope mass fraction
print("1:", M_05k1)
print("f_env1:", f_env_05k1)
if (i == 1) and (f_env_05k1 < 0.):
# then I am still in the first RK iteration, and the initial step size was likely too large
# set step size to minumum
dt = min_step_size
break
Mdot2 = mass_loss_rate_forward_LO14(t+0.5*dt, epsilon, K_on, beta_on, planet_object, f_env_05k1, track_dict)
k2 = (dt*Myr_to_sec * Mdot2)/M_earth
M_05k2 = M + 0.5 * k2
M_env_05k2 = M_05k2 - M_core
f_env_05k2 = (M_env_05k2/M_05k2) * 100
print("2:", M_05k2)
print("f_env2:", f_env_05k2)
if (i == 1) and (f_env_05k2 < 0.):
# then I am still in the first RK iteration, and the initial step size was likely too large
# set step size to minumum
dt = min_step_size
break
Mdot3 = mass_loss_rate_forward_LO14(t+0.5*dt, epsilon, K_on, beta_on, planet_object, f_env_05k2, track_dict)
k3 = (dt*Myr_to_sec * Mdot3)/M_earth
M_k3 = M + k3
M_env_k3 = M_k3 - M_core
f_env_k3 = (M_env_k3/M_k3) * 100
print("3:", M_k3)
print("f_env3:", f_env_k3)
if (i == 1) and (f_env_k3 < 0.):
# then I am still in the first RK iteration, and the initial step size was likely too large
# set step size to minumum
dt = min_step_size
break
Mdot4 = mass_loss_rate_forward_LO14(t+dt, epsilon, K_on, beta_on, planet_object, f_env_k3, track_dict)
k4 = (dt*Myr_to_sec * Mdot4)/M_earth
# new envelope mass fraction:
f_env = (M_env/M)*100 # in %
#print("f_env = ", f_env)
# calculate new radius with new planet mass/envelope mass fraction & one time step later
R = plmo14.calculate_planet_radius(M_core, f_env, t, flux_at_planet_earth(planet_object.Lbol,
planet_object.distance),
planet_object.metallicity)
R_arr.append(R)
i = i+1 # update step to i+1
else:
# all atmosphere is gone -> terminate
M_lost = (k1 + 2*k2 + 2*k3 + k4)/6. # after time-step dt
print(M_lost)
###########################################################################################
# this part is new compared to the one used in the PAPER (there we used a fixed step size!)
###########################################################################################
# if radius change is too drastic, decrease step size
M_new = M + M_lost
M_env_new = M_new - M_core
if (M_env_new <= 0.) and (dt == min_step_size):
# if this condition is not satisfied, then we assume the RK iteration would
# remove all atmosphere and only the rocky core is left at t_i+dt;
# this terminates the code and returns the final planet properties
envelope_left = False # no atmosphere left; this breaks the while loop
print("Atmosphere has evaportated! Only bare rocky core left! STOP this madness!")
# print("params=", t, M, f_env, R)
# print("Lx", Lx_evo(t=t_arr[-1]+dt, track_dict=track_dict))
# if the stop criterium is reached, I add the core mass & core radius to the array at t_i+1
print("t = ", t_arr[-1]+dt)
# since the the stop criterium is reached, we assume at t_i+1 the planet only consists of the bare rocky
# core with the planet mass equal to the core mass and the planet radius equal to the core radius
t_arr.append(t_arr[-1]+dt)
M_arr.append(M_core)
print("M_core = ", M_core)
R_arr.append(R_core)
print("R_core = ", R_core)
Lx_arr.append(Lx_evo(t=t_arr[-1]+dt, track_dict=track_dict))
Lx_arr.append(Lx_evo(t=t_arr[-1]+dt, track_dict=track_dict))
#print("t = ", t_arr[-1]+dt)
#print("M_core = ", M_core)
#print("R_core = ", R_core)
return np.array(t_arr), np.array(M_arr), np.array(R_arr), np.array(Lx_arr)
print("Done!")
elif (M_env_new <= 0.) and (dt != min_step_size):
# planet close to evaporation
# set step size ti minimal value and rund RK iteration again
dt = min_step_size
break
# if you're still in the while loop at this point, then calculate new radius and check how
# drastic the radius change would be; adjust the step size if too drastic or too little
f_env_new = (M_env_new/M_new)*100 # in %
R_new = plmo14.calculate_planet_radius(M_core, f_env_new, t, flux_at_planet_earth(planet_object.Lbol,
planet_object.distance),
planet_object.metallicity)
# if radius change is larger than 1%, make step size smaller by factor 10
# or if radius change is smaller than 0.01%, make step size bigger by factor 10, but not bigger than 10 Myr!
# and don't write anything to file, instead do RK iteration again with new step size
if (abs((R-R_new)/R)*100 >= 1.) & (t < track_dict["t_curr"]) & (dt > min_step_size):
dt = dt / 10.
break
elif (abs((R-R_new)/R)*100 < (0.01)) & (t < track_dict["t_curr"]) & (dt < max_step_size):
dt = dt * 10.
break
# in principle I can adjust the code such that these hardcoded parameters are different for early planet
# evolution where much is happening, and late planet evolution where almost no change is occurring anymore
elif (abs((R-R_new)/R)*100 >= 1.) & (t >= track_dict["t_curr"]) & (dt > min_step_size):
dt = dt / 10.
break
elif (abs((R-R_new)/R)*100 < (0.01)) & (t >= track_dict["t_curr"]) & (dt < max_step_size):
dt = dt * 10
break
else: # radius change is ok
# sanity check: is new planet mass is still greater than the core mass? -> then the planet still has some atmosphere left
# in this case update parameters and go into next RK iteration
if ((M + M_lost) - M_core) >= 0:
M = M + M_lost # new planet mass (M_lost is negative)
t = t_arr[-1] + dt # updated time value t_i_plus_1
M_arr.append(M)
t_arr.append(t)
Lx_arr.append(Lx_evo(t=t, track_dict=track_dict)) # Lx at new time t
# calculate new envelope mass fraction:
M_env = M - M_core # new envelope mass
f_env = (M_env/M)*100 # in %
# calculate new radius with new planet mass/envelope mass fraction & one time step later
R = plmo14.calculate_planet_radius(M_core, f_env, t, flux_at_planet_earth(planet_object.Lbol,
planet_object.distance),
planet_object.metallicity)
R_arr.append(R)
i = i+1 # update step to i+1
break
else:
print("sth went wrong of you see this!")
break
#print("Done!")
return np.array(t_arr), np.array(M_arr), np.array(R_arr), np.array(Lx_arr)
##################################################################################################################################
##################################################################################################################################
......
{
"cells": [
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"a = 5\n",
"b = True"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"step 1\n",
"step 2\n",
"final result\n",
"update planet params & check radius change\n"
]
}
],
"source": [
"for a in [0,5,4,2]:\n",
" while (b==True):\n",
" print(\"step 1\")\n",
" if a == 1:\n",
" print(\"break out of 1\")\n",
" break\n",
" print(\"step 2\")\n",
" if a == 2:\n",
" print(\"break out of 2\")\n",
" break\n",
" print(\"final result\")\n",
" print(\"update planet params & check radius change\")\n",
" b = False\n",
" break"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment