From 68d7a446fe23e8d0c0513ccea2db02d1a86ea3e4 Mon Sep 17 00:00:00 2001 From: sciarinl Date: Wed, 18 Dec 2024 10:59:40 +0100 Subject: [PATCH 1/4] Commit updates to the mesa_r15140 interface Update of the mesa_r15140 interface, adding getters required to make it compatible with TRES: - get_core_radius< - get_convective_envelope_mass - get_convective_envelope_radius - get_wind_mass_loss_rate - get_apsidal_motion_constant - get_gyration_radius Update of the evolve_until subroutine in mesa_interface.f90 solving timestep issues occuring in the previous version when the subroutine was called multiple times in a row. --- src/amuse/community/mesa_r15140/interface.f90 | 292 +++++++++++++++--- src/amuse/community/mesa_r15140/interface.py | 121 +++++++- .../community/mesa_r15140/mesa_interface.f90 | 190 +++++++----- 3 files changed, 473 insertions(+), 130 deletions(-) diff --git a/src/amuse/community/mesa_r15140/interface.f90 b/src/amuse/community/mesa_r15140/interface.f90 index 3d95a4115b..9773f39183 100644 --- a/src/amuse/community/mesa_r15140/interface.f90 +++ b/src/amuse/community/mesa_r15140/interface.f90 @@ -3,7 +3,7 @@ module amuse_support character (len=4096) :: AMUSE_inlist_path character (len=4096) :: AMUSE_mesa_dir,AMUSE_mesa_data_dir ! Normally $MESA_DIR and $MESA_DIR/data character (len=4096) :: AMUSE_local_data_dir ! Used for output starting_models - character (len=4096) :: AMUSE_gyre_in_file + character (len=4096) :: AMUSE_gyre_in_file character (len=4096) :: AMUSE_temp_dir ! Used for mesa_temp_caches support double precision :: AMUSE_mass double precision :: AMUSE_metallicity = 0.02d0 @@ -43,7 +43,7 @@ integer function set_MESA_paths(AMUSE_inlist_path_in, & AMUSE_mesa_dir_in, AMUSE_mesa_data_dir_in, & AMUSE_local_data_dir_in, AMUSE_gyre_in_file_in,& AMUSE_temp_dir_in) - + character(*), intent(in) :: AMUSE_inlist_path_in, & AMUSE_mesa_dir_in, AMUSE_local_data_dir_in, AMUSE_gyre_in_file_in, & AMUSE_temp_dir_in, AMUSE_mesa_data_dir_in @@ -60,7 +60,7 @@ end function set_MESA_paths ! Initialize the stellar evolution code integer function initialize_code() - + initialize_code = 0 end function initialize_code @@ -82,11 +82,11 @@ end function cleanup_code integer function new_particle(AMUSE_id, AMUSE_value) integer, intent(out) :: AMUSE_id double precision, intent(in) :: AMUSE_value - + new_particle = new_zams_particle(AMUSE_id, AMUSE_value) end function new_particle -! load zams model +! load zams model integer function new_zams_particle(AMUSE_id, AMUSE_value) integer, intent(out) :: AMUSE_id double precision, intent(in) :: AMUSE_value @@ -115,7 +115,7 @@ integer function new_zams_particle(AMUSE_id, AMUSE_value) end function new_zams_particle -! load an existing stellar model +! load an existing stellar model integer function load_model(AMUSE_id, AMUSE_filename) integer, intent(out) :: AMUSE_id character(len=*), intent(in) :: AMUSE_filename @@ -123,7 +123,7 @@ integer function load_model(AMUSE_id, AMUSE_filename) load_model = -1 - call allocate_star(AMUSE_id, ierr) + call allocate_star(AMUSE_id, ierr) if(ierr/=MESA_SUCESS) return number_of_particles = AMUSE_id @@ -206,7 +206,7 @@ integer function load_photo(AMUSE_id, filename) call allocate_star(id, ierr) if(ierr/=MESA_SUCESS) return - AMUSE_mass = 1.0 + AMUSE_mass = 1.0 AMUSE_id = id number_of_particles = AMUSE_id @@ -222,7 +222,7 @@ integer function load_photo(AMUSE_id, filename) load_photo = 0 end function load_photo - + subroutine set_amuse_options(AMUSE_id) ! Dont call this directly as variables will be reset during initlization ! Instead it must be used as a calback function in finish_init_star() @@ -245,7 +245,7 @@ subroutine set_amuse_options(AMUSE_id) end subroutine set_amuse_options -! Remove a particle +! Remove a particle integer function delete_star(AMUSE_id) integer, intent(in) :: AMUSE_id integer :: ierr @@ -332,12 +332,81 @@ integer function get_core_mass(AMUSE_id, AMUSE_value) endif end function get_core_mass +! Return the current core radius of the star, where hydrogen abundance is <= h1_boundary_limit + integer function get_core_radius(AMUSE_id, AMUSE_value) + integer, intent(in) :: AMUSE_id + double precision, intent(out) :: AMUSE_value + integer :: ierr + + get_core_radius = 0 + call get_history_value(AMUSE_id,'he_core_radius', AMUSE_value, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_core_radius = -1 + endif + end function get_core_radius + +! Return the current size (in mass) of the convective envelope of the star +! It actually returns the mass size of the first convective region from surface to center +! Not appropriate for a star with a radiative envelope. + integer function get_convective_envelope_mass(AMUSE_id, AMUSE_value) + integer, intent(in) :: AMUSE_id + double precision, intent(out) :: AMUSE_value + integer :: ierr + double precision :: conv_env_mtop,conv_env_mbot + + get_convective_envelope_mass = 0 + call get_history_value(AMUSE_id,'conv_mx1_top', conv_env_mtop, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_convective_envelope_mass = -1 + endif + + call get_history_value(AMUSE_id,'conv_mx1_bot', conv_env_mbot, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_convective_envelope_mass = -1 + endif + AMUSE_value = conv_env_mtop - conv_env_mbot + + end function get_convective_envelope_mass + +! Return the current size (in radius) of the convective envelope of the star +! It actually returns the size in radius of the first convective region from surface to center +! Not appropriate for a star with a radiative envelope. + integer function get_convective_envelope_radius(AMUSE_id, AMUSE_value) + integer, intent(in) :: AMUSE_id + double precision, intent(out) :: AMUSE_value + integer :: ierr + double precision :: conv_env_rtop,conv_env_rbot + + get_convective_envelope_radius = 0 + call get_history_value(AMUSE_id,'conv_mx1_top_r', conv_env_rtop, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_convective_envelope_radius = -1 + endif + + call get_history_value(AMUSE_id,'conv_mx1_bot_r', conv_env_rbot, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_convective_envelope_radius = -1 + endif + AMUSE_value = conv_env_rtop - conv_env_rbot + + end function get_convective_envelope_radius + ! Return the current mass loss rate of the star integer function get_mass_loss_rate(AMUSE_id, AMUSE_value) integer, intent(in) :: AMUSE_id double precision, intent(out) :: AMUSE_value integer :: ierr - + get_mass_loss_rate = 0 call get_history_value(AMUSE_id,'star_mdot', AMUSE_value, ierr) @@ -348,6 +417,135 @@ integer function get_mass_loss_rate(AMUSE_id, AMUSE_value) end function get_mass_loss_rate +! Return the current mass loss rate of the star +! Note: this getter is identical to get_mass_loss_rate. +! In the SeBa interface, the corresponding getter is get_wind_mass_loss_rate, and is called +! in triples simulations with TRES (which uses SeBa by default). +! In order to run TRES with MESA, a getter with the same name is required. + integer function get_wind_mass_loss_rate(AMUSE_id, AMUSE_value) + integer, intent(in) :: AMUSE_id + double precision, intent(out) :: AMUSE_value + integer :: ierr + + get_wind_mass_loss_rate = 0 + call get_history_value(AMUSE_id,'star_mdot', AMUSE_value, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_wind_mass_loss_rate = -1 + endif + + end function get_wind_mass_loss_rate + +! Return the current apsidal motion constant k2, assuming a spherical star without accounting for rotation + integer function get_apsidal_motion_constant(AMUSE_id, AMUSE_value) + integer, intent(in) :: AMUSE_id + double precision, intent(out) :: AMUSE_value + integer :: ierr + + get_apsidal_motion_constant = 0 + call get_history_value(AMUSE_id,'apsidal_constant_k2', AMUSE_value, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_apsidal_motion_constant = -1 + endif + + end function get_apsidal_motion_constant + +! Return the current gyration radius + integer function get_gyration_radius(AMUSE_id, AMUSE_value) + integer, intent(in) :: AMUSE_id + double precision, intent(out) :: AMUSE_value + integer :: ierr, n_zone, i + double precision :: I_tot, R_star, M_tot, dm, rmid, mass1, mass2, n_zone_double + character :: rotating_star + ! NOTE (LS 2024): for a rotating star, the moment of inertia can simply be retrieved by + ! calling get_history_value with variable name 'i_rot_total' + ! + ! However, when running non-rotating MESA stars, i_rot_total is not computed in MESA. + ! In the case where the moment of inertia of a non-rotating star is needed, it is + ! directly computed in the interface adding the contribution of spherical layers. + ! In this case, the formula i_rot = 2/3 dm r^2 for each layer is used. + ! For r, we use rmid, the radius value in middle of each shell. + + ! If rotating star, set_initial_surface_rotation_v must be 'T' + call get_star_job_nml(AMUSE_id, 'set_initial_surface_rotation_v', rotating_star, ierr) + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_gyration_radius = -1 + endif + + ! If non-rotating star, calculate moment of inertia manually + if (rotating_star == 'F') then + + call get_history_value(AMUSE_id,'num_zones', n_zone_double, ierr) + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_gyration_radius = -1 + endif + n_zone = int(n_zone_double) + + get_gyration_radius = 0 + I_tot = 0.d0 + + i = 1 + do while (i < n_zone) + call get_profile_value_zone(AMUSE_id, "mass", i, mass1, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_gyration_radius = -1 + endif + + call get_profile_value_zone(AMUSE_id, "mass", i+1, mass2, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_gyration_radius = -1 + endif + + dm = mass1 - mass2 + + call get_profile_value_zone(AMUSE_id, "rmid", i, rmid, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_gyration_radius = -1 + endif + + I_tot = I_tot + (2.d0 * dm * rmid**2.d0)/3.d0 + i = i + 1 + end do + + else ! If rotating star, directly use i_rot_total + call get_history_value(AMUSE_id,'i_rot_total', I_tot, ierr) + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_gyration_radius = -1 + write(*,*)'Issue in accessing i_rot_total' + endif + endif + + call get_history_value(AMUSE_id,'radius', R_star, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_gyration_radius = -1 + endif + + call get_history_value(AMUSE_id,'star_mass', M_tot, ierr) + + if (ierr /= MESA_SUCESS) then + AMUSE_value = -1.0 + get_gyration_radius = -1 + endif + + ! Gyration radius = sqrt(I/MR^2) + AMUSE_value = (I_tot/(M_tot * R_star**2.d0))**(0.5d0) + + end function get_gyration_radius + ! Return the current user-specified mass transfer rate of the star integer function get_manual_mass_transfer_rate(AMUSE_id, AMUSE_value) integer, intent(in) :: AMUSE_id @@ -374,7 +572,7 @@ integer function set_manual_mass_transfer_rate(AMUSE_id, AMUSE_value) character(len=128) :: tmp integer :: ierr set_manual_mass_transfer_rate = 0 - + write(tmp , *) AMUSE_value ierr = set_opt(AMUSE_id, CONTROL_NML,'mass_change', tmp) @@ -598,18 +796,18 @@ integer function get_stellar_type(AMUSE_id, AMUSE_value) call get_history_value(AMUSE_id,'center he3',che3, ierr) if(ierr/=MESA_SUCESS) return call get_history_value(AMUSE_id,'center he4',che4, ierr) - if(ierr/=MESA_SUCESS) return + if(ierr/=MESA_SUCESS) return call get_history_value(AMUSE_id,'center c12',cc12, ierr) - if(ierr/=MESA_SUCESS) return + if(ierr/=MESA_SUCESS) return call get_history_value(AMUSE_id,'center ne20',cne20, ierr) - if(ierr/=MESA_SUCESS) return + if(ierr/=MESA_SUCESS) return call get_history_value(AMUSE_id,'log_LH',lgLH, ierr) - if(ierr/=MESA_SUCESS) return + if(ierr/=MESA_SUCESS) return call get_history_value(AMUSE_id,'log_LHe',lgLHe, ierr) - if(ierr/=MESA_SUCESS) return + if(ierr/=MESA_SUCESS) return call get_history_value(AMUSE_id,'log_L',lgL, ierr) - if(ierr/=MESA_SUCESS) return + if(ierr/=MESA_SUCESS) return call get_history_value(AMUSE_id,'average h1',ah1, ierr) if(ierr/=MESA_SUCESS) return @@ -631,12 +829,12 @@ integer function get_stellar_type(AMUSE_id, AMUSE_value) AMUSE_value = 0 ! Convective low mass star else AMUSE_value = 1 ! Main sequence star - endif - else + endif + else if(lgLHe - lgL < -3) then AMUSE_value = 3 ! Red giant branch else - if(che4 > 1d-4) then + if(che4 > 1d-4) then AMUSE_value = 4 ! Core He burning if(ah1 > 1d-5) then if(ahe3+ahe4 < 0.75*ah1) then @@ -741,7 +939,7 @@ integer function get_number_of_species(AMUSE_id, AMUSE_value) if (ierr /= MESA_SUCESS) then AMUSE_value = -1 get_number_of_species = -1 - endif + endif AMUSE_value = int(val) end function @@ -753,13 +951,13 @@ integer function get_mass_of_species(AMUSE_id, AMUSE_species, AMUSE_value) double precision, intent(out) :: AMUSE_value integer :: ierr get_mass_of_species = 0 - + call get_mass_number_species(AMUSE_id, AMUSE_species, AMUSE_value, ierr) if (ierr /= MESA_SUCESS) then AMUSE_value = -1 get_mass_of_species = ierr - endif + endif end function ! Return the name of chemical abundance given by 'AMUSE_species' of the star @@ -768,7 +966,7 @@ integer function get_name_of_species(AMUSE_id, AMUSE_species, AMUSE_value) character(len=*), intent(out) :: AMUSE_value integer :: ierr - + get_name_of_species = 0 call get_species_name(AMUSE_id, AMUSE_species, AMUSE_value, ierr) @@ -776,18 +974,18 @@ integer function get_name_of_species(AMUSE_id, AMUSE_species, AMUSE_value) if (ierr /= MESA_SUCESS) then AMUSE_value = '' get_name_of_species = ierr - endif + endif end function - - ! Return the chem_id of the species given by AMUSE_species + + ! Return the chem_id of the species given by AMUSE_species integer function get_id_of_species(AMUSE_id, AMUSE_species, AMUSE_value) integer, intent(in) :: AMUSE_id character(len=*), intent(in) :: AMUSE_species integer, intent(out) :: AMUSE_value integer :: ierr - + get_id_of_species = 0 call get_species_id(AMUSE_id, AMUSE_species, AMUSE_value, ierr) @@ -795,7 +993,7 @@ integer function get_id_of_species(AMUSE_id, AMUSE_species, AMUSE_value) if (ierr /= MESA_SUCESS) then AMUSE_value = -1 get_id_of_species= ierr - endif + endif end function @@ -812,7 +1010,7 @@ integer function get_nuclear_network(AMUSE_id, AMUSE_value) if (ierr /= MESA_SUCESS) then AMUSE_value = '' get_nuclear_network= ierr - endif + endif end function get_nuclear_network @@ -828,7 +1026,7 @@ integer function set_nuclear_network(AMUSE_id, AMUSE_value) if (ierr /= MESA_SUCESS) then set_nuclear_network = ierr - endif + endif end function set_nuclear_network @@ -843,7 +1041,7 @@ function evolve_one_step(AMUSE_id) evolve_one_step = -1 call do_evolve_one_step(AMUSE_id, res, ierr) - if (ierr /= MESA_SUCESS) return + if (ierr /= MESA_SUCESS) return evolve_one_step = res @@ -862,7 +1060,7 @@ integer function evolve_for(AMUSE_id, AMUSE_delta_t) evolve_for = -1 call flush() return - endif + endif result = get_age(AMUSE_id, target_times(AMUSE_id)) end function evolve_for @@ -920,7 +1118,7 @@ end function get_maximum_number_of_stars ! This expects all arrays to be in MESA order (1==surface n==center) ! xq = 1-q (where q is mass fraction at zone i) ! star_mass needs to be an array for the interface to work but we only need the total star mass in star_mass(1) - + ! Negative return value indicates an error while positive is the new id integer function new_stellar_model(star_mass, xq, rho, temperature, & XH1,XHe3,XHe4,XC12,XN14,XO16,XNe20,XMg24,XSi28,XS32, & @@ -1107,20 +1305,20 @@ integer function new_stellar_model(star_mass, xq, rho, temperature, & new_model_defined = .true. - contains - + contains + subroutine set_comp(iso,comp) character(len=*) :: iso double precision,dimension(n) :: comp integer :: net_id, ierr - + call get_species_id(id, iso, net_id, ierr) if(ierr/=MESA_SUCESS) then new_stellar_model = ierr return end if xa(net_id,:) = max(1d-50,comp) - + end subroutine set_comp end function new_stellar_model @@ -1301,7 +1499,7 @@ integer function get_gyre(AMUSE_id, mode_l, & call get_gyre_data(AMUSE_ID, mode_l, & add_center_point, keep_surface_point, add_atmosphere, & - fileout, ierr) + fileout, ierr) get_gyre = ierr @@ -1318,7 +1516,7 @@ integer function solve_one_step(AMUSE_ID, first_try, result) integer, intent(in) :: AMUSE_ID logical, intent(in) :: first_try integer, intent(out) :: result - integer :: ierr + integer :: ierr ierr = 0 @@ -1330,7 +1528,7 @@ end function solve_one_step integer function solve_one_step_pre(AMUSE_ID, result) integer, intent(in) :: AMUSE_ID integer, intent(out) :: result - integer :: ierr + integer :: ierr call do_solve_one_step_pre(AMUSE_ID, result, ierr) solve_one_step_pre = ierr @@ -1341,7 +1539,7 @@ end function solve_one_step_pre integer function solve_one_step_post(AMUSE_ID, result) integer, intent(in) :: AMUSE_ID integer, intent(out) :: result - integer :: ierr + integer :: ierr call do_solve_one_step_post(AMUSE_ID, result, ierr) solve_one_step_post = ierr @@ -1355,7 +1553,7 @@ integer function prepare_retry_step(AMUSE_ID, dt, result) integer, intent(in) :: AMUSE_ID double precision,intent(in) :: dt integer, intent(out) :: result - integer :: ierr + integer :: ierr call set_current_dt(AMUSE_ID, dt, ierr) ! Sets dt not dt_next @@ -1372,10 +1570,10 @@ end function prepare_retry_step integer function prepare_redo_step(AMUSE_ID, result) ! Retake the same timestep with the same dt. Useful when mdot changes - ! Does not actualy redo the step + ! Does not actualy redo the step integer, intent(in) :: AMUSE_ID integer, intent(out) :: result - integer :: ierr + integer :: ierr result = do_prepare_for_redo(AMUSE_id) prepare_redo_step = 0 @@ -1383,4 +1581,4 @@ integer function prepare_redo_step(AMUSE_ID, result) end function prepare_redo_step -end module AMUSE_mesa \ No newline at end of file +end module AMUSE_mesa diff --git a/src/amuse/community/mesa_r15140/interface.py b/src/amuse/community/mesa_r15140/interface.py index 980c1b3ef6..b2e15b0001 100644 --- a/src/amuse/community/mesa_r15140/interface.py +++ b/src/amuse/community/mesa_r15140/interface.py @@ -238,14 +238,91 @@ def get_core_mass(): , description="The current core mass of the star, where hydrogen abundance is <= h1_boundary_limit") function.result_type = 'int32' return function + + @legacy_function + def get_core_radius(): + """ + Retrieve the current core radius of the star, where hydrogen abundance is <= h1_boundary_limit + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter('index_of_the_star', dtype='int32', direction=function.IN + , description="The index of the star to get the value of") + function.addParameter('core_radius', dtype='float64', direction=function.OUT + , description="The current core radius of the star, where hydrogen abundance is <= h1_boundary_limit") + function.result_type = 'int32' + return function + @legacy_function + def get_convective_envelope_mass(): + """ + Retrieve the size (in mass) of the convective envelope of the star + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter('index_of_the_star', dtype='int32', direction=function.IN + , description="The index of the star to get the value of") + function.addParameter('convective_envelope_mass', dtype='float64', direction=function.OUT + , description="The current convective envelope mass of the star") + function.result_type = 'int32' + return function + + @legacy_function + def get_convective_envelope_radius(): + """ + Retrieve the size (in radius) of the convective envelope of the star + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter('index_of_the_star', dtype='int32', direction=function.IN + , description="The index of the star to get the value of") + function.addParameter('convective_envelope_radius', dtype='float64', direction=function.OUT + , description="The current convective envelope radius of the star") + function.result_type = 'int32' + return function + @remote_function(can_handle_array=True) def get_mass_loss_rate(index_of_the_star='i'): """ Retrieve the current mass loss rate of the star. (positive for winds, negative for accretion) """ returns (mass_change='d' | units.MSun/units.julianyr) - + + @remote_function(can_handle_array=True) + def get_wind_mass_loss_rate(index_of_the_star='i'): + """ + Retrieve the current mass loss rate of the star. (only winds) + """ + returns (mass_change='d' | units.MSun/units.julianyr) + + @legacy_function + def get_apsidal_motion_constant(): + """ + Retrieve the apsidal motion constant k2, assuming a spherical star without accounting for rotation + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter('index_of_the_star', dtype='int32', direction=function.IN + , description="The index of the star to get the value of") + function.addParameter('apsidal_motion_constant', dtype='float64', direction=function.OUT + , description="The current apsidal motion constant k2 of the star") + function.result_type = 'int32' + return function + + @legacy_function + def get_gyration_radius(): + """ + Retrieve the gyration radius I/MR**2 + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter('index_of_the_star', dtype='int32', direction=function.IN + , description="The index of the star to get the value of") + function.addParameter('gyration_radius', dtype='float64', direction=function.OUT + , description="The current gyration radius of the star") + function.result_type = 'int32' + return function + @remote_function(can_handle_array=True) def get_manual_mass_transfer_rate(index_of_the_star='i'): """ @@ -1306,7 +1383,13 @@ def define_particle_sets(self, handler): handler.add_getter(particle_set_name, 'get_mass', names=('mass',)) handler.add_setter(particle_set_name, 'set_mass', names=('mass',)) handler.add_getter(particle_set_name, 'get_core_mass', names=('core_mass',)) + handler.add_getter(particle_set_name, 'get_core_radius', names=('core_radius',)) + handler.add_getter(particle_set_name, 'get_convective_envelope_mass', names = ('convective_envelope_mass',)) + handler.add_getter(particle_set_name, 'get_convective_envelope_radius', names = ('convective_envelope_radius',)) handler.add_getter(particle_set_name, 'get_mass_loss_rate', names=('wind',)) + handler.add_getter(particle_set_name, 'get_wind_mass_loss_rate', names=('wind_mass_loss_rate',)) + handler.add_getter(particle_set_name, 'get_apsidal_motion_constant', names=('apsidal_motion_constant',)) + handler.add_getter(particle_set_name, 'get_gyration_radius', names=('gyration_radius',)) handler.add_getter(particle_set_name, 'get_age', names=('age',)) handler.add_setter(particle_set_name, 'set_age', names=('age',)) handler.add_getter(particle_set_name, 'get_time_step', names=('time_step',)) @@ -1487,11 +1570,41 @@ def define_methods(self, handler): (handler.INDEX,), (units.MSun, handler.ERROR_CODE,) ) + handler.add_method( + "get_core_radius", + (handler.INDEX,), + (units.RSun, handler.ERROR_CODE,) + ) + handler.add_method( + "get_convective_envelope_mass", + (handler.INDEX,), + (units.MSun, handler.ERROR_CODE,) + ) + handler.add_method( + "get_convective_envelope_radius", + (handler.INDEX,), + (units.RSun, handler.ERROR_CODE,) + ) handler.add_method( "get_mass_loss_rate", (handler.INDEX,), (units.MSun / units.julianyr, handler.ERROR_CODE,) ) + handler.add_method( + "get_wind_mass_loss_rate", + (handler.INDEX,), + (units.MSun / units.julianyr, handler.ERROR_CODE,) + ) + handler.add_method( + "get_apsidal_motion_constant", + (handler.INDEX,), + (handler.NO_UNIT, handler.ERROR_CODE,) + ) + handler.add_method( + "get_gyration_radius", + (handler.INDEX,), + (handler.NO_UNIT, handler.ERROR_CODE,) + ) handler.add_method( "get_manual_mass_transfer_rate", (handler.INDEX,), @@ -1559,6 +1672,12 @@ def define_methods(self, handler): (handler.INDEX, units.julianyr), (handler.ERROR_CODE,) ) + + handler.add_method( + "evolve_one_step", + (handler.INDEX), + (handler.ERROR_CODE,) + ) handler.add_method( "get_max_age_stop_condition", diff --git a/src/amuse/community/mesa_r15140/mesa_interface.f90 b/src/amuse/community/mesa_r15140/mesa_interface.f90 index a3c4d804c7..b823c603e9 100644 --- a/src/amuse/community/mesa_r15140/mesa_interface.f90 +++ b/src/amuse/community/mesa_r15140/mesa_interface.f90 @@ -117,7 +117,7 @@ subroutine load_saved_model(id, filename, ierr) s% job% saved_model_name = trim(filename) s% job% create_pre_main_sequence_model = .false. - + end subroutine load_saved_model @@ -132,7 +132,7 @@ subroutine load_zams_model(id, ierr) s% job% load_saved_model = .false. s% job% create_pre_main_sequence_model = .false. - + end subroutine load_zams_model @@ -167,7 +167,7 @@ subroutine save_mesa_model(id, filename, ierr) integer, intent(out) :: ierr call star_write_model(id, filename, ierr) - + end subroutine save_mesa_model @@ -210,7 +210,7 @@ end subroutine init_callback call star_ptr(id, s, ierr) if (failed('star_ptr',ierr)) return - call init_callback(id) ! Call here and in evolve_controls as there are options that + call init_callback(id) ! Call here and in evolve_controls as there are options that ! need to be set before and after the other inlist options get set id_from_read_star_job = id @@ -232,7 +232,7 @@ subroutine extras_controls(id, ierr) ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return - + ! This is a local copy of extras_controls, ! Plese add your changes to the version in run_star_extras.f90 not this one @@ -242,30 +242,30 @@ subroutine extras_controls(id, ierr) s% job% check_after_step_timing = 0 s% job% time0_initial = 0 s% calculate_Brunt_N2 = .true. - + call init_callback(id) ! Set zbase if its not been set yet - if(s%kap_rq% zbase == -1) s%kap_rq% zbase = s%initial_z + if(s%kap_rq% zbase == -1) s%kap_rq% zbase = s%initial_z ! Override photo_directory otherwise we look in photos/ instead of ./ - ! This does nothing for 15140 but later versions can use this to set + ! This does nothing for 15140 but later versions can use this to set ! the folder for the restart photo to the current directory. if(len_trim(restart_name(id)) > 0) then s%photo_directory = '.' end if call rse_extras_controls(id, ierr) - - + + end subroutine extras_controls - + end subroutine finish_init_star subroutine remove_star(id, ierr) integer, intent(in) :: id - integer, intent(out) :: ierr + integer, intent(out) :: ierr call free_star(id, ierr) if (failed('free_star',ierr)) return @@ -276,7 +276,7 @@ subroutine max_num_stars(num) integer, intent(out) :: num num = max_star_handles - + end subroutine max_num_stars @@ -323,12 +323,12 @@ subroutine set_init_options(id, & !call chdir(output_folder) call star_ptr(id, s, ierr) - if (failed('star_ptr',ierr)) return + if (failed('star_ptr',ierr)) return mesa_dir = mesa_dir_in call const_init(mesa_dir_in, ierr) - if (failed('const_init',ierr)) return + if (failed('const_init',ierr)) return mesa_data_dir = mesa_data_dir_in @@ -350,7 +350,7 @@ subroutine set_init_options(id, & s% terminal_interval = 1 s% write_header_frequency = 100 - mesa_temp_caches_dir = trim(temp_dir) + mesa_temp_caches_dir = trim(temp_dir) s% report_ierr = .false. @@ -364,7 +364,7 @@ end subroutine set_init_options subroutine update_star_job(id, ierr) integer, intent(in) :: id integer, intent(out) :: ierr - type (star_info), pointer :: s + type (star_info), pointer :: s ierr = MESA_SUCESS @@ -382,7 +382,7 @@ subroutine get_net_name(id, net_name, ierr) integer, intent(in) :: id character(len=*), intent(out) :: net_name integer, intent(out) :: ierr - type (star_info), pointer :: s + type (star_info), pointer :: s ierr = MESA_SUCESS @@ -398,7 +398,7 @@ subroutine set_net_name(id, net_name, ierr) integer, intent(in) :: id character(len=*), intent(in) :: net_name integer, intent(out) :: ierr - type (star_info), pointer :: s + type (star_info), pointer :: s ierr = MESA_SUCESS @@ -502,7 +502,7 @@ end subroutine set_kap_nml subroutine do_evolve_one_step(id, result, ierr) integer, intent(in) :: id integer, intent(out) :: result, ierr - type (star_info), pointer :: s + type (star_info), pointer :: s integer :: model_number logical :: continue_evolve_loop, first_try @@ -516,33 +516,33 @@ subroutine do_evolve_one_step(id, result, ierr) call before_step_loop(s% id, ierr) if (failed('before_step_loop',ierr)) return - result = s% extras_start_step(id) - if (result /= keep_going) return + result = s% extras_start_step(id) + if (result /= keep_going) return first_try = .true. - + step_loop: do ! may need to repeat this loop - + if (stop_is_requested(s)) then continue_evolve_loop = .false. result = terminate exit end if - + result = star_evolve_step(id, first_try) if (result == keep_going) result = star_check_model(id) if (result == keep_going) result = s% extras_check_model(id) - if (result == keep_going) result = star_pick_next_timestep(id) + if (result == keep_going) result = star_pick_next_timestep(id) if (result == keep_going) exit step_loop - + model_number = get_model_number(id, ierr) if (failed('get_model_number',ierr)) return - + if (result == retry .and. s% job% report_retries) then write(*,'(i6,3x,a,/)') model_number, & 'retry reason ' // trim(result_reason_str(s% result_reason)) end if - + if (result == redo) then result = star_prepare_to_redo(id) end if @@ -560,7 +560,7 @@ subroutine do_evolve_one_step(id, result, ierr) call after_step_loop(s% id, s% inlist_fname, & dbg, result, ierr) if (failed('after_step_loop',ierr)) return - + call do_saves(id, ierr) if (failed('do_saves',ierr)) return @@ -583,7 +583,7 @@ end function do_prepare_for_retry subroutine do_solve_one_step_pre(id, result, ierr) integer, intent(in) :: id integer, intent(out) :: result, ierr - type (star_info), pointer :: s + type (star_info), pointer :: s integer :: model_number logical :: continue_evolve_loop, first_try @@ -597,8 +597,8 @@ subroutine do_solve_one_step_pre(id, result, ierr) call before_step_loop(s% id, ierr) if (failed('before_step_loop',ierr)) return - result = s% extras_start_step(id) - if (result /= keep_going) return + result = s% extras_start_step(id) + if (result /= keep_going) return end subroutine do_solve_one_step_pre @@ -608,7 +608,7 @@ subroutine do_solve_one_step(id, first_try, result, ierr) integer, intent(in) :: id logical, intent(in) :: first_try integer, intent(out) :: result, ierr - type (star_info), pointer :: s + type (star_info), pointer :: s integer :: model_number logical :: continue_evolve_loop @@ -623,7 +623,7 @@ subroutine do_solve_one_step(id, first_try, result, ierr) result = star_evolve_step(id, first_try) if (result == keep_going) result = star_check_model(id) if (result == keep_going) result = s% extras_check_model(id) - if (result == keep_going) result = star_pick_next_timestep(id) + if (result == keep_going) result = star_pick_next_timestep(id) if (result == keep_going) return end subroutine do_solve_one_step @@ -632,7 +632,7 @@ end subroutine do_solve_one_step subroutine do_solve_one_step_post(id, result, ierr) integer, intent(in) :: id integer, intent(out) :: result, ierr - type (star_info), pointer :: s + type (star_info), pointer :: s integer :: model_number logical :: continue_evolve_loop, first_try @@ -649,7 +649,7 @@ subroutine do_solve_one_step_post(id, result, ierr) call after_step_loop(s% id, s% inlist_fname, & dbg, result, ierr) if (failed('after_step_loop',ierr)) return - + call do_saves(id, ierr) if (failed('do_saves',ierr)) return @@ -662,41 +662,67 @@ subroutine evolve_until(id, delta_t, ierr) integer, intent(in) :: id real(dp), intent(in) :: delta_t integer, intent(out) :: ierr - type (star_info), pointer :: s - integer :: result + type (star_info), pointer :: s + integer :: result, count logical,parameter :: dbg=.false. - real(dp) :: old_age + real(dp) :: end_time, dt_save call star_ptr(id, s, ierr) if (failed('star_ptr',ierr)) return - old_age = s% max_age + ! (LS 2024) + ! Update of evolve_until in case subroutine is called several times subsequently. - s% max_age = s% star_age + delta_t - s% num_adjusted_dt_steps_before_max_age = -1 + ! In the previous version, the timestep was automatically reduced at the end of + ! the evolve loop so that the star age reach max_age. In a subsequent call to + ! evolve_until, the evolution would start again with this reduced timestep, which + ! could in some cases prevent the timestep from increasing. - evolve_loop: do + ! End time of evolution + end_time = s% star_age + delta_t + + ! To check whether the evolve loop was entered + count = 0 + + ! Evolve loop while end_time not reached + ! Use of epsilon(0.0d0) to avoid precision issues + do while (s% star_age + s% dt_next/secyer < end_time*(1.0d0 - epsilon(0.0d0))) call do_evolve_one_step(id, result, ierr) + count = count + 1 if (result /= keep_going) then if (s% result_reason == result_reason_normal) then - exit evolve_loop + exit else ierr = -1 - exit evolve_loop + exit end if end if s% doing_first_model_of_run = .false. - end do evolve_loop + end do - ! In case you want to evolve further - s% max_age = old_age - if(s% dt_next==0) s% dt_next = s% dt + ! Save the next timestep MESA would have required if it was not imposed by end_time + dt_save = s% dt_next - s% doing_first_model_of_run = .false. + ! Set the timestep required to reach end_time + s% dt_next = (end_time - s% star_age)*secyer - end subroutine evolve_until + ! Evolve one last step to reach end_time + call do_evolve_one_step(id, result, ierr) + if (result /= keep_going) then + if (s% result_reason /= result_reason_normal) then + ierr = -1 + end if + end if + ! If several evolve steps were performed, set the timestep for the future evolution + ! to the value it would have reached without the condition to finish at end_time + if (count /= 0) then + s% dt_next = dt_save + endif + s% doing_first_model_of_run = .false. + + end subroutine evolve_until ! *********************************************************************** @@ -743,7 +769,7 @@ subroutine set_current_dt(id, dt, ierr) integer, intent(in) :: id real(dp), intent(in) :: dt integer, intent(out) :: ierr - type (star_info), pointer :: s + type (star_info), pointer :: s call star_ptr(id, s, ierr) if (failed('set_dt',ierr)) return @@ -785,17 +811,17 @@ end subroutine set_new_metallicity subroutine relax_to_new_comp(id, xa, xqs, ierr) integer, intent(in) :: id real(dp), intent(in) :: xa(:,:), xqs(:) - type (star_info), pointer :: s + type (star_info), pointer :: s integer, intent(out) :: ierr integer :: k real(dp) :: max_age - ierr=0 + ierr=0 call star_ptr(id, s, ierr) if (failed('star_ptr',ierr)) return if(s%species /= size(xa,dim=1)) then ierr = -50 - return + return end if call star_relax_composition(id, s% job% num_steps_to_relax_composition, & @@ -808,7 +834,7 @@ end subroutine relax_to_new_comp subroutine relax_to_new_entropy(id, xqs, temperature, rho, ierr) integer, intent(in) :: id real(dp), intent(in) :: xqs(:), temperature(:), rho(:) - type (star_info), pointer :: s + type (star_info), pointer :: s integer, intent(out) :: ierr real(dp), allocatable, dimension(:) :: entropy real(dp), dimension(num_eos_basic_results) :: res,d_dlnd, d_dlnT, d_dabar, d_dzbar @@ -1039,7 +1065,7 @@ subroutine get_next_timestep(id, dt_next, ierr) end subroutine get_next_timestep - + integer function reverse_zone_id(id, zone, ierr) integer, intent(in) :: id, zone integer, intent(out) :: ierr @@ -1063,7 +1089,7 @@ subroutine get_species_name(id, net_id, name, ierr) integer, intent(in) :: id, net_id character(*), intent(out) :: name integer, intent(out) :: ierr - type (star_info), pointer :: s + type (star_info), pointer :: s name='' @@ -1087,7 +1113,7 @@ subroutine get_species_id(id, iso_name, net_id, ierr) character(len=*), intent(in) :: iso_name integer, intent(out) :: ierr type (star_info), pointer :: s - integer :: k + integer :: k call star_ptr(id, s, ierr) if (failed('star_ptr',ierr)) return @@ -1108,8 +1134,8 @@ subroutine get_mass_number_species(id, net_id, mass, ierr) integer, intent(in) :: net_id real(dp), intent(out) :: mass integer, intent(out) :: ierr - type (star_info), pointer :: s - integer :: species_id + type (star_info), pointer :: s + integer :: species_id call star_ptr(id, s, ierr) if (failed('star_ptr',ierr)) return @@ -1129,8 +1155,8 @@ subroutine get_total_mass_species(id, net_id, mass, ierr) integer, intent(in) :: net_id real(dp), intent(out) :: mass integer, intent(out) :: ierr - type (star_info), pointer :: s - integer :: species_id + type (star_info), pointer :: s + integer :: species_id call star_ptr(id, s, ierr) if (failed('star_ptr',ierr)) return @@ -1140,7 +1166,7 @@ subroutine get_total_mass_species(id, net_id, mass, ierr) return end if - mass = dot_product(s% xa(net_id,1:s% nz),s% dq(1:s% nz))/sum(s% dq(1:s% nz)) + mass = dot_product(s% xa(net_id,1:s% nz),s% dq(1:s% nz))/sum(s% dq(1:s% nz)) mass = mass*s% xmstar/Msun end subroutine get_total_mass_species @@ -1150,8 +1176,8 @@ subroutine get_species_at_zone(id, net_id, zone, mass, ierr) integer, intent(in) :: net_id real(dp), intent(out) :: mass integer, intent(out) :: ierr - type (star_info), pointer :: s - integer :: species_id + type (star_info), pointer :: s + integer :: species_id mass = -1 @@ -1161,7 +1187,7 @@ subroutine get_species_at_zone(id, net_id, zone, mass, ierr) if(zone<0 .or. zone> s%nz) then ierr = -2 return - end if + end if if(net_id < 1 .or. net_id > s% species) then ierr =-3 @@ -1176,12 +1202,12 @@ subroutine time_of_step(id, time, ierr) integer, intent(in) :: id integer, intent(out) :: time integer, intent(out) :: ierr - type (star_info), pointer :: s + type (star_info), pointer :: s call star_ptr(id, s, ierr) if (failed('star_ptr',ierr)) return - time = s% job% after_step_timing + time = s% job% after_step_timing end subroutine time_of_step @@ -1192,7 +1218,7 @@ subroutine get_value(id, name, val, k, ierr) character(len=*), intent(in) :: name real(dp), intent(out) :: val integer, intent(out) :: ierr - type (star_info), pointer :: s + type (star_info), pointer :: s ierr = 0 call star_ptr(id, s, ierr) @@ -1204,7 +1230,7 @@ subroutine get_value(id, name, val, k, ierr) end if select case(name) - case('extra_heat') + case('extra_heat') val = s% extra_heat(k) case default ierr = -3 @@ -1219,7 +1245,7 @@ subroutine set_value(id, name, val, k, ierr) character(len=*), intent(in) :: name real(dp), intent(in) :: val integer, intent(out) :: ierr - type (star_info), pointer :: s + type (star_info), pointer :: s ierr = 0 call star_ptr(id, s, ierr) @@ -1231,7 +1257,7 @@ subroutine set_value(id, name, val, k, ierr) end if select case(name) - case('extra_heat') + case('extra_heat') s% extra_heat(k) = val case default ierr = -3 @@ -1256,15 +1282,15 @@ subroutine setup_gyre(gyre_in) call gyre_init(gyre_in) ! Set constants - + call gyre_set_constant('G_GRAVITY', standard_cgrav) call gyre_set_constant('C_LIGHT', clight) call gyre_set_constant('A_RADIATION', crad) - + call gyre_set_constant('M_SUN', Msun) call gyre_set_constant('R_SUN', Rsun) call gyre_set_constant('L_SUN', Lsun) - + call gyre_set_constant('GYRE_DIR', TRIM(mesa_dir)//'/gyre/gyre') @@ -1278,7 +1304,7 @@ subroutine get_gyre_data(id, mode_l, & integer, intent(in) :: id, mode_l integer, intent(out) :: ierr logical, intent(in) :: add_center_point, keep_surface_point, add_atmosphere - type (star_info), pointer :: s + type (star_info), pointer :: s real(dp), allocatable :: global_data(:) real(dp), allocatable :: point_data(:,:) integer :: ipar(1) @@ -1310,18 +1336,18 @@ subroutine get_gyre_data(id, mode_l, & contains subroutine process_mode_ (md, ipar, rpar, retcode) - + type(mode_t), intent(in) :: md integer, intent(inout) :: ipar(:) real(dp), intent(inout) :: rpar(:) integer, intent(out) :: retcode - + integer :: k, unit complex(dp) :: cfreq real(dp) :: freq, growth type(grid_t) :: gr - - retcode= 0 + + retcode= 0 ipar(1) = 0 cfreq = md% freq('HZ') gr = md%grid() From df19ad02a40c8a96f0c45a4ccc7f9fee6f5c6baf Mon Sep 17 00:00:00 2001 From: sciarinl Date: Wed, 18 Dec 2024 14:08:35 +0100 Subject: [PATCH 2/4] Commit adding test scripts Two test scripts are provided. There are located in src/amuse/community/mesa_r15140/Tests test_new_getters.py tests the functionality of the getters added to the interace. test_updated_evolve_until.py tests the functionality of the updated evolve_until function in mesa_interface.f90. --- .../mesa_r15140/Tests/test_new_getters.py | 98 +++++++++++++++++++ .../Tests/test_updated_evolve_until.py | 85 ++++++++++++++++ 2 files changed, 183 insertions(+) create mode 100644 src/amuse/community/mesa_r15140/Tests/test_new_getters.py create mode 100644 src/amuse/community/mesa_r15140/Tests/test_updated_evolve_until.py diff --git a/src/amuse/community/mesa_r15140/Tests/test_new_getters.py b/src/amuse/community/mesa_r15140/Tests/test_new_getters.py new file mode 100644 index 0000000000..3e8d0689ff --- /dev/null +++ b/src/amuse/community/mesa_r15140/Tests/test_new_getters.py @@ -0,0 +1,98 @@ +from amuse.lab import * +import numpy as np +import matplotlib.pyplot as plt + + +# This script evolves a 10 Msun single star with MESA +# To test the getters added to the mesa interface + +stars = Particles(1) + +stars[0].mass = 10 | units.MSun +stars[0].metallicity = 0.02 + + +stellar_code = MESA(redirection="none") + +evo_stars = stellar_code.particles.add_particles(stars) + +# Standard MESA setup +stellar_code.particles[0].set_control('initial_z',0.02) +stellar_code.particles[0].set_kap('use_Type2_opacities',True) +stellar_code.particles[0].set_kap('Zbase', 0.02) + +# Add wind mass loss +stellar_code.particles[0].set_control('hot_wind_scheme','Dutch') +stellar_code.particles[0].set_control('Dutch_scaling_factor',1.0) + +# To store data +times = np.array([]) + +core_radius = np.array([]) +conv_envelope_mass = np.array([]) +conv_envelope_radius = np.array([]) +wind_mass_loss_rate = np.array([]) +apsidal_motion_constant = np.array([]) +gyration_radius = np.array([]) + + + + +time = 0 | units.Myr + + +# Evolve and store the data retrieved from the getters +while time < 0.1 | units.Myr: + stellar_code.evolve_for(1,20|units.kyr) + time = (stellar_code.particles[0].age) + + core_radius_temp = stellar_code.particles[0].core_radius + conv_envelope_mass_temp = stellar_code.particles[0].convective_envelope_mass + conv_envelope_radius_temp = stellar_code.particles[0].convective_envelope_radius + wind_mass_loss_rate_temp = stellar_code.particles[0].wind_mass_loss_rate + apsidal_motion_constant_temp = stellar_code.particles[0].apsidal_motion_constant + gyration_radius_temp = stellar_code.particles[0].gyration_radius + + + times = np.append(times,time.value_in(units.kyr)) + + core_radius = np.append(core_radius,core_radius_temp.value_in(units.RSun)) + conv_envelope_mass = np.append(conv_envelope_mass,conv_envelope_mass_temp.value_in(units.MSun)) + conv_envelope_radius = np.append(conv_envelope_radius,conv_envelope_radius_temp.value_in(units.RSun)) + wind_mass_loss_rate = np.append(wind_mass_loss_rate,wind_mass_loss_rate_temp.value_in(units.MSun / units.yr)) + apsidal_motion_constant = np.append(apsidal_motion_constant,apsidal_motion_constant_temp) + gyration_radius = np.append(gyration_radius,gyration_radius_temp) + + + +# After the evolution, plot the data as function of time + +plt.plot(times,core_radius,linewidth=2) +plt.xlabel(r'$t$ [kyr]',fontsize=13) +plt.ylabel(r'$R_{core}$ [R$_\odot$]',fontsize=13) +plt.show() + +plt.plot(times,conv_envelope_mass,linewidth=2) +plt.xlabel(r'$t$ [kyr]',fontsize=13) +plt.ylabel(r'$M_{env}$ [M$_\odot$]',fontsize=13) +plt.show() + +plt.plot(times,conv_envelope_radius,linewidth=2) +plt.xlabel(r'$t$ [kyr]',fontsize=13) +plt.ylabel(r'$R_{env}$ [R$_\odot$]',fontsize=13) +plt.show() + +plt.plot(times,wind_mass_loss_rate,linewidth=2) +plt.xlabel(r'$t$ [kyr]',fontsize=13) +plt.ylabel(r'$\dot M$ [M$_\odot$/yr]',fontsize=13) +plt.show() + +plt.plot(times,apsidal_motion_constant,linewidth=2) +plt.xlabel(r'$t$ [kyr]',fontsize=13) +plt.ylabel(r'$k_2$',fontsize=13) +plt.show() + +plt.plot(times,gyration_radius,linewidth=2) +plt.xlabel(r'$t$ [kyr]',fontsize=13) +plt.ylabel(r'$r_g$',fontsize=13) +plt.show() \ No newline at end of file diff --git a/src/amuse/community/mesa_r15140/Tests/test_updated_evolve_until.py b/src/amuse/community/mesa_r15140/Tests/test_updated_evolve_until.py new file mode 100644 index 0000000000..7b7ac5f397 --- /dev/null +++ b/src/amuse/community/mesa_r15140/Tests/test_updated_evolve_until.py @@ -0,0 +1,85 @@ +from amuse.lab import * +import numpy as np +import matplotlib.pyplot as plt + + +# This script evolves a 10 Msun single star with MESA +# To test the updated evolve function (evolve_until in mesa_interface.f90) + +stars = Particles(1) + +stars[0].mass = 10 | units.MSun +stars[0].metallicity = 0.02 + + +stellar_code = MESA(redirection="none") + +evo_stars = stellar_code.particles.add_particles(stars) + +# Standard MESA setup +stellar_code.particles[0].set_control('initial_z',0.02) +stellar_code.particles[0].set_kap('use_Type2_opacities',True) +stellar_code.particles[0].set_kap('Zbase', 0.02) + + +# To store time and timestep +times = np.array([]) +time_step = np.array([]) + + +time = 0 | units.Myr + + +# Subsequent calls to evolve function (the updated function is evolve_until in mesa_interface.f90) +# To check that the timestep can increase +# In the previous version it would be prevented from increasing when the evolve function was called +# multiple times in a row + +# With the update, it increases freely up to the maximum allowed value, i.e. the value given as argument +# to evolve_for + + +while time < 0.02 | units.Myr: + time = (stellar_code.particles[0].age) + dt = stellar_code.particles[0].time_step + + times = np.append(times,time.value_in(units.kyr)) + time_step = np.append(time_step,dt.value_in(units.kyr)) + + stellar_code.evolve_for(1,1|units.kyr) + + + +while time < 0.06 | units.Myr: + time = (stellar_code.particles[0].age) + dt = stellar_code.particles[0].time_step + + times = np.append(times,time.value_in(units.kyr)) + time_step = np.append(time_step,dt.value_in(units.kyr)) + + stellar_code.evolve_for(1,5|units.kyr) + + +while time < 0.12 | units.Myr: + time = (stellar_code.particles[0].age) + dt = stellar_code.particles[0].time_step + + times = np.append(times,time.value_in(units.kyr)) + time_step = np.append(time_step,dt.value_in(units.kyr)) + + stellar_code.evolve_for(1,10|units.kyr) + + + +# After the evolution, the timesteps are plotted as function of time. The three visible plateau +# correspond to the value given as argument to evolve_for (1 kyr, 5 kyrs, 10 kyrs), +# and indicate that the timestep was able to reach the maximum allowed value. + +# When running the same script with the previous evolve_until subroutine, the time step is cut at +# the end of each call to evolve_until and can not increase freely to reach the maximum allowed value, +# which is why the plateau are not seen in this case. + +plt.plot(times,time_step,linewidth=2) +plt.xlabel(r'$t$ [kyr]',fontsize=13) +plt.ylabel(r'd$t$ [kyr]',fontsize=13) +plt.show() From 386d5794735edca923dae293fdcbcdb2daa343b1 Mon Sep 17 00:00:00 2001 From: sciarinl Date: Thu, 16 Jan 2025 11:53:49 +0100 Subject: [PATCH 3/4] Update the test scripts Commit to update the test scripts. The test scripts added in the PR, which were located in src/amuse/community/mesa_r15140/Tests have been removed (they were not in the intended format and were plotting their results). The script test_mesa_15140.py located in src/amuse/test/suite/codes_tests has been extended to include additional tests (test26 and test27) to test the updated evolve_until (in test26) and the functionality of the added getters (in test27). --- .../mesa_r15140/Tests/test_new_getters.py | 98 ------------------- .../Tests/test_updated_evolve_until.py | 85 ---------------- .../test/suite/codes_tests/test_mesa_15140.py | 90 +++++++++++++++++ 3 files changed, 90 insertions(+), 183 deletions(-) delete mode 100644 src/amuse/community/mesa_r15140/Tests/test_new_getters.py delete mode 100644 src/amuse/community/mesa_r15140/Tests/test_updated_evolve_until.py diff --git a/src/amuse/community/mesa_r15140/Tests/test_new_getters.py b/src/amuse/community/mesa_r15140/Tests/test_new_getters.py deleted file mode 100644 index 3e8d0689ff..0000000000 --- a/src/amuse/community/mesa_r15140/Tests/test_new_getters.py +++ /dev/null @@ -1,98 +0,0 @@ -from amuse.lab import * -import numpy as np -import matplotlib.pyplot as plt - - -# This script evolves a 10 Msun single star with MESA -# To test the getters added to the mesa interface - -stars = Particles(1) - -stars[0].mass = 10 | units.MSun -stars[0].metallicity = 0.02 - - -stellar_code = MESA(redirection="none") - -evo_stars = stellar_code.particles.add_particles(stars) - -# Standard MESA setup -stellar_code.particles[0].set_control('initial_z',0.02) -stellar_code.particles[0].set_kap('use_Type2_opacities',True) -stellar_code.particles[0].set_kap('Zbase', 0.02) - -# Add wind mass loss -stellar_code.particles[0].set_control('hot_wind_scheme','Dutch') -stellar_code.particles[0].set_control('Dutch_scaling_factor',1.0) - -# To store data -times = np.array([]) - -core_radius = np.array([]) -conv_envelope_mass = np.array([]) -conv_envelope_radius = np.array([]) -wind_mass_loss_rate = np.array([]) -apsidal_motion_constant = np.array([]) -gyration_radius = np.array([]) - - - - -time = 0 | units.Myr - - -# Evolve and store the data retrieved from the getters -while time < 0.1 | units.Myr: - stellar_code.evolve_for(1,20|units.kyr) - time = (stellar_code.particles[0].age) - - core_radius_temp = stellar_code.particles[0].core_radius - conv_envelope_mass_temp = stellar_code.particles[0].convective_envelope_mass - conv_envelope_radius_temp = stellar_code.particles[0].convective_envelope_radius - wind_mass_loss_rate_temp = stellar_code.particles[0].wind_mass_loss_rate - apsidal_motion_constant_temp = stellar_code.particles[0].apsidal_motion_constant - gyration_radius_temp = stellar_code.particles[0].gyration_radius - - - times = np.append(times,time.value_in(units.kyr)) - - core_radius = np.append(core_radius,core_radius_temp.value_in(units.RSun)) - conv_envelope_mass = np.append(conv_envelope_mass,conv_envelope_mass_temp.value_in(units.MSun)) - conv_envelope_radius = np.append(conv_envelope_radius,conv_envelope_radius_temp.value_in(units.RSun)) - wind_mass_loss_rate = np.append(wind_mass_loss_rate,wind_mass_loss_rate_temp.value_in(units.MSun / units.yr)) - apsidal_motion_constant = np.append(apsidal_motion_constant,apsidal_motion_constant_temp) - gyration_radius = np.append(gyration_radius,gyration_radius_temp) - - - -# After the evolution, plot the data as function of time - -plt.plot(times,core_radius,linewidth=2) -plt.xlabel(r'$t$ [kyr]',fontsize=13) -plt.ylabel(r'$R_{core}$ [R$_\odot$]',fontsize=13) -plt.show() - -plt.plot(times,conv_envelope_mass,linewidth=2) -plt.xlabel(r'$t$ [kyr]',fontsize=13) -plt.ylabel(r'$M_{env}$ [M$_\odot$]',fontsize=13) -plt.show() - -plt.plot(times,conv_envelope_radius,linewidth=2) -plt.xlabel(r'$t$ [kyr]',fontsize=13) -plt.ylabel(r'$R_{env}$ [R$_\odot$]',fontsize=13) -plt.show() - -plt.plot(times,wind_mass_loss_rate,linewidth=2) -plt.xlabel(r'$t$ [kyr]',fontsize=13) -plt.ylabel(r'$\dot M$ [M$_\odot$/yr]',fontsize=13) -plt.show() - -plt.plot(times,apsidal_motion_constant,linewidth=2) -plt.xlabel(r'$t$ [kyr]',fontsize=13) -plt.ylabel(r'$k_2$',fontsize=13) -plt.show() - -plt.plot(times,gyration_radius,linewidth=2) -plt.xlabel(r'$t$ [kyr]',fontsize=13) -plt.ylabel(r'$r_g$',fontsize=13) -plt.show() \ No newline at end of file diff --git a/src/amuse/community/mesa_r15140/Tests/test_updated_evolve_until.py b/src/amuse/community/mesa_r15140/Tests/test_updated_evolve_until.py deleted file mode 100644 index 7b7ac5f397..0000000000 --- a/src/amuse/community/mesa_r15140/Tests/test_updated_evolve_until.py +++ /dev/null @@ -1,85 +0,0 @@ -from amuse.lab import * -import numpy as np -import matplotlib.pyplot as plt - - -# This script evolves a 10 Msun single star with MESA -# To test the updated evolve function (evolve_until in mesa_interface.f90) - -stars = Particles(1) - -stars[0].mass = 10 | units.MSun -stars[0].metallicity = 0.02 - - -stellar_code = MESA(redirection="none") - -evo_stars = stellar_code.particles.add_particles(stars) - -# Standard MESA setup -stellar_code.particles[0].set_control('initial_z',0.02) -stellar_code.particles[0].set_kap('use_Type2_opacities',True) -stellar_code.particles[0].set_kap('Zbase', 0.02) - - -# To store time and timestep -times = np.array([]) -time_step = np.array([]) - - -time = 0 | units.Myr - - -# Subsequent calls to evolve function (the updated function is evolve_until in mesa_interface.f90) -# To check that the timestep can increase -# In the previous version it would be prevented from increasing when the evolve function was called -# multiple times in a row - -# With the update, it increases freely up to the maximum allowed value, i.e. the value given as argument -# to evolve_for - - -while time < 0.02 | units.Myr: - time = (stellar_code.particles[0].age) - dt = stellar_code.particles[0].time_step - - times = np.append(times,time.value_in(units.kyr)) - time_step = np.append(time_step,dt.value_in(units.kyr)) - - stellar_code.evolve_for(1,1|units.kyr) - - - -while time < 0.06 | units.Myr: - time = (stellar_code.particles[0].age) - dt = stellar_code.particles[0].time_step - - times = np.append(times,time.value_in(units.kyr)) - time_step = np.append(time_step,dt.value_in(units.kyr)) - - stellar_code.evolve_for(1,5|units.kyr) - - -while time < 0.12 | units.Myr: - time = (stellar_code.particles[0].age) - dt = stellar_code.particles[0].time_step - - times = np.append(times,time.value_in(units.kyr)) - time_step = np.append(time_step,dt.value_in(units.kyr)) - - stellar_code.evolve_for(1,10|units.kyr) - - - -# After the evolution, the timesteps are plotted as function of time. The three visible plateau -# correspond to the value given as argument to evolve_for (1 kyr, 5 kyrs, 10 kyrs), -# and indicate that the timestep was able to reach the maximum allowed value. - -# When running the same script with the previous evolve_until subroutine, the time step is cut at -# the end of each call to evolve_until and can not increase freely to reach the maximum allowed value, -# which is why the plateau are not seen in this case. - -plt.plot(times,time_step,linewidth=2) -plt.xlabel(r'$t$ [kyr]',fontsize=13) -plt.ylabel(r'd$t$ [kyr]',fontsize=13) -plt.show() diff --git a/src/amuse/test/suite/codes_tests/test_mesa_15140.py b/src/amuse/test/suite/codes_tests/test_mesa_15140.py index eb07c96796..f2cd4e964b 100644 --- a/src/amuse/test/suite/codes_tests/test_mesa_15140.py +++ b/src/amuse/test/suite/codes_tests/test_mesa_15140.py @@ -1377,3 +1377,93 @@ def test25(self): ) self.assertAlmostEqual(composition[:, -1], [0.75, 0, 0, 0, 0, 0, 0, 0.25]) instance.stop() + + def test26(self): + print("Testing basic operations: evolve...") + instance = self.new_instance_of_an_optional_code(MESAInterface) + if instance is None: + print("MESA was not built. Skipping test.") + return + set_mesa_paths_instance(instance) + instance.initialize_code() + (index_of_the_star, error) = instance.new_particle(1.0) + self.assertEqual(0, error) + self.assertEqual(index_of_the_star, 1) + self.assertEqual(0, instance.commit_particles()) + + + # Subsequent calls to evolve function (the updated function is evolve_until in mesa_interface.f90) + # To check that the timestep can increase + # In the previous version it would be prevented from increasing when the evolve function was called + # multiple times in a row + + for i in range(10): + self.assertEqual(0,instance.evolve_for(index_of_the_star, 1.0e6)) + + # With the update, it increases freely up to the maximum allowed value, i.e. the value given as argument + # to evolve_for. In this example, evolve_for is called several time in a row with value 1 Myr + # After a few (6) calls to evolve_for, the MESA timestep has reached the maximum value it can, i.e. 1 Myr + # After the last call, the next timestep should be 1.2*previous_timestep, i.e. 1.2 Myr, if the timestep can + # increase freely. With the previous version of evolve_until, next_timestep would not be 1.2 Myr. + + desired_next_timestep_yr = 1.2e6 #value in yr + desired_next_timestep_second = desired_next_timestep_yr * 365.25*24*3600 #value in seconds + + obtained_next_time_step = instance.get_time_step(index_of_the_star).values()[0] + + #check desired_next_timestep has the intended value (1.2 Myr) + self.assertEqual(desired_next_timestep_second, obtained_next_time_step) + + instance.stop() + + def test27(self): + print("Testing basic operations: evolve...") + instance = self.new_instance_of_an_optional_code(MESAInterface) + if instance is None: + print("MESA was not built. Skipping test.") + return + set_mesa_paths_instance(instance) + instance.initialize_code() + (index_of_the_star, error) = instance.new_particle(1.0) + self.assertEqual(0, error) + self.assertEqual(index_of_the_star, 1) + self.assertEqual(0, instance.commit_particles()) + + initial_dt = 1.0e5 + + instance.set_time_step(index_of_the_star, initial_dt) + self.assertEqual([initial_dt, 0], list(instance.get_time_step(index_of_the_star).values())) + + self.assertEqual(0, instance.evolve_for(index_of_the_star, initial_dt)) + self.assertEqual( + [initial_dt, 0], + list(instance.get_age(index_of_the_star).values()) + ) + + target_end_time = 3.0e5 # (years) + self.assertEqual( + 0, instance.evolve_for( + index_of_the_star, target_end_time-initial_dt + ) + ) + + self.assertTrue( + instance.get_age(index_of_the_star)['age'] >= target_end_time + ) + #Tests the functionality of new getters + (ML_of_the_star, error) = instance.get_wind_mass_loss_rate(index_of_the_star) + self.assertEqual(0, error) + self.assertAlmostEqual(ML_of_the_star, 0, -4) + (gyr_of_the_star, error) = instance.get_gyration_radius(index_of_the_star) + self.assertEqual(0, error) + self.assertAlmostEqual(gyr_of_the_star, 0.3, 1) + (k2_of_the_star, error) = instance.get_apsidal_motion_constant(index_of_the_star) + self.assertEqual(0, error) + self.assertAlmostEqual(k2_of_the_star, 0.024, 2) + (Renv_of_the_star, error) = instance.get_convective_envelope_radius(index_of_the_star) + self.assertEqual(0, error) + self.assertAlmostEqual(Renv_of_the_star, 0.12, 2) + (Menv_of_the_star, error) = instance.get_convective_envelope_mass(index_of_the_star) + self.assertEqual(0, error) + self.assertAlmostEqual(Menv_of_the_star, 0.08, 2) + instance.stop() From 50711778eed82e9a04c6d03f73c71bbeec7fa7c0 Mon Sep 17 00:00:00 2001 From: sciarinl Date: Fri, 17 Jan 2025 20:47:27 +0100 Subject: [PATCH 4/4] Update of evolve_until in mesa_interface.f90 Commit to update the evolve_until improving the compactness of the code and adding a few comments. --- .../community/mesa_r15140/mesa_interface.f90 | 48 ++++++++++--------- 1 file changed, 25 insertions(+), 23 deletions(-) diff --git a/src/amuse/community/mesa_r15140/mesa_interface.f90 b/src/amuse/community/mesa_r15140/mesa_interface.f90 index b823c603e9..8344e5a0f6 100644 --- a/src/amuse/community/mesa_r15140/mesa_interface.f90 +++ b/src/amuse/community/mesa_r15140/mesa_interface.f90 @@ -663,7 +663,7 @@ subroutine evolve_until(id, delta_t, ierr) real(dp), intent(in) :: delta_t integer, intent(out) :: ierr type (star_info), pointer :: s - integer :: result, count + integer :: result, count_calls_to_evolve_one_step logical,parameter :: dbg=.false. real(dp) :: end_time, dt_save @@ -681,14 +681,28 @@ subroutine evolve_until(id, delta_t, ierr) ! End time of evolution end_time = s% star_age + delta_t - ! To check whether the evolve loop was entered - count = 0 + ! Counts the number of calls to evolve_one_step + count_calls_to_evolve_one_step = 0 + + ! If MESA delta_t (dt_next) >= AMUSE delta_t (delta_t), only one call to + ! evolve_one_step is needed (setting MESA_delta_t to AMUSE_delta_t). In this case + ! after the loop count_calls_to_evolve_one_step = 1 + ! If MESA delta_t < AMUSE delta_t, several calls to evolve_one_step. At last step + ! MESA delta_t is set so that star_age reaches end_time. In this case after the loop + ! count_calls_to_evolve_one_step > 1 and MESA time_step is set to the value it was + ! before being set to reach end_time. ! Evolve loop while end_time not reached ! Use of epsilon(0.0d0) to avoid precision issues - do while (s% star_age + s% dt_next/secyer < end_time*(1.0d0 - epsilon(0.0d0))) + do while (s% star_age < end_time*(1.0d0 - epsilon(0.0d0))) + if (s% star_age + s% dt_next/secyer > end_time*(1.0d0 - epsilon(0.0d0))) then + ! Last step: save the next timestep MESA would have required if it was not imposed by end_time + dt_save = s% dt_next + + ! Set the last timestep so that star_age reaches end_time + s% dt_next = (end_time - s% star_age)*secyer + end if call do_evolve_one_step(id, result, ierr) - count = count + 1 if (result /= keep_going) then if (s% result_reason == result_reason_normal) then exit @@ -698,27 +712,15 @@ subroutine evolve_until(id, delta_t, ierr) end if end if s% doing_first_model_of_run = .false. + count_calls_to_evolve_one_step = count_calls_to_evolve_one_step + 1 end do - ! Save the next timestep MESA would have required if it was not imposed by end_time - dt_save = s% dt_next - - ! Set the timestep required to reach end_time - s% dt_next = (end_time - s% star_age)*secyer - - ! Evolve one last step to reach end_time - call do_evolve_one_step(id, result, ierr) - if (result /= keep_going) then - if (s% result_reason /= result_reason_normal) then - ierr = -1 - end if - end if - - ! If several evolve steps were performed, set the timestep for the future evolution - ! to the value it would have reached without the condition to finish at end_time - if (count /= 0) then + ! In case evolve_until is called several times subsequently. + ! If there was more than one call to evolve_one_step, sets the next timestep + ! to the value it was before it was set for star_age to reach end_time. + if (count_calls_to_evolve_one_step > 1) then s% dt_next = dt_save - endif + end if s% doing_first_model_of_run = .false.