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..8344e5a0f6 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,43 +662,71 @@ 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_calls_to_evolve_one_step 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 - - s% max_age = s% star_age + delta_t - s% num_adjusted_dt_steps_before_max_age = -1 - - evolve_loop: do + ! (LS 2024) + ! Update of evolve_until in case subroutine is called several times subsequently. + + ! 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. + + ! End time of evolution + end_time = s% star_age + delta_t + + ! 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 < 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) 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 + count_calls_to_evolve_one_step = count_calls_to_evolve_one_step + 1 + end do - ! In case you want to evolve further - s% max_age = old_age - if(s% dt_next==0) s% dt_next = s% dt + ! 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 + end if s% doing_first_model_of_run = .false. end subroutine evolve_until - - ! *********************************************************************** ! Routines for modifying a star ! *********************************************************************** @@ -743,7 +771,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 +813,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 +836,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 +1067,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 +1091,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 +1115,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 +1136,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 +1157,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 +1168,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 +1178,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 +1189,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 +1204,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 +1220,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 +1232,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 +1247,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 +1259,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 +1284,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 +1306,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 +1338,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() 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()