Tests using physical simulations

.github/workflows/tests.yaml

@@ -22,76 +22,149 @@ jobs:
       matrix:
         os: ['ubuntu-latest', 'macos-latest']
         python-version: ['3.11', '3.13']
+        include:
+          - os: ubuntu-latest
+            INSTALL_DEPS: sudo apt-get install libnetcdff-dev netcdf-bin
+            CC: gcc
+            CXX: g++
+            FC: gfortran
+          - os: macos-latest
+            INSTALL_DEPS: brew install gfortran netcdf netcdf-fortran
+            CC: gcc
+            CXX: g++
+            FC: gfortran
 
     env:
-      FWL_DATA: $HOME/work/fwl_data
-      RAD_DIR: ./SOCRATES
-      SOCRATES: ./SOCRATES  # https://github.com/FormingWorlds/JANUS/issues/51
+      FWL_DATA: ${{ github.workspace }}/fwl_data
+      RAD_DIR: ${{ github.workspace }}/SOCRATES
 
     runs-on: ${{ matrix.os }}
     steps:
-      - uses: actions/checkout@v4
+
+      # https://stackoverflow.com/a/65356209
+      - name: Install system dependencies
+        run: ${{ matrix.INSTALL_DEPS }}
+
+      # MacOS only: create symbolic link for gfortran
+      - name: Symlink gfortran
+        if: runner.os == 'macOS'
+        run: |
+          sudo ln -s /opt/homebrew/bin/gfortran-13 /opt/homebrew/bin/gfortran
+          sudo ln -s /opt/homebrew/Cellar/gcc/12.*/lib/gcc/12/*.dylib /opt/homebrew/lib/
+          which gfortran
+
+      # Checkout PROTEUS
+      - name: Checkout PROTEUS
+        uses: actions/checkout@v4
+
+      # Checkout SOCRATES
+      - name: Checkout SOCRATES radiative transfer model
+        uses: actions/checkout@v4
+        with:
+          repository: 'nichollsh/SOCRATES'
+          path: ${{ env.RAD_DIR }}
+
+      # Try to restore SOCRATES from the cache
+      - name: Restore SOCRATES from cache
+        uses: actions/cache@v4
+        id: cache-socrates
         with:
-          submodules: 'false'
+          path: |
+            ${{ env.RAD_DIR }}/bin
+            ${{ env.RAD_DIR }}/sbin
+            ${{ env.RAD_DIR }}/set_rad_env
+          key: socrates-${{ hashFiles('SOCRATES/version') }}-${{ runner.os }}
 
+      # Build SOCRATES if we can't restore it from the cache
+      - name: Build SOCRATES
+        if: steps.cache-socrates.outputs.cache-hit != 'true'
+        run: |
+          export LD_LIBRARY_PATH=""
+          echo $RAD_DIR
+          cd $RAD_DIR
+          ./configure
+          ./build_code
+          cd ..
+
+      # Setup Python using the version defined in the matrix
       - name: Set up Python ${{ matrix.python-version }}
         uses: actions/setup-python@v4
         with:
           python-version: ${{ matrix.python-version }}
 
-      - uses: actions/cache@v4
+      # Try to restore the Python environment from the cache
+      - name: Restore Python environment from cache
+        uses: actions/cache@v4
         id: cache-virtualenv
         with:
           path: ${{ env.pythonLocation }}
           key: ${{ env.pythonLocation }}-${{ hashFiles('pyproject.toml') }}
 
-      - name: Install dependencies
-        # if: steps.cache-virtualenv.outputs.cache-hit != 'true'
-        run: |
-          python -m pip install -e .[develop]
+      # Setup Aragog
+      - name: Checkout Aragog interior model
+        uses: actions/checkout@v4
+        with:
+          repository: 'ExPlanetology/aragog'
+          path: 'aragog'
 
-      - uses: actions/cache@v4
+      - name: Install Aragog
+        run: python -m pip install -e aragog/.
+
+      - name: Install Python dependencies
+        run: python -m pip install -e .[develop]
+
+      # Get lookup data for PROTEUS
+      - name: Get FWL data from cache
+        uses: actions/cache@v4
         id: cache-fwl-data
         with:
           path: ${{ env.FWL_DATA }}
           key: fwl-data-2
 
-      - name: Pre-download fwl data
+      # Get FWL data
+      - name: Get additional FWL data
         if: steps.cache-fwl-data.cache-hit != 'true'
         run: |
           proteus get stellar
           proteus get spectral --name Frostflow --bands 48
 
+      # Run PROTEUS tests
       - name: Test with pytest
-        run: |
-          coverage run -m pytest
+        run: coverage run -m pytest
 
+      # Upload output run images if tests fail
       - name: Upload PNG plots as artifacts
-        if: failure()  # Only upload if tests fail
+        if: failure()
         uses: actions/upload-artifact@v3
         with:
           name: png-plots
           path: |
-            output/dummy/plot_elements.png
-            output/dummy/plot_escape.png
-            output/dummy/plot_fluxes_global.png
-            output/dummy/plot_global_lin.png
-            output/dummy/plot_global_log.png
-            output/dummy/plot_observables.png
-            output/dummy/plot_sflux.png
+            output/physical/plot_atmosphere.png
+            output/physical/plot_escape.png
+            output/physical/plot_global_log.png
+            output/physical/plot_observables.png
+            output/physical/plot_stacked.png
+            output/physical/plot_elements.png
+            output/physical/plot_fluxes_atmosphere.png
+            output/physical/plot_interior_cmesh.png
+            output/physical/plot_sflux_cross.png
+            output/physical/plot_emission.png
+            output/physical/plot_interior.png
+            output/physical/plot_sflux.png
 
       - name: Report coverage
         run: |
           coverage json
           export TOTAL=$(python -c "import json;print(json.load(open('coverage.json'))['totals']['percent_covered_display'])")
+          echo "Total coverage: $TOTAL"
           echo "total=$TOTAL" >> $GITHUB_ENV
           echo "### Total coverage: ${TOTAL}%" >> $GITHUB_STEP_SUMMARY
           echo $'\n```' >> $GITHUB_STEP_SUMMARY
           coverage report >> $GITHUB_STEP_SUMMARY
           echo $'\n```' >> $GITHUB_STEP_SUMMARY
 
       - name: Make coverage badge
-        if: ${{ github.ref == 'refs/heads/main' && matrix.python-version == '3.10' }}
+        if: ${{ github.ref == 'refs/heads/main' && matrix.python-version == '3.13' }}
         uses: schneegans/dynamic-badges-action@v1.7.0
         with:
           auth: ${{ secrets.GIST_TOKEN }}

docs/installation.md

@@ -123,6 +123,8 @@ You will need to setup Python (>=3.10) on your system. This can be done via brew
 
 7. Setup numerical computing library (**PETSc**)
 
+    You will need to do this step in an environment with Python <=3.12.
+
     ```console
     ./tools/get_petsc.sh
     ```

input/aragog.toml

@@ -26,7 +26,7 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
     # output files
     [params.out]
         path        = "dummy_aragog"
-        logging     = "INFO"
+        logging     = "DEBUG"
         plot_mod    = 1      # Plotting frequency, 0: wait until completion | n: every n iterations
         plot_fmt    = "png"  # Plotting image file format, "png" or "pdf" recommended
 
@@ -123,7 +123,7 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
 
 # Atmosphere - physics table
 [atmos_clim]
-    prevent_warming = false      # do not allow the planet to heat up
+    prevent_warming = true      # do not allow the planet to heat up
     surface_d       = 0.01      # m, conductive skin thickness
     surface_k       = 2.0       # W m-1 K-1, conductive skin thermal conductivity
     cloud_enabled   = false     # enable water cloud radiative effects

input/hd63433d.toml

@@ -90,7 +90,7 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
     mass    = 0.99      # M_sun
     age_ini = 0.100     # Gyr, model initialisation/start age
 
-    module  = "dummy"
+    module  = "mors"
     [star.mors]
         rot_pctle = 50.0    # rotation percentile
         tracks  = "spada" # evolution tracks: spada | baraffe
@@ -141,10 +141,10 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         p_top           = 1.0e-5        # bar, top of atmosphere grid pressure
         spectral_group  = "Honeyside"   # which gas opacities to include
         spectral_bands  = "256"         # how many spectral bands?
-        num_levels      = 55            # Number of atmospheric grid levels
+        num_levels      = 60            # Number of atmospheric grid levels
         chemistry       = "none"        # "none" | "eq"
         surf_material   = "greybody"    # surface material file for scattering
-        solution_atol   = 2e-3          # solver absolute tolerance
+        solution_atol   = 1e-2          # solver absolute tolerance
         solution_rtol   = 5e-2          # solver relative tolerance
 
     [atmos_clim.janus]
@@ -183,7 +183,7 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         tolerance       = 1.0e-10   # solver tolerance
         tsurf_atol      = 20.0      # tsurf_poststep_change
         tsurf_rtol      = 0.01      # tsurf_poststep_change_frac
-        ini_entropy     = 3200.0    # Surface entropy conditions [J K-1 kg-1]
+        ini_entropy     = 3300.0    # Surface entropy conditions [J K-1 kg-1]
         ini_dsdr        = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]
 
     [interior.aragog]
@@ -192,7 +192,7 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
 
 # Outgassing - physics table
 [outgas]
-    fO2_shift_IW    = 2         # log10(ΔIW), atmosphere/interior boundary oxidation state
+    fO2_shift_IW    = 4         # log10(ΔIW), atmosphere/interior boundary oxidation state
 
     module = "calliope"         # Which outgassing module to use
 

input/trappist1c.toml

@@ -141,10 +141,10 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         p_top           = 1.0e-5        # bar, top of atmosphere grid pressure
         spectral_group  = "Honeyside"   # which gas opacities to include
         spectral_bands  = "256"         # how many spectral bands?
-        num_levels      = 55            # Number of atmospheric grid levels
+        num_levels      = 60            # Number of atmospheric grid levels
         chemistry       = "none"        # "none" | "eq"
         surf_material   = "greybody"    # surface material file for scattering
-        solution_atol   = 2e-3          # solver absolute tolerance
+        solution_atol   = 1e-2          # solver absolute tolerance
         solution_rtol   = 5e-2          # solver relative tolerance
 
     [atmos_clim.janus]
@@ -183,7 +183,7 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
         tolerance       = 1.0e-10   # solver tolerance
         tsurf_atol      = 20.0      # tsurf_poststep_change
         tsurf_rtol      = 0.01      # tsurf_poststep_change_frac
-        ini_entropy     = 3200.0    # Surface entropy conditions [J K-1 kg-1]
+        ini_entropy     = 3300.0    # Surface entropy conditions [J K-1 kg-1]
         ini_dsdr        = -4.698e-6 # Interior entropy gradient [J K-1 kg-1 m-1]
 
     [interior.aragog]
@@ -192,7 +192,7 @@ author  = "Harrison Nicholls, Tim Lichtenberg"
 
 # Outgassing - physics table
 [outgas]
-    fO2_shift_IW    = 2         # log10(ΔIW), atmosphere/interior boundary oxidation state
+    fO2_shift_IW    = 4         # log10(ΔIW), atmosphere/interior boundary oxidation state
 
     module = "calliope"         # Which outgassing module to use
 

pyproject.toml

@@ -49,6 +49,7 @@ dependencies = [
     "fwl-mors>=24.10.27",
     "fwl-calliope>=24.9.10",
     "fwl-zephyrus>=24.10.15",
+    # "aragog",
     "cmcrameri",
     "juliacall",
     "matplotlib",

src/proteus/atmos_clim/agni.py

@@ -72,13 +72,13 @@ def _construct_voldict(hf_row:dict, config:Config, dirs:dict):
 
     # get from hf_row
     vol_dict = {}
+    vol_sum = 0.0
     for vol in gas_list:
-        vmr = hf_row[vol+"_vmr"]
-        if vmr > 1e-40:
-            vol_dict[vol] = vmr
+        vol_dict[vol] = hf_row[vol+"_vmr"]
+        vol_sum += vol_dict[vol]
 
-    # check values
-    if len(vol_dict) == 0:
+    # Check that the total VMR is not zero
+    if vol_sum < 1e-4:
         UpdateStatusfile(dirs, 20)
         raise Exception("All volatiles have a volume mixing ratio of zero")
 
@@ -142,9 +142,6 @@ def init_agni_atmos(dirs:dict, config:Config, hf_row:dict):
         # set all gases as condensates, except the least abundant gas
         condensates = [v[0] for v in vol_sorted[1:]]
 
-        # set top two gases to be condensible
-        # condensates = [v[0] for v in vol_sorted[-2:]]
-
     # Chemistry
     chem_type = config.atmos_clim.agni.chemistry
     include_all = False

src/proteus/atmos_clim/common.py

@@ -13,9 +13,28 @@
 
 log = logging.getLogger("fwl."+__name__)
 
-def read_ncdf_profile(nc_fpath:str):
-    """
-    Read temperature, pressure, height data from NetCDF file.
+def read_ncdf_profile(nc_fpath:str, extra_keys:list=[]):
+    """Read data from atmosphere NetCDF output file.
+
+    Automatically reads pressure (p), temperature (t), height (z) arrays with
+    cell-centre (N) and cell-edge (N+1) values interleaved into a single combined array of
+    length (2*N+1).
+
+    Extra keys can be read-in using the extra_keys parameter. These will be stored with
+    the same dimensions as in the NetCDF file.
+
+    Parameters
+    ----------
+        nc_fpath : str
+            Path to NetCDF file.
+
+        extra_keys : list
+            List of extra keys (strings) to read from the file.
+
+    Returns
+    ----------
+        out : dict
+            Dictionary containing numpy arrays of data from the file.
     """
 
     # open file
@@ -30,14 +49,14 @@ def read_ncdf_profile(nc_fpath:str):
     z = np.array(ds.variables["z"][:])
     zl = np.array(ds.variables["zl"][:])
 
-    nlev = len(p)
+    nlev_c = len(p)
 
     # read pressure, temperature, height data into dictionary values
     out = {}
     out["p"] = [pl[0]]
     out["t"] = [tl[0]]
     out["z"] = [zl[0]]
-    for i in range(nlev):
+    for i in range(nlev_c):
         out["p"].append(p[i])
         out["p"].append(pl[i+1])
 
@@ -47,12 +66,17 @@ def read_ncdf_profile(nc_fpath:str):
         out["z"].append(z[i])
         out["z"].append(zl[i+1])
 
+    # Read extra keys
+    for key in extra_keys:
+        if key in ds.variables.keys():
+            out[key] = np.array(ds.variables[key][:])
+
     # close file
     ds.close()
 
     # convert to np arrays
-    for k in out.keys():
-        out[k] = np.array(out[k], dtype=float)
+    for key in out.keys():
+        out[key] = np.array(out[key], dtype=float)
 
     return out
 

src/proteus/interior/aragog.py

@@ -206,7 +206,7 @@ def GetAragogOutput(hf_row:dict):
     output["M_mantle"] = aragog_output.mantle_mass
     output["T_magma"] = aragog_output.solution_top_temperature
     output["Phi_global"] = float(aragog_output.melt_fraction_global[-1])
-    output["RF_depth"] = aragog_output.rheological_front
+    output["RF_depth"] = float(aragog_output.rheological_front)
     output["F_int"] = aragog_output.convective_heat_flux_basic[-1,-1] # Need to be revised for consistency
 
     output["M_mantle_liquid"] = output["M_mantle"] * output["Phi_global"]