Merge pull request #85 from FormingWorlds/article

Plot updates. Import optimisations. Config updates.

input/jgr_grid.cfg

@@ -39,7 +39,7 @@ time_planet                 = 0.        # yr, time since MO start
 time_target                 = 4.567e+9  # yr, target time for MO evolution
 
 # SOCRATES spectral file to use (relative to JANUS folder)
-spectral_file               = JANUS/spectral_files/shared/Dayspring256/Dayspring.sf
+spectral_file               = JANUS/spectral_files/Dayspring/256/Dayspring.sf
 
 # Stellar heating toggle, 0: disabled | 1: enabled
 stellar_heating             = 1
@@ -60,14 +60,14 @@ dt_initial                  = 8e2    # Inital step size
 # Flux convergence scheme and tolerances for surface equilibration
 flux_convergence            = 0         # 0: off | 1: on | 2: restart
 F_atm_bc                    = 0         # Boundary condition choice for F_atm, 0: TOA | 1: Surface
-F_crit                      = 0.1       # Critical flux, below which attempt to stabilise the evolution (0: disabled)
+F_crit                      = 0.0      # Critical flux, below which attempt to stabilise the evolution (0: disabled)
 F_eps                       = 1.0       # W m-2
 F_diff                      = 0.1       # flux fraction
 RF_crit                     = 0.01      # depth fraction
 dTs_atm                     = 30        # K
 skin_d                      = 0.01      # m
 skin_k                      = 2.0       # W m-1 K-1
-prevent_warming             = 0         # Require that the planet only cool down over time, 0: disabled | 1: enabled
+prevent_warming             = 1         # Require that the planet only cool down over time, 0: disabled | 1: enabled
 
 # Break at solidification?
 solid_stop                  = 1
@@ -84,15 +84,15 @@ atmosphere_solve_energy     = 0         # Enable time-stepped atmosphere solutio
 atmosphere_surf_state       = 2         # Atmosphere bottom edge boundary condition, 0: free | 1: fixed at T_surf | 2: conductive skin
 
 # Number of levels
-atmosphere_nlev             = 210
+atmosphere_nlev             = 200
 
 # Temperature limits throughout atmosphere [K]
 min_temperature             = 0.5   
 max_temperature             = 5000.0  
 
 # Clouds
-water_cloud                 = 1         # enable water cloud radiative effects? (1: yes, 0: no)
-alpha_cloud                 = 1.0       # condensate retention fraction (1 -> fully retained)
+water_cloud                 = 0         # enable water cloud radiative effects? (1: yes, 0: no)
+alpha_cloud                 = 0.0       # condensate retention fraction (1 -> fully retained)
 
 # Tropopause type, 0: none | 1: skin temperature | 2: dynamic
 tropopause                  = 1

plot/cpl_atmosphere_cbar.py

@@ -34,7 +34,7 @@ def _read_nc(nc_fpath):
 
 
 # Plotting function
-def plot_atmosphere_cbar(output_dir, plot_both=False):
+def plot_atmosphere_cbar(output_dir):
 
     print("Plot atmosphere temperatures colourbar")
 
@@ -63,15 +63,11 @@ def plot_atmosphere_cbar(output_dir, plot_both=False):
 
     # Initialise plot
     scale = 1.1
-    if plot_both:
-        fig,(ax1,ax2) = plt.subplots(2,1, sharex=True, figsize=(5*scale,8*scale))
-        ax1.set_ylabel("Height [km]")
-    else:
-        fig,ax2 = plt.subplots(1,1,figsize=(5*scale,4*scale))
-    ax2.set_ylabel("Pressure [bar]")
-    ax2.set_xlabel("Temperature [K]")
-    ax2.invert_yaxis()
-    ax2.set_yscale("log")
+    fig,ax = plt.subplots(1,1,figsize=(5*scale,4*scale))
+    ax.set_ylabel("Pressure [bar]")
+    ax.set_xlabel("Temperature [K]")
+    ax.invert_yaxis()
+    ax.set_yscale("log")
 
     # Colour mapping
     norm = mpl.colors.Normalize(vmin=sorted_times[0], vmax=sorted_times[-1])
@@ -82,20 +78,21 @@ def plot_atmosphere_cbar(output_dir, plot_both=False):
     for i in range(nfiles):
         a = 0.7
         c = sm.to_rgba(sorted_times[i])
-        if plot_both:
-            ax1.plot(sorted_t[i], sorted_z[i], color=c, alpha=a, zorder=3)
-        ax2.plot(sorted_t[i], sorted_p[i], color=c, alpha=a, zorder=3)
+        ax.plot(sorted_t[i], sorted_p[i], color=c, alpha=a, zorder=3)
 
     # Grid
-    if plot_both:
-        ax1.grid(alpha=0.2, zorder=2)
-    ax2.grid(alpha=0.2, zorder=2)
-    ax2.set_xlim(0,np.amax(sorted_t)+100)
+    ax.grid(alpha=0.2, zorder=2)
+    ax.set_xlim(0,np.amax(sorted_t)+100)
+    ax.xaxis.set_minor_locator(MultipleLocator(base=250))
+
+    ax.set_ylim(np.amax(sorted_p)*3, np.amin(sorted_p)/3)
+    ax.yaxis.set_major_locator(LogLocator())
 
     # Plot colourbar
-    cax = inset_axes(ax2, width="50%", height="8%", loc='upper right', borderpad=1.0) 
-    cbar = fig.colorbar(sm, cax=cax, orientation='horizontal') 
-    cbar.ax.set_xticks([round(v,1) for v in np.linspace(sorted_times[0] , sorted_times[-1], 4)])
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes('right', size='5%', pad=0.05)
+    cbar = fig.colorbar(sm, cax=cax, orientation='vertical') 
+    cbar.ax.set_yticks([round(v,1) for v in np.linspace(sorted_times[0] , sorted_times[-1], 8)])
     cbar.set_label("Time [Myr]") 
 
     # Save plot

