From 0830f9132e53c0e45a413d7d862b98068386b20b Mon Sep 17 00:00:00 2001 From: Laurent Soucasse <135331272+lsoucasse@users.noreply.github.com> Date: Mon, 18 Nov 2024 10:29:19 +0100 Subject: [PATCH 1/6] Implement Rosenbrock method with fixed time step. --- src/mors/rotevo.py | 39 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/src/mors/rotevo.py b/src/mors/rotevo.py index 3868969..bcb7fc8 100644 --- a/src/mors/rotevo.py +++ b/src/mors/rotevo.py @@ -420,6 +420,12 @@ def EvolveRotationStep(Mstar=None,Age=None,OmegaEnv=None,OmegaCore=None,dAgeMax= OmegaCore=OmegaCore , dAge=dAge , dAgeMax=dAgeMax , params=params , StarEvo=StarEvo ) + elif ( params['TimeIntegrationMethod'] == 'RosenbrockFixed' ): + + dAge , dAgeNew , OmegaEnv , OmegaCore = _EvolveRotationStepRBFixed( Mstar=Mstar , Age=Age , OmegaEnv=OmegaEnv , + OmegaCore=OmegaCore , dAge=dAge , dAgeMax=dAgeMax , + params=params , StarEvo=StarEvo ) + else: raise Exception("invalid value of TimeIntegrationMethod in parameters") @@ -528,7 +534,7 @@ def _EvolveRotationStepRKF(Mstar=None,Age=None,OmegaEnv=None,OmegaCore=None,dAge return dAge , dAgeNew , OmegaEnv , OmegaCore def _EvolveRotationStepRB(Mstar=None,Age=None,OmegaEnv=None,OmegaCore=None,dAgeMax=None,dAge=None,params=params.paramsDefault,StarEvo=None): - """Takes basic stellar parameters, evolves by timestep using Runge-Kutta-Fehlberg method.""" + """Takes basic stellar parameters, evolves by timestep using Rosenbrock method.""" # Get coefficients for solver CoefficientsRB = params['CoefficientsRB'] @@ -587,6 +593,37 @@ def _EvolveRotationStepRB(Mstar=None,Age=None,OmegaEnv=None,OmegaCore=None,dAgeM return dAge , dAgeNew , OmegaEnv , OmegaCore +def _EvolveRotationStepRBFixed(Mstar=None,Age=None,OmegaEnv=None,OmegaCore=None,dAgeMax=None,dAge=None,params=params.paramsDefault,StarEvo=None): + """Takes basic stellar parameters, evolves by timestep using Rosenbrock method with fixed time grid.""" + + # Get coefficients for solver + CoefficientsRB = params['CoefficientsRB'] + + # Get time step + dAge = min(0.1 * np.power(Age, 0.75), dAgeMax) + dAgeNew = dAge + + # Setup array to hold integration values + X = np.array([OmegaEnv,OmegaCore]) + nVar = len(X) + + # Get Jacobian and assume it is a constant over timestep + Jac = _JacobianRB( Mstar , Age+dAge , X , nVar , params , StarEvo ) + + # Calculate the k coefficients + kCoeff = _kCoeffRB( Mstar , Age , dAge , X , Jac , nVar , CoefficientsRB , params , StarEvo ) + + # Get updated values + Xnew = copy.deepcopy(X) + for i in range(0,CoefficientsRB['s']): + Xnew += CoefficientsRB['b'][i] * kCoeff[:,i] + + # Save new estimate + OmegaEnv = Xnew[0] + OmegaCore = Xnew[1] + + return dAge , dAgeNew , OmegaEnv , OmegaCore + def _JacobianRB(Mstar,Age,X,nVar,params,StarEvo): """Calculates Jacobian d(dOmega/dt)/dOmega terms for Rosenbrock solver.""" From 0ae676f06538604b86a95263ee808e410b2f8cc0 Mon Sep 17 00:00:00 2001 From: Laurent Soucasse <135331272+lsoucasse@users.noreply.github.com> Date: Mon, 18 Nov 2024 13:17:41 +0100 Subject: [PATCH 2/6] Attempt to test new implementation. --- src/mors/parameters.py | 2 +- tests/{test_spada.py => test_spada_FE.py} | 4 +++- tests/test_spada_RB.py | 23 +++++++++++++++++++++++ 3 files changed, 27 insertions(+), 2 deletions(-) rename tests/{test_spada.py => test_spada_FE.py} (82%) create mode 100644 tests/test_spada_RB.py diff --git a/src/mors/parameters.py b/src/mors/parameters.py index c1d0de0..9bbdafa 100644 --- a/src/mors/parameters.py +++ b/src/mors/parameters.py @@ -18,7 +18,7 @@ def SetDefaultParameters(paramsDefault): """Sets up the default parameters.""" # Time integration solver parameters - paramsDefault['TimeIntegrationMethod'] = 'ForwardEuler' # options are 'ForwardEuler', 'RungeKutta4', 'RungeKuttaFehlberg', 'Rosenbrock' + paramsDefault['TimeIntegrationMethod'] = 'RosenbrockFixed' # options are 'ForwardEuler', 'RungeKutta4', 'RungeKuttaFehlberg', 'Rosenbrock' paramsDefault['nStepMax'] = 10**6 # maximum number of timesteps to allow before error paramsDefault['DeltaDesired'] = 1.0e-5 # desired Delta value to use in numerical solvers such as Runge-Kutta-Fehlberg and Rosenbrock paramsDefault['AgeMinDefault'] = 1.0 # Myr - default age to start evolutionary tracks diff --git a/tests/test_spada.py b/tests/test_spada_FE.py similarity index 82% rename from tests/test_spada.py rename to tests/test_spada_FE.py index 4959c0e..ffce165 100644 --- a/tests/test_spada.py +++ b/tests/test_spada_FE.py @@ -13,7 +13,9 @@ def test_spada(inp, expected): mors.DownloadEvolutionTracks('Spada') - star = mors.Star(Mstar=inp[0], Omega=inp[1]) + params = mors.NewParams() + params['TimeIntegrationMethod'] = 'ForwardEuler' + star = mors.Star(Mstar=inp[0], Omega=inp[1], params = params) ret = ( star.Value(inp[2], 'Rstar'), star.Value(inp[2], 'Lbol'), diff --git a/tests/test_spada_RB.py b/tests/test_spada_RB.py new file mode 100644 index 0000000..04f01cd --- /dev/null +++ b/tests/test_spada_RB.py @@ -0,0 +1,23 @@ +import mors +import pytest +from numpy.testing import assert_allclose + +TEST_DATA = ( +# ((0.128, 45.2, 8.5e1),(0.21263373, 1.68612614e+31, 1.73275380e+28)), +# ((1.113, 17.7, 3.2e3),(1.16581827, 6.81769926e+33, 7.94292986e+28)), +# ((0.995, 1.005, 1.0e4),(1.48904952e+00, 7.80667546e+33, 3.22933432e+28)), + ((1.000, 1.00, 1.0e4),(1.52623102e+00, 8.13153012e+33, 3.23080690e+28)), +) + +@pytest.mark.parametrize("inp,expected", TEST_DATA) +def test_spada(inp, expected): + + mors.DownloadEvolutionTracks('Spada') + star = mors.Star(Mstar=inp[0], Omega=inp[1]) + ret = ( + star.Value(inp[2], 'Rstar'), + star.Value(inp[2], 'Lbol'), + star.Value(inp[2], 'Leuv'), + ) + + assert_allclose(ret, expected, rtol=1e-5, atol=0) From 90a8dcdd6d22d167dbec39da71db52d2ca793283 Mon Sep 17 00:00:00 2001 From: Laurent Soucasse <135331272+lsoucasse@users.noreply.github.com> Date: Mon, 18 Nov 2024 13:56:22 +0100 Subject: [PATCH 3/6] Add test values and reduce tolerance. --- tests/test_spada_FE.py | 2 +- tests/test_spada_RB.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_spada_FE.py b/tests/test_spada_FE.py index ffce165..e26ccc8 100644 --- a/tests/test_spada_FE.py +++ b/tests/test_spada_FE.py @@ -22,4 +22,4 @@ def test_spada(inp, expected): star.Value(inp[2], 'Leuv'), ) - assert_allclose(ret, expected, rtol=1e-5, atol=0) + assert_allclose(ret, expected, rtol=1e-6, atol=0) diff --git a/tests/test_spada_RB.py b/tests/test_spada_RB.py index 04f01cd..ea68a10 100644 --- a/tests/test_spada_RB.py +++ b/tests/test_spada_RB.py @@ -3,9 +3,9 @@ from numpy.testing import assert_allclose TEST_DATA = ( -# ((0.128, 45.2, 8.5e1),(0.21263373, 1.68612614e+31, 1.73275380e+28)), -# ((1.113, 17.7, 3.2e3),(1.16581827, 6.81769926e+33, 7.94292986e+28)), -# ((0.995, 1.005, 1.0e4),(1.48904952e+00, 7.80667546e+33, 3.22933432e+28)), + ((0.128, 45.2, 8.5e1),(0.21264087, 1.68629162e+31, 1.74475018e+28)), + ((1.113, 17.7, 3.2e3),(1.16581971, 6.81767904e+33, 7.94621466e+28)), + ((0.995, 1.005, 1.0e4),(1.48938741e+00, 7.80623475e+33, 3.11812348e+28)), ((1.000, 1.00, 1.0e4),(1.52623102e+00, 8.13153012e+33, 3.23080690e+28)), ) @@ -20,4 +20,4 @@ def test_spada(inp, expected): star.Value(inp[2], 'Leuv'), ) - assert_allclose(ret, expected, rtol=1e-5, atol=0) + assert_allclose(ret, expected, rtol=1e-6, atol=0) From 7648536786d0654898735d797f1854d69b9deef3 Mon Sep 17 00:00:00 2001 From: Laurent Soucasse <135331272+lsoucasse@users.noreply.github.com> Date: Mon, 18 Nov 2024 14:08:40 +0100 Subject: [PATCH 4/6] Cap to a smaller maximum time step of 50 Myr. --- src/mors/parameters.py | 4 ++-- src/mors/rotevo.py | 4 +++- tests/test_spada_RB.py | 4 ++-- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/mors/parameters.py b/src/mors/parameters.py index 9bbdafa..a2001f4 100644 --- a/src/mors/parameters.py +++ b/src/mors/parameters.py @@ -18,7 +18,7 @@ def SetDefaultParameters(paramsDefault): """Sets up the default parameters.""" # Time integration solver parameters - paramsDefault['TimeIntegrationMethod'] = 'RosenbrockFixed' # options are 'ForwardEuler', 'RungeKutta4', 'RungeKuttaFehlberg', 'Rosenbrock' + paramsDefault['TimeIntegrationMethod'] = 'RosenbrockFixed'# options are 'ForwardEuler', 'RungeKutta4', 'RungeKuttaFehlberg', 'Rosenbrock' and 'RosenbrockFixed' paramsDefault['nStepMax'] = 10**6 # maximum number of timesteps to allow before error paramsDefault['DeltaDesired'] = 1.0e-5 # desired Delta value to use in numerical solvers such as Runge-Kutta-Fehlberg and Rosenbrock paramsDefault['AgeMinDefault'] = 1.0 # Myr - default age to start evolutionary tracks @@ -26,7 +26,7 @@ def SetDefaultParameters(paramsDefault): paramsDefault['deltaJac'] = 1.0e-5 # pertubation to use when calculating Jacobian for Rosenbrock solver paramsDefault['CoefficientsRB'] = _CoefficientsRB() # all coefficents needed for the Rosenbrock solver paramsDefault['dAgeMin'] = 1.0e-5 # Myr - minimum timestep to allow - paramsDefault['dAgeMax'] = 5000.0 # Myr - maximum timestep to allow + paramsDefault['dAgeMax'] = 50.0 # Myr - maximum timestep to allow # Physical processes to include paramsDefault['CoreEnvelopeDecoupling'] = True # should core envelope decoupling be included diff --git a/src/mors/rotevo.py b/src/mors/rotevo.py index bcb7fc8..adb5c2c 100644 --- a/src/mors/rotevo.py +++ b/src/mors/rotevo.py @@ -600,7 +600,9 @@ def _EvolveRotationStepRBFixed(Mstar=None,Age=None,OmegaEnv=None,OmegaCore=None, CoefficientsRB = params['CoefficientsRB'] # Get time step - dAge = min(0.1 * np.power(Age, 0.75), dAgeMax) + dAge = 0.1 * np.power(Age, 0.75) + dAge = min( dAge , dAgeMax ) + dAge = min( dAge , params['dAgeMax'] ) dAgeNew = dAge # Setup array to hold integration values diff --git a/tests/test_spada_RB.py b/tests/test_spada_RB.py index ea68a10..ff25adc 100644 --- a/tests/test_spada_RB.py +++ b/tests/test_spada_RB.py @@ -5,8 +5,8 @@ TEST_DATA = ( ((0.128, 45.2, 8.5e1),(0.21264087, 1.68629162e+31, 1.74475018e+28)), ((1.113, 17.7, 3.2e3),(1.16581971, 6.81767904e+33, 7.94621466e+28)), - ((0.995, 1.005, 1.0e4),(1.48938741e+00, 7.80623475e+33, 3.11812348e+28)), - ((1.000, 1.00, 1.0e4),(1.52623102e+00, 8.13153012e+33, 3.23080690e+28)), + ((0.995, 1.005, 1.0e4),(1.48910765e+00, 7.80659378e+33, 3.16581356e+28)), + ((1.000, 1.00, 1.0e4),(1.52588718e+00, 8.13197065e+33, 3.29155751e+28)), ) @pytest.mark.parametrize("inp,expected", TEST_DATA) From 83b169d3ff84475ffb4a62850240e217aec25ce4 Mon Sep 17 00:00:00 2001 From: Laurent Soucasse <135331272+lsoucasse@users.noreply.github.com> Date: Mon, 18 Nov 2024 15:26:31 +0100 Subject: [PATCH 5/6] Update test names. --- tests/test_spada_FE.py | 2 +- tests/test_spada_RB.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_spada_FE.py b/tests/test_spada_FE.py index e26ccc8..0c3c9b8 100644 --- a/tests/test_spada_FE.py +++ b/tests/test_spada_FE.py @@ -10,7 +10,7 @@ ) @pytest.mark.parametrize("inp,expected", TEST_DATA) -def test_spada(inp, expected): +def test_spada_FE(inp, expected): mors.DownloadEvolutionTracks('Spada') params = mors.NewParams() diff --git a/tests/test_spada_RB.py b/tests/test_spada_RB.py index ff25adc..0b460f6 100644 --- a/tests/test_spada_RB.py +++ b/tests/test_spada_RB.py @@ -10,7 +10,7 @@ ) @pytest.mark.parametrize("inp,expected", TEST_DATA) -def test_spada(inp, expected): +def test_spada_RB(inp, expected): mors.DownloadEvolutionTracks('Spada') star = mors.Star(Mstar=inp[0], Omega=inp[1]) From f603efaeb753ec99d6b3ebf8adb798672cfa490e Mon Sep 17 00:00:00 2001 From: Laurent Soucasse <135331272+lsoucasse@users.noreply.github.com> Date: Mon, 18 Nov 2024 15:27:09 +0100 Subject: [PATCH 6/6] Release 24.11.18 --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 84950e8..c262f60 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "fwl-mors" -version = "24.10.27" +version = "24.11.18" description = "Stellar rotation and activity evolution model" readme = "README.md" authors = [ @@ -70,7 +70,7 @@ testpaths = ["tests"] [tool.bumpversion] # https://callowayproject.github.io/bump-my-version/howtos/calver/ -current_version = "24.10.27" +current_version = "24.11.18" parse = """(?x) # Verbose mode (?P # The release part (?:[1-9][0-9])\\. # YY.