Automating Solar Model Calibration with PySEP

This Python script is dedicated to automating the calibration of solar models using data from various sources, including atmospheric data, opacity files, and composition data. The script utilizes the PySEP library to calibrate stellar models based on input parameters and output the calibrated models into a specified directory. It includes functionality for handling command-line arguments, reading input files, adjusting model parameters, and performing calibration with error handling.

from pysep.dsep import calibrate
from pysep.dsep import stellarModel as sm
from pysep.dm import filetypes as ft
 
from pysep.io.nml.control.defaults import solar
from pysep.io.nml.physics.defaults import phys1
 
from pysep.opac.tops import open_and_parse
 
from pysep.newpoly.popFileToPolyModel import generate_prems_models
 
import os
 
import pandas as pd
 
import pickle as pkl
 
if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser(description="Run Solar Calibration for models")
    parser.add_argument("inputs", help="path to csv with model file paths listed", type=str)
    parser.add_argument("outputDir", help="path to save calibrated results too", type=str)
    parser.add_argument("-t", "--tau", help="atmospheric fitting point", type=str, default="010")
 
    args = parser.parse_args()
 
    if os.path.exists(args.outputDir):
        raise RuntimeError("Error! Will not run calibration if outputDir exists")
    os.mkdir(args.outputDir)
 
    modelsToBuild = pd.read_csv(os.path.join(args.inputs, "composition.dat"), comment='#')
 
    for idx, modelFiles in modelsToBuild.iterrows():
        atmDirPath = os.path.join(args.inputs, modelFiles['atm'])
        files = os.listdir(atmDirPath)
        atmFilePath = next(filter(lambda x: f"tau{args.tau}.sbc.nc" in x, files))
        atmFile = ft.bcphx95_file(atmFilePath)
 
        opacHFile = ft.opacf_file(os.path.join(args.inputs, modelFiles['opacH']))
 
 
        lowTempOpacFileNames = os.listdir(os.path.join(args.inputs, modelFiles['opacL']))
        lowTempOpacPaths = [os.path.join(args.inputs, modelFiles['opacL'],  x) for x in lowTempOpacFileNames]
 
        parsedAbunFile = open_and_parse(os.path.join(args.inputs, modelFiles['comp']))
 
        controlNML = solar.copy()
        physicsNML = phys1.copy()
 
        controlNML.update_low_temp_opacity_files(lowTempOpacPaths)
 
        premsModelPath = generate_prems_models(parsedAbunFile, [1.0], "premsTMP/", compName=modelFiles['name'], DDAGE=1000, fnameT="temp_m%e%m")
        premsFile = ft.premsf_file(os.path.join('premsTMP/', premsModelPath[1.0]))
 
        modelOutputPath = os.path.join(args.outputDir, modelFiles['name'])
        os.mkdir(modelOutputPath)
 
        model = sm(modelOutputPath, controlNML, physicsNML, opacHFile, premsFile)
        model.update_composition(X=parsedAbunFile['AbundanceRatio']['X'], Z=parsedAbunFile['AbundanceRatio']['Z'])
 
 
 
        try:
            calibResult = calibrate.LRZX_calib_solar(model, 0.7, 0.02, 1.9)
 
            with open(os.path.join(args.outputDir, f"{modelFiles['name']}_calib.pkl"), 'wb') as f:
                pkl.dump(calibResult,f)
        except:
            print(f"Error Calibrating {idx}, moving on!")