plot/cpl_fluxes.py

@@ -57,7 +57,7 @@ def plot_fluxes_global(output_dir, COUPLER_options, t0=100.0):
     log.info("Plotting global fluxes")
 
     # Get values
-    df = pd.read_csv(output_dir+"/runtime_helpfile.csv", sep="\t")
+    df = pd.read_csv(output_dir+"/runtime_helpfile.csv", sep=r"\s+")
 
     df_atm = df.loc[df['Input']=='Atmosphere'].drop_duplicates(subset=['Time'], keep='last')
     df_atm = df_atm.loc[df_atm['Time']>t0]
@@ -79,7 +79,8 @@ def plot_fluxes_global(output_dir, COUPLER_options, t0=100.0):
 
     # Create plot
     mpl.use('Agg')
-    fig,ax = plt.subplots(figsize=(5,4))
+    scale = 1.1
+    fig,ax = plt.subplots(1,1,figsize=(5*scale,4*scale))
     lw = 2.0
     al = 0.96
 
@@ -100,7 +101,9 @@ def plot_fluxes_global(output_dir, COUPLER_options, t0=100.0):
     ax.set_xscale("log")
     ax.set_xlabel("Time [yr]")
     ax.set_xlim(t_atm[0], t_atm[-1])
-    ax.legend(loc='center left')
+    legend = ax.legend(loc='lower left')
+    legend.get_frame().set_alpha(None)
+    legend.get_frame().set_facecolor((1, 1, 1, 0.99))
     ax.grid(color='black', alpha=0.05)
 
     plt.close()

proteus.py

@@ -67,7 +67,7 @@ def main():
 
                     "steady": 0,           # Number of iterations passed since steady-state declared
                     "steady_loops": 3,     # Number of iterations to perform post-steady state
-                    "steady_check": 25     # Number of iterations to look backwards when checking steady state
+                    "steady_check": 15     # Number of iterations to look backwards when checking steady state
                     }
 
     # Model has completed?

tools/GridPROTEUS.py

@@ -3,7 +3,7 @@
 # Run PROTEUS for a grid of parameters
 
 # Prepare
-import os, itertools, time, subprocess, shutil, sys, multiprocessing, logging
+import os, itertools, time, subprocess, shutil, sys, multiprocessing, logging, gc
 from datetime import datetime
 import numpy as np
 COUPLER_DIR=os.getenv('COUPLER_DIR')
@@ -330,6 +330,7 @@ def run(self,num_threads:int,test_run:bool=False):
                 # Ensure data is written to disk
                 hdl.flush()
                 os.fsync(hdl.fileno()) 
+        gc.collect()
 
 
          # Thread targget
@@ -350,8 +351,8 @@ def _thread_target(cfg_path):
             # Check cfg exists
             cfgexists = False
             waitfor   = 4.0 
-            for _ in range(int(waitfor*10)):  # do n=100*wait_for checks
-                time.sleep(0.1)
+            for _ in range(int(waitfor*100)):  # do n=100*wait_for checks
+                time.sleep(0.01)
                 if os.path.exists(cfgfile):
                     cfgexists = True
                     break
@@ -388,6 +389,7 @@ def _thread_target(cfg_path):
                         status[i] = 2
                         threads[i].join()   # make everything wait until it's completed
                         threads[i].close()  # release resources
+                        gc.collect()
 
             # Print info
             count_que = np.count_nonzero(status == 0)
@@ -461,8 +463,8 @@ def _thread_target(cfg_path):
     # -----
 
     cfg_base = os.path.join(os.getenv('COUPLER_DIR'),"input","jgr_grid.cfg")
-    symlink  = "/network/group/aopp/planetary/RTP035_NICHOLLS_PROTEUS/outputs/jgr_4"
-    pg = Pgrid("jgr_4", cfg_base, symlink_dir=symlink)
+    symlink  = "/network/group/aopp/planetary/RTP035_NICHOLLS_PROTEUS/outputs/CHANGEME"
+    pg = Pgrid("CHANGEME", cfg_base, symlink_dir=symlink)
 
     # pg.add_dimension("Planet")
     # pg.set_dimension_hyper("Planet")

utils/modules_ext.py

@@ -4,7 +4,6 @@
 # https://docs.python.org/3/tutorial/modules.html
 
 # System and file management
-import importlib.util
 import argparse
 import logging
 import pathlib
@@ -27,9 +26,9 @@
 
 import matplotlib.ticker as ticker
 from matplotlib import cm
-import matplotlib.transforms as transforms
 from matplotlib.colors import LinearSegmentedColormap
 from mpl_toolkits.axes_grid1 import make_axes_locatable
+from matplotlib.ticker import LogLocator, LinearLocator, MultipleLocator
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
 import matplotlib.font_manager as fm
 
@@ -42,7 +41,7 @@
 from scipy.interpolate import RectBivariateSpline, interp1d, PchipInterpolator
 from scipy import interpolate
 from scipy.integrate import solve_ivp, odeint
-from scipy.optimize import newton, fsolve, curve_fit
+from scipy.optimize import fsolve
 from scipy import stats
 
 

utils/spider.py

@@ -987,7 +987,7 @@ def _try_spider( time_dict, dirs, COUPLER_options, loop_counter, runtime_helpfil
 
             # Calculate number of macro steps for SPIDER to perform within
             # this time-step of PROTEUS, which sets the number of json files.
-            nsteps = 2
+            nsteps = 1
             dtmacro = math.ceil(dtswitch / nsteps)   # Ensures that dtswitch is divisible by nsteps
             dtswitch = nsteps * dtmacro