Automated Stellar Evolution Model Calibration Script
This script automatically calibrates stellar evolution models to the Sun using the Dartmouth Stellar Evolution Program (DSEP). It includes functionality for setting up the environment, running simulations with different initial conditions, and analyzing the results to find the best-fit model. It leverages libraries such as numpy, pandas, and scipy for numerical operations, and f90nml for Fortran namelist file manipulation.
import os import numpy as np import argparse import f90nml import shutil import pandas as pd import re from scipy import optimize as op def config_env(args): saveDir = f"output_{args.namelist.split('.')[0]}_test" local = "/home/tboudreaux/d/Astronomy/GraduateSchool/SecondYearProject/DSEP/dsep_Thomas" nml = f"{local}/run/SolarCalibration/{saveDir}" out = f"{local}/run/SolarCalibration/{saveDir}/output" prog = f"{local}/dsepX" zams = f"{local}/zams" prems = f"{local}/prems" opac = f"{local}/opac" atm = f"{local}/surfBC" opal = f"{opac}/OPAL" phx = f"{opac}/phoenix" bcphx = f"{atm}/phoenix/GS98" bckur = f"{atm}/atmk" fort = { 13: f"{local}/run/phys1.low.nml", 14: f"{nml}/", 15: f"{opac}/FERMI.TAB", 38: f"{bckur}/atmk1990p00.tab", 48: f"{opal}/GS98hz", 95: f"{bcphx}/z_p0d0.afe_p0d0.dat", 96: f"{bcphx}/z_m0d5.afe_p0d0.dat", 97: f"{bcphx}/z_m0d7.afe_p0d0.dat", 98: f"{bcphx}/z_m1d0.afe_p0d0.dat", 99: f"{bcphx}/z_m1d5.afe_p0d0.dat", 12: f"{prems}/m100.GS98", 19: f"{out}/", 20: f"{out}/", 37: f"{out}/", 11: f"{out}/", 24: f"{out}/", 25: f"{out}/", 26: f"{out}/" } return saveDir, fort, prog def load_iso(path): header = ["Age", "Log_T", "Log_g", "Log_L", "Log_R", "Y_core", "Z_core", "(Z/X)_surf", "L_H", "L_He", "M_He_core", "M_CO_core"] iso = pd.read_csv(path, delim_whitespace=True, names=header, comment='#') with open(path, 'r') as f: metadata_str = f.readline() Namematch = [x.strip('=') for x in re.findall("[^0-9=.\s]+=", metadata_str)] Valuematch = re.findall("-?\d+(?:\.\d*)?(?:[eE][+\-]?\d+)?", metadata_str) metadata = dict() for i, (name, value) in enumerate(zip(Namematch, Valuematch)): if name != 'COMP': metadata[name] = float(value) return iso, metadata def get_calibration_result(iso, age): return iso.iloc[(iso['Age']-age).abs().idxmin()] def write_namelist(nml, path): nml.write(path) def clear_fort_files(path): for file in os.listdir(path): if os.path.islink(file): if file.startswith('fort'): print(f'deleting: {file}') os.remove(file) def symlink_file(src, dst): if not os.path.exists(src): with open(src, 'w') as f: pass os.symlink(src, dst) def config_fort_files(fort, run): clear_fort_files('.') runFort = fort.copy() runFort[14] += run runFort[19] += run.replace('.nml', '.track') runFort[20] += run.replace('.nml', '.short') runFort[37] += run.replace('.nml', '.iso') runFort[11] += run.replace('.nml', '.last') runFort[24] += run.replace('.nml', '.pmod') runFort[25] += run.replace('.nml', '.penv') runFort[26] += run.replace('.nml', '.patm') for fortKey in fort: print(f'linking {runFort[fortKey]} -> fort.{fortKey}') symlink_file(runFort[fortKey], f'fort.{fortKey}') def run_dsep(ic, fort, nml, saveDir, prog): print(f"Running DSEP with: {ic}") nml['control']['rsclx'] = ic[0] nml['control']['rsclz'] = ic[1] nml['control']['cmixla'] = [ic[2], ic[2]] nml['control']['endage'] = 4.57e9 filename = f"{ic[0]}-{ic[1]}-{ic[2]}.nml" if os.path.exists(f"{saveDir}/{filename}"): os.remove(f"{saveDir}/{filename}") write_namelist(nml, f"{saveDir}/{filename}") runFort = config_fort_files(fort, filename) os.system(f"{prog}/dsepX.master") iso, meta = load_iso('fort.37') cal_results = get_calibration_result(iso, 4.57e9) R = abs(cal_results['Log_R']) L = abs(cal_results['Log_L']) ZX = abs(cal_results['(Z/X)_surf'] - 0.0229) res = np.array([R, L, ZX]) print(f"Completed Run With: \n\tR:{res[0]}\n\tL:{res[1]}\n\tZX:{res[2]}") return res if __name__ == "__main__": parser = argparse.ArgumentParser(description="Automatically Calibrate Stellar Evolution Models to the Sun") parser.add_argument('namelist', type=str, help="namelist file to base calibraton off of") parser.add_argument("-X", type=float, help="Input X for RESCLX") parser.add_argument("-Z", type=float, help="Input Z for RESCLZ") parser.add_argument("-A", type=float, help="Input alpha for CMIXLA") parser.add_argument("-o", "--output", type=str, help="Directory to save DSEP result too", default="output") args = parser.parse_args() saveDir, fort, prog = config_env(args) nml = f90nml.read(args.namelist) if os.path.exists(saveDir): shutil.rmtree(saveDir) os.mkdir(saveDir) os.mkdir(f"{saveDir}/output") res = run_dsep([args.X, args.Z, args.A], fort, nml, saveDir, prog) with open(os.path.join(saveDir, "calibResults.csv"), "w") as f: f.write(f"{args.X},{args.Z},{args.A},{res[0]},{res[1]},{res[2]}\n") # X = np.linspace(0.7816, 0.72, 400) # Z = np.linspace(0.017, 0.018, 400) # A = np.linspace(1.8, 1.9, 400) # # z = Z[0] # a = A[0] # with open("SolarCalibrate-X-Results.csv", 'w') as f: # for x in X: # res = run_dsep([x, z, a], fort, nml, saveDir, prog) # f.write(f"{x},{z},{a},{res[0]},{res[1]},{res[2]}\n") # x = X[0] # a = A[0] # with open("SolarCalibrate-Z-Results.csv", 'w') as f: # for z in Z: # res = run_dsep([x, z, a], fort, nml, saveDir, prog) # f.write(f"{x},{z},{a},{res[0]},{res[1]},{res[2]}\n") # x = X[0] # z = Z[0] # with open("SolarCalibrate-A-Results.csv", 'w') as f: # for a in A: # res = run_dsep([x, z, a], fort, nml, saveDir, prog) # f.write(f"{x},{z},{a},{res[0]},{res[1]},{res[2]}\n")