From 0f4a2db9e7f0903d53a624188b629fcf81bf289f Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Mon, 15 Apr 2024 15:56:20 +0200 Subject: [PATCH 01/10] update phantom --- src/amuse/community/phantom/interface.f90 | 116 +++++++++++++++-- src/amuse/community/phantom/interface.py | 152 +++++++++++++++++++++- 2 files changed, 256 insertions(+), 12 deletions(-) diff --git a/src/amuse/community/phantom/interface.f90 b/src/amuse/community/phantom/interface.f90 index 0c79f4a67d..9aee823f6e 100644 --- a/src/amuse/community/phantom/interface.f90 +++ b/src/amuse/community/phantom/interface.f90 @@ -1,6 +1,10 @@ -!MODULE PhantomInterface -! -!CONTAINS +!module amuseinterface +! use allocutils, only:allocate_array +! implicit none +! integer, allocatable:: sph_particle_map(:) +! +! call allocate_array('sph_particle_map', sph_particle_map, maxp) +!contains function initialize_code() use StoppingConditions @@ -8,7 +12,7 @@ function initialize_code() integer :: initialize_code integer :: error double precision :: polyk - call amuse_initialize_code() + call amuse_initialize_code_new() !error = set_support_for_condition(TIMEOUT_DETECTION) !error = set_support_for_condition(NUMBER_OF_STEPS_DETECTION) @@ -202,6 +206,15 @@ function get_internal_energy(index_of_the_particle, u) get_internal_energy=0 end function +function get_particle_type(index_of_the_particle, itype) + implicit none + integer :: index_of_the_particle + integer :: itype + integer :: get_particle_type + !call amuse_get_particle_type(index_of_the_particle, itype) + get_particle_type = -1 +end function + function get_index_of_first_particle(index_of_the_particle) implicit none integer :: index_of_the_particle @@ -241,6 +254,7 @@ function evolve_model(tmax) integer :: is_density_limit_detection_enabled, stopping_index integer :: error double precision :: minimum_density_parameter, maximum_density_parameter, rho, radius + double precision :: time error = reset_stopping_conditions() error = is_stopping_condition_enabled(& @@ -248,11 +262,12 @@ function evolve_model(tmax) error = get_stopping_condition_minimum_density_parameter(minimum_density_parameter) error = get_stopping_condition_maximum_density_parameter(maximum_density_parameter) + call amuse_get_time(time) call amuse_evolve_model(tmax) if (is_density_limit_detection_enabled > 0) then call amuse_get_number_of_sph_particles(nmax) - do i=1, nmax - call amuse_get_radius(i, radius) + do i=1, nmax ! need to loop over sph particles, excluding 'dead' particles! + call amuse_get_radius(i, radius) ! if hsmooth <= 0, particle is dead if (radius > 0) then call amuse_get_density(i, rho) if (& @@ -265,6 +280,8 @@ function evolve_model(tmax) error = set_stopping_condition_particle_index(stopping_index, 0, i) endif endif + ! else + ! nmax += 1 endif enddo endif @@ -475,6 +492,24 @@ function get_radius(index_of_the_particle, radius) get_radius=0 end function +function get_sink_temperature(index_of_the_particle, temperature) + implicit none + integer :: index_of_the_particle + double precision :: temperature + integer :: get_sink_temperature + call amuse_get_sink_temperature(index_of_the_particle, temperature) + get_sink_temperature=0 +end function + +function get_sink_luminosity(index_of_the_particle, luminosity) + implicit none + integer :: index_of_the_particle + double precision :: luminosity + integer :: get_sink_luminosity + call amuse_get_sink_luminosity(index_of_the_particle, luminosity) + get_sink_luminosity=0 +end function + function get_state_star(index_of_the_particle, mass, x, y, z, vx, vy, vz, & tform, radius) implicit none @@ -491,6 +526,14 @@ function set_begin_time(time) set_begin_time=-1 end function +integer function set_phantom_option(name, value, match) + implicit none + character(*), intent(inout):: name, value + logical, intent(inout):: match + call amuse_set_phantom_option(name, value, match) + set_phantom_option = 0 +end function + function new_dm_particle(index_of_the_particle, mass, x, y, z, vx, vy, vz) implicit none integer :: index_of_the_particle @@ -507,6 +550,7 @@ function new_sink_particle(index_of_the_particle, mass, x, y, z, vx, vy, vz, & integer :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz, radius, h_smooth integer :: new_sink_particle + write(*,*) 'particle ', index_of_the_particle, ': radius set to ', radius call amuse_new_sink_particle(index_of_the_particle, mass, x, y, z, & vx, vy, vz, radius, h_smooth) new_sink_particle=0 @@ -521,6 +565,24 @@ function set_radius(index_of_the_particle, radius) set_radius=0 end function +function set_sink_temperature(index_of_the_particle, temperature) + implicit none + integer :: index_of_the_particle + double precision :: temperature + integer :: set_sink_temperature + call amuse_set_sink_temperature(index_of_the_particle, temperature) + set_sink_temperature=0 +end function + +function set_sink_luminosity(index_of_the_particle, luminosity) + implicit none + integer :: index_of_the_particle + double precision :: luminosity + integer :: set_sink_luminosity + call amuse_set_sink_luminosity(index_of_the_particle, luminosity) + set_sink_luminosity=0 +end function + function recommit_parameters() implicit none integer :: recommit_parameters @@ -538,11 +600,23 @@ function get_potential_energy(potential_energy) function get_state_sph(index_of_the_particle, mass, x, y, z, vx, vy, vz, u, & h_smooth) implicit none - integer :: index_of_the_particle + integer :: index_of_the_particle, n double precision :: mass, x, y, z, vx, vy, vz, u, h_smooth integer :: get_state_sph - call amuse_get_state_gas(index_of_the_particle, mass, x, y, z, vx, vy, vz, u, h_smooth) - get_state_sph=0 + get_state_sph = -1 + call amuse_get_number_of_sph_particles(n) + if (index_of_the_particle < 1) then + get_state_sph = -2 + elseif (index_of_the_particle > n) then + !write(*,*) index_of_the_particle, n, "error?" + call amuse_get_state_gas(index_of_the_particle, mass, x, y, z, vx, vy, vz, u, h_smooth) + !write(*,*) mass, x, y, z + !get_state_sph = -3 + get_state_sph = 0 + else + call amuse_get_state_gas(index_of_the_particle, mass, x, y, z, vx, vy, vz, u, h_smooth) + get_state_sph=0 + endif end function function get_gravity_at_point(eps, x, y, z, ax, ay, az, npoints) @@ -601,6 +675,7 @@ function set_position(index_of_the_particle, x, y, z) function commit_parameters() implicit none integer :: commit_parameters + call amuse_commit_parameters() commit_parameters=0 end function @@ -1133,4 +1208,25 @@ function get_constant_planckh(planckh) get_constant_planckh=0 end function -!END MODULE +!integer function get_sonic_type(sonic_type) +! implicit none +! integer :: sonic_type +! call amuse_get_sonic_type(sonic_type) +! return 0 +!end function + +!integer function create_wind(time_step) +! use inject, only: inject_particles +! ! use inject_wind, only: inject_particles +! implicit none +! double precision :: time_step +! +! ! Normally, Phantom would use/fill the position/velocity arrays of +! ! gas/sink particles for this subroutine. +! ! We can just create a new array instead, and return the values to AMUSE +! !call inject_particles(time, dtlast, xyzh, vxyzu, xyzmh_ptmass, vxyz_ptmass,& +! ! npart, npart_old, npartoftype, dtinject) +! create_wind = 0 +!end function + +!end module diff --git a/src/amuse/community/phantom/interface.py b/src/amuse/community/phantom/interface.py index 309f84f626..00ab99874d 100644 --- a/src/amuse/community/phantom/interface.py +++ b/src/amuse/community/phantom/interface.py @@ -42,6 +42,25 @@ def __init__(self, **options): **options) LiteratureReferencesMixIn.__init__(self) + @legacy_function + def set_phantom_option(): + function = LegacyFunctionSpecification() + function.addParameter( + 'name', dtype='string', direction=function.IN + ) + function.addParameter( + 'value', dtype='string', direction=function.IN + ) + function.addParameter( + 'match', dtype='bool', direction=function.OUT + ) + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + -1 - ERROR + """ + return function + @legacy_function def new_dm_particle(): function = LegacyFunctionSpecification() @@ -232,7 +251,11 @@ def get_state_sph(): 0 - OK particle was removed from the model -1 - ERROR - particle could not be found + something went wrong + -2 - ERROR + wrong particle type + -3 - ERROR + particle index too large """ return function @@ -398,6 +421,44 @@ def set_state_sink(): """ return function + @legacy_function + def set_sink_temperature(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN, + description='', + ) + function.addParameter( + 'temperature', dtype='float64', direction=function.IN, + description='', + ) + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + -1 - ERROR + """ + return function + + @legacy_function + def set_sink_luminosity(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN, + description='', + ) + function.addParameter( + 'luminosity', dtype='float64', direction=function.IN, + description='', + ) + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + -1 - ERROR + """ + return function + @legacy_function def set_internal_energy(): function = LegacyFunctionSpecification() @@ -664,6 +725,44 @@ def get_co_abundance(): """ return function + @legacy_function + def get_sink_temperature(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN, + description='' + ) + function.addParameter( + 'temperature', dtype='float64', direction=function.OUT, + description="The effective temperature of the sink particle", + ) + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + -1 - ERROR + """ + return function + + @legacy_function + def get_sink_luminosity(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN, + description='' + ) + function.addParameter( + 'luminosity', dtype='float64', direction=function.OUT, + description="The luminosity of the sink particle", + ) + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + -1 - ERROR + """ + return function + @legacy_function def get_internal_energy(): function = LegacyFunctionSpecification() @@ -1541,7 +1640,8 @@ def __init__( # phantom_time = 60 * 60 * 24 * 365.25 * 1e6 | units.s # phantom_length = (phantom_time**2 * phantom_gg * phantom_mass)**(1/3) phantom_mass = 1.0 * phantom_solarm - phantom_length = 0.1 * phantom_pc + # phantom_length = 0.1 * phantom_pc + phantom_length = 1.0 | units.au phantom_time = (phantom_length**3 / (phantom_gg*phantom_mass))**0.5 unit_converter = ConvertBetweenGenericAndSiUnits( @@ -1926,6 +2026,10 @@ def define_particle_sets(self, handler): handler.add_getter('sink_particles', 'get_acceleration') handler.add_getter('sink_particles', 'get_radius') handler.add_setter('sink_particles', 'set_radius') + handler.add_getter('sink_particles', 'get_sink_temperature') + handler.add_setter('sink_particles', 'set_sink_temperature') + handler.add_getter('sink_particles', 'get_sink_luminosity') + handler.add_setter('sink_particles', 'set_sink_luminosity') handler.add_getter('sink_particles', 'get_smoothing_length') handler.add_setter('sink_particles', 'set_smoothing_length') @@ -2100,6 +2204,50 @@ def define_methods(self, handler): ) ) + handler.add_method( + "set_sink_temperature", + ( + handler.INDEX, + units.K, + ), + ( + handler.ERROR_CODE, + ) + ) + + handler.add_method( + "get_sink_temperature", + ( + handler.INDEX, + ), + ( + units.K, + handler.ERROR_CODE, + ) + ) + + handler.add_method( + "set_sink_luminosity", + ( + handler.INDEX, + units.LSun, + ), + ( + handler.ERROR_CODE, + ) + ) + + handler.add_method( + "get_sink_luminosity", + ( + handler.INDEX, + ), + ( + units.LSun, + handler.ERROR_CODE, + ) + ) + handler.add_method( "get_density", ( From 73f690c2f8d4f6fc1900d43dc4b3dd738d5a0f92 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Mon, 15 Apr 2024 15:57:47 +0200 Subject: [PATCH 02/10] update phantom 2 --- src/amuse/community/phantom/download.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/amuse/community/phantom/download.py b/src/amuse/community/phantom/download.py index 959efcc131..5f5a3d63f7 100755 --- a/src/amuse/community/phantom/download.py +++ b/src/amuse/community/phantom/download.py @@ -70,7 +70,8 @@ def new_option_parser(): result = OptionParser() result.add_option( "--version", - default="ca6907f5f9a95f866b1d7e520cc356ab8cec8dd0", + # default="ca6907f5f9a95f866b1d7e520cc356ab8cec8dd0", + default="1e7bd070b79f2408e333797f544c3b1daa047c8a", dest="version", help="version number to download", type="string" From cd7926e6a4865c56c479ac4e74dfe9d64123b8cd Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Tue, 4 Jun 2024 14:01:08 +0200 Subject: [PATCH 03/10] distinction between effective radius and accretion radius --- src/amuse/community/phantom/interface.f90 | 37 +++++++++--- src/amuse/community/phantom/interface.py | 74 +++++++++++++++++++++++ 2 files changed, 101 insertions(+), 10 deletions(-) diff --git a/src/amuse/community/phantom/interface.f90 b/src/amuse/community/phantom/interface.f90 index 9aee823f6e..a339964bd1 100644 --- a/src/amuse/community/phantom/interface.f90 +++ b/src/amuse/community/phantom/interface.f90 @@ -371,24 +371,24 @@ function synchronize_model() end function function set_state_sink(index_of_the_particle, mass, x, y, z, & - vx, vy, vz, radius, h_smooth) + vx, vy, vz, radius, accretion_radius, h_smooth) implicit none integer :: index_of_the_particle - double precision :: mass, x, y, z, vx, vy, vz, radius, h_smooth + double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius, h_smooth integer :: set_state_sink call amuse_set_state_sink(index_of_the_particle, mass, x, y, z, & - vx, vy, vz, radius, h_smooth) + vx, vy, vz, radius, accretion_radius, h_smooth) set_state_sink=0 end function function get_state_sink(index_of_the_particle, mass, x, y, z, & - vx, vy, vz, radius, h_smooth) + vx, vy, vz, radius, accretion_radius, h_smooth) implicit none integer :: index_of_the_particle - double precision :: mass, x, y, z, vx, vy, vz, radius, h_smooth + double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius, h_smooth integer :: get_state_sink call amuse_get_state_sink(index_of_the_particle, mass, x, y, z, & - vx, vy, vz, radius, h_smooth) + vx, vy, vz, radius, accretion_radius, h_smooth) get_state_sink=0 end function @@ -492,6 +492,15 @@ function get_radius(index_of_the_particle, radius) get_radius=0 end function +function get_sink_accretion_radius(index_of_the_particle, radius) + implicit none + integer :: index_of_the_particle + double precision :: radius + integer :: get_sink_accretion_radius + call amuse_get_sink_accretion_radius(index_of_the_particle, radius) + get_sink_accretion_radius=0 +end function + function get_sink_temperature(index_of_the_particle, temperature) implicit none integer :: index_of_the_particle @@ -545,14 +554,13 @@ function new_dm_particle(index_of_the_particle, mass, x, y, z, vx, vy, vz) end function function new_sink_particle(index_of_the_particle, mass, x, y, z, vx, vy, vz, & - radius, h_smooth) + radius, accretion_radius, h_smooth) implicit none integer :: index_of_the_particle - double precision :: mass, x, y, z, vx, vy, vz, radius, h_smooth + double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius, h_smooth integer :: new_sink_particle - write(*,*) 'particle ', index_of_the_particle, ': radius set to ', radius call amuse_new_sink_particle(index_of_the_particle, mass, x, y, z, & - vx, vy, vz, radius, h_smooth) + vx, vy, vz, radius, accretion_radius, h_smooth) new_sink_particle=0 end function @@ -565,6 +573,15 @@ function set_radius(index_of_the_particle, radius) set_radius=0 end function +function set_sink_accretion_radius(index_of_the_particle, radius) + implicit none + integer :: index_of_the_particle + double precision :: radius + integer :: set_sink_accretion_radius + call amuse_set_sink_accretion_radius(index_of_the_particle, radius) + set_sink_accretion_radius=0 +end function + function set_sink_temperature(index_of_the_particle, temperature) implicit none integer :: index_of_the_particle diff --git a/src/amuse/community/phantom/interface.py b/src/amuse/community/phantom/interface.py index 00ab99874d..ec1776cc79 100644 --- a/src/amuse/community/phantom/interface.py +++ b/src/amuse/community/phantom/interface.py @@ -102,6 +102,9 @@ def new_sink_particle(): function.addParameter(x, dtype='float64', direction=function.IN) function.addParameter( 'radius', dtype='float64', direction=function.IN, default=0.01, + ) + function.addParameter( + 'accretion_radius', dtype='float64', direction=function.IN, default=0.01, # default should be h_acc ) function.addParameter( @@ -190,6 +193,9 @@ def get_state_sink(): description="The current velocity vector of the particle") function.addParameter( 'radius', dtype='float64', direction=function.OUT, + description="The effective radius of the particle") + function.addParameter( + 'accretion_radius', dtype='float64', direction=function.OUT, description="The accretion radius of the particle") function.addParameter( 'h_smooth', dtype='float64', direction=function.OUT, @@ -404,6 +410,9 @@ def set_state_sink(): description="The new velocity vector of the particle") function.addParameter( 'radius', dtype='float64', direction=function.IN, + description="The effective radius of the particle") + function.addParameter( + 'accretion_radius', dtype='float64', direction=function.IN, description="The accretion radius of the particle") function.addParameter( 'h_smooth', dtype='float64', direction=function.IN, @@ -421,6 +430,25 @@ def set_state_sink(): """ return function + @legacy_function + def set_sink_accretion_radius(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN, + description='', + ) + function.addParameter( + 'accretion_radius', dtype='float64', direction=function.IN, + description='', + ) + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + -1 - ERROR + """ + return function + @legacy_function def set_sink_temperature(): function = LegacyFunctionSpecification() @@ -725,6 +753,25 @@ def get_co_abundance(): """ return function + @legacy_function + def get_sink_accretion_radius(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN, + description='' + ) + function.addParameter( + 'accretion_radius', dtype='float64', direction=function.OUT, + description="The accretion radius of the sink particle", + ) + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + -1 - ERROR + """ + return function + @legacy_function def get_sink_temperature(): function = LegacyFunctionSpecification() @@ -2026,6 +2073,8 @@ def define_particle_sets(self, handler): handler.add_getter('sink_particles', 'get_acceleration') handler.add_getter('sink_particles', 'get_radius') handler.add_setter('sink_particles', 'set_radius') + handler.add_getter('sink_particles', 'get_sink_accretion_radius') + handler.add_setter('sink_particles', 'set_sink_accretion_radius') handler.add_getter('sink_particles', 'get_sink_temperature') handler.add_setter('sink_particles', 'set_sink_temperature') handler.add_getter('sink_particles', 'get_sink_luminosity') @@ -2086,6 +2135,7 @@ def define_methods(self, handler): generic_unit_system.speed, generic_unit_system.length, generic_unit_system.length, + generic_unit_system.length, ), ( handler.INDEX, @@ -2181,6 +2231,7 @@ def define_methods(self, handler): generic_unit_system.speed, generic_unit_system.length, generic_unit_system.length, + generic_unit_system.length, handler.ERROR_CODE, ) ) @@ -2198,12 +2249,35 @@ def define_methods(self, handler): generic_unit_system.speed, generic_unit_system.length, generic_unit_system.length, + generic_unit_system.length, + ), + ( + handler.ERROR_CODE, + ) + ) + + handler.add_method( + "set_sink_accretion_radius", + ( + handler.INDEX, + generic_unit_system.length, ), ( handler.ERROR_CODE, ) ) + handler.add_method( + "get_sink_accretion_radius", + ( + handler.INDEX, + ), + ( + generic_unit_system.length, + handler.ERROR_CODE, + ) + ) + handler.add_method( "set_sink_temperature", ( From d080934aade3a18c87ed01b665ed455e1390a28b Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Mon, 12 Aug 2024 17:20:10 +0200 Subject: [PATCH 04/10] unit fixes --- src/amuse/community/phantom/interface.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/amuse/community/phantom/interface.py b/src/amuse/community/phantom/interface.py index ec1776cc79..89d481cfa0 100644 --- a/src/amuse/community/phantom/interface.py +++ b/src/amuse/community/phantom/interface.py @@ -1680,6 +1680,7 @@ def __init__( # Not doing this *really* complicates things as we'd need to change the # internal units used in Phantom as well. phantom_solarm = 1.9891e33 | units.g + phantom_au = 1.496e13 | units.cm phantom_pc = 3.086e18 | units.cm phantom_gg = 6.672041e-8 | units.cm**3 * units.g**-1 * units.s**-2 # phantom_mass = 1.0 | units.MSun @@ -1688,7 +1689,7 @@ def __init__( # phantom_length = (phantom_time**2 * phantom_gg * phantom_mass)**(1/3) phantom_mass = 1.0 * phantom_solarm # phantom_length = 0.1 * phantom_pc - phantom_length = 1.0 | units.au + phantom_length = 1.0 * phantom_au phantom_time = (phantom_length**3 / (phantom_gg*phantom_mass))**0.5 unit_converter = ConvertBetweenGenericAndSiUnits( @@ -1731,6 +1732,8 @@ def define_state(self, handler): handler.add_transition('END', 'INITIALIZED', 'initialize_code', False) handler.add_method('END', 'initialize_code') + handler.add_method('INITIALIZED', 'set_phantom_option') + handler.add_method('EDIT', 'set_phantom_option') handler.add_transition('RUN', 'UPDATE', 'new_sph_particle', False) handler.add_method('EDIT', 'new_sph_particle') handler.add_method('UPDATE', 'new_sph_particle') @@ -2304,7 +2307,7 @@ def define_methods(self, handler): "set_sink_luminosity", ( handler.INDEX, - units.LSun, + generic_unit_system.energy / generic_unit_system.time, ), ( handler.ERROR_CODE, @@ -2317,7 +2320,7 @@ def define_methods(self, handler): handler.INDEX, ), ( - units.LSun, + generic_unit_system.energy / generic_unit_system.time, handler.ERROR_CODE, ) ) From f8eed19d6458220094ba684b223b73ff6fa96761 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Mon, 14 Oct 2024 13:27:44 +0200 Subject: [PATCH 05/10] updates to phantom --- src/amuse/community/phantom/interface.py | 59 ++++++++++++++++++++++-- 1 file changed, 54 insertions(+), 5 deletions(-) diff --git a/src/amuse/community/phantom/interface.py b/src/amuse/community/phantom/interface.py index 89d481cfa0..a4f9486f90 100644 --- a/src/amuse/community/phantom/interface.py +++ b/src/amuse/community/phantom/interface.py @@ -5,6 +5,7 @@ LegacyFunctionSpecification, legacy_function, LiteratureReferencesMixIn, + CodeWithDataDirectories, ) from amuse.community.interface.gd import ( @@ -22,11 +23,12 @@ class PhantomInterface( - CodeInterface, - LiteratureReferencesMixIn, - GravitationalDynamicsInterface, - StoppingConditionInterface, - # SinglePointGravityFieldInterface, + CodeInterface, + GravitationalDynamicsInterface, + LiteratureReferencesMixIn, + StoppingConditionInterface, + # GravityFieldInterface, + CodeWithDataDirectories, ): """ The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. @@ -34,6 +36,7 @@ class PhantomInterface( References: .. [#] ADS:2018PASA...35...31P (Price et al., 2018, PASA, Volume 35, id.e031 82 pp) """ + use_modules = ["StoppingConditions", "AmuseInterface"] def __init__(self, **options): CodeInterface.__init__( @@ -73,6 +76,14 @@ def new_dm_particle(): function.result_type = 'int32' return function + @legacy_function + def get_maximum_particle_index(): + function = LegacyFunctionSpecification() + function.addParameter( + 'norig', dtype='int32', direction=function.OUT) + function.result_type = 'int32' + return function + def new_particle(self, mass, x, y, z, vx, vy, vz): return self.new_dm_particle(mass, x, y, z, vx, vy, vz) @@ -1667,6 +1678,28 @@ def set_unit_time(): """ return function + @legacy_function + def get_number_of_new_gas_particles(): + """ + Returns the number of gas particles added since the last call to + update_particles. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "number_of_particles", dtype="int32", direction=function.OUT + ) + function.result_type = "int32" + return function + + @legacy_function + def reset_number_of_new_gas_particles(): + """ + Resets the counter for new gas particles to zero. + """ + function = LegacyFunctionSpecification() + function.result_type = "int32" + return function + class Phantom(GravitationalDynamics, GravityFieldCode): __interface__ = PhantomInterface @@ -1724,6 +1757,22 @@ def initialize_code(self): # self.set_unit_time(self.unit_converter.to_si(generic_unit_system.time).value_in(units.s)) return result + def update_gas_particle_set(self): + """ + Updates the gas particle set after changes in the code (e.g. added wind + particles). + """ + number_of_new_gas_particles = self.get_number_of_new_gas_particles() + print(f"new particles: {number_of_new_gas_particles}") + if number_of_new_gas_particles: + nmax = self.get_maximum_particle_index() + self.gas_particles._private.attribute_storage._add_indices( + list(range(1 + nmax - number_of_new_gas_particles, 1 + nmax)) + ) + self.reset_number_of_new_gas_particles() + deleted_particles = self.gas_particles[self.gas_particles.mass.number == 0] + self.gas_particles.remove_particles(deleted_particles) + def define_state(self, handler): GravitationalDynamics.define_state(self, handler) GravityFieldCode.define_state(self, handler) From 9dc4319f6a92a8245f7dd23f0dc8598aa22dda74 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Mon, 14 Oct 2024 13:28:46 +0200 Subject: [PATCH 06/10] updates to Phantom --- src/amuse/community/phantom/Makefile | 27 +++++---- src/amuse/community/phantom/interface.f90 | 70 +++++++++++++++++------ 2 files changed, 70 insertions(+), 27 deletions(-) diff --git a/src/amuse/community/phantom/Makefile b/src/amuse/community/phantom/Makefile index dbcd5da12c..d115a1c67c 100644 --- a/src/amuse/community/phantom/Makefile +++ b/src/amuse/community/phantom/Makefile @@ -15,9 +15,10 @@ OBJS = interface.o CODELIBDIR = src/phantom/bin CODELIB = src/phantom/bin/libphantom-amuse.a -FCFLAGS += -L${CODELIBDIR} -lphantom-amuse +#FCFLAGS += -L$(CODELIBDIR) -lphantom-amuse CODEDIR = src/phantom +BUILDDIR = $(CODEDIR)/build SETUP = amuse SYSTEM = gfortran @@ -27,28 +28,32 @@ DOWNLOAD_FROM_WEB = $(PYTHON) ./download.py all: phantom_worker -$(CODEDIR)/Makefile: - make -C . download +$(CODEDIR)/Makefile: download + #make -C . download download: - $(RM) -Rf src - mkdir src - $(DOWNLOAD_FROM_WEB) + echo "not downloading" + #$(RM) -Rf src + #mkdir src + #$(DOWNLOAD_FROM_WEB) clean: $(RM) -f *.so *.o *.pyc worker_code.cc worker_code.h $(RM) *~ worker_code worker_code.f90 - $(RM) ${CODELIB} - make -C ${CODEDIR} clean + $(RM) $(CODELIB) + make -C $(CODEDIR) clean -$(CODELIB):$(CODEDIR)/Makefile download - make -C ${CODEDIR} libphantom-amuse SETUP=$(SETUP) SYSTEM=$(SYSTEM) OPENMP=yes +$(CODELIB): $(CODEDIR)/Makefile + make -C $(CODEDIR) libphantom-amuse SETUP=$(SETUP) SYSTEM=$(SYSTEM) OPENMP=yes FC="$(FC)" worker_code.f90: interface.py $(CODE_GENERATOR) --type=f90 interface.py PhantomInterface -o $@ +interface.o: interface.f90 + $(MPIFC) $(FCFLAGS) $(SC_FLAGS) -I$(BUILDDIR) -c -o $@ $< + phantom_worker: worker_code.f90 $(CODELIB) $(OBJS) - $(MPIFC) $(FCFLAGS) $(FS_FLAGS) $< $(OBJS) $(CODELIB) $(SC_FCLIBS) $(FS_LIBS) -fopenmp -o $@ + $(MPIFC) $(FCFLAGS) $(SC_FLAGS) $(FS_FLAGS) $(LDFLAGS) -I$(BUILDDIR) $< interface.o -o $@ $(CODELIB) $(SC_FCLIBS) $(FS_LIBS) -fopenmp %.o: %.f90 $(FC) $(SC_FLAGS) $(FCFLAGS) -c -o $@ $< diff --git a/src/amuse/community/phantom/interface.f90 b/src/amuse/community/phantom/interface.f90 index 620ad188df..dc94961ad5 100644 --- a/src/amuse/community/phantom/interface.f90 +++ b/src/amuse/community/phantom/interface.f90 @@ -1,3 +1,8 @@ +module AmuseInterface + use AmusePhantom + use StoppingConditions + +contains !module amuseinterface ! use allocutils, only:allocate_array ! implicit none @@ -6,20 +11,19 @@ ! call allocate_array('sph_particle_map', sph_particle_map, maxp) !contains -function initialize_code() - use StoppingConditions +integer function initialize_code() result(ret) implicit none - integer :: initialize_code + !integer :: initialize_code integer :: error double precision :: polyk - call amuse_initialize_code_new() + call amuse_initialize_code() !error = set_support_for_condition(TIMEOUT_DETECTION) !error = set_support_for_condition(NUMBER_OF_STEPS_DETECTION) !error = set_support_for_condition(OUT_OF_BOX_DETECTION) error = set_support_for_condition(DENSITY_LIMIT_DETECTION) !error = set_support_for_condition(INTERNAL_ENERGY_LIMIT_DETECTION) - initialize_code=0 + ret=0 end function function cleanup_code() @@ -29,6 +33,14 @@ function cleanup_code() cleanup_code=0 end function +function reset_new_particles_counter() + use AmusePhantom, only: new_particles_since_last_update + implicit none + integer :: reset_new_particles_counter + new_particles_since_last_update = 0 + reset_new_particles_counter = 0 +end function + function commit_particles() implicit none integer :: commit_particles @@ -337,6 +349,23 @@ function get_index_of_next_particle(index_of_the_particle, & get_index_of_next_particle=-1 end function +function get_number_of_new_gas_particles(number) + use AmusePhantom, only: new_particles_since_last_update + implicit none + integer :: get_number_of_new_gas_particles + integer :: number + number = new_particles_since_last_update + get_number_of_new_gas_particles = 0 +end function + +function reset_number_of_new_gas_particles() + use AmusePhantom, only: new_particles_since_last_update + implicit none + integer :: reset_number_of_new_gas_particles + new_particles_since_last_update = 0 + reset_number_of_new_gas_particles = 0 +end function + function new_sph_particle(index_of_the_particle, mass, x, y, z, & vx, vy, vz, u, h_smooth) implicit none @@ -377,7 +406,7 @@ function set_state_sink(index_of_the_particle, mass, x, y, z, & double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius, h_smooth integer :: set_state_sink call amuse_set_state_sink(index_of_the_particle, mass, x, y, z, & - vx, vy, vz, radius, accretion_radius, h_smooth) + vx, vy, vz, radius, accretion_radius) set_state_sink=0 end function @@ -388,7 +417,7 @@ function get_state_sink(index_of_the_particle, mass, x, y, z, & double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius, h_smooth integer :: get_state_sink call amuse_get_state_sink(index_of_the_particle, mass, x, y, z, & - vx, vy, vz, radius, accretion_radius, h_smooth) + vx, vy, vz, radius, accretion_radius) get_state_sink=0 end function @@ -453,6 +482,14 @@ function get_thermal_energy(thermal_energy) get_thermal_energy=0 end function +function get_maximum_particle_index(i) + implicit none + integer :: i + integer :: get_maximum_particle_index + call amuse_get_norig(i) + get_maximum_particle_index=0 +end function + function get_number_of_particles(n) implicit none integer :: n @@ -538,10 +575,10 @@ integer function set_phantom_option(name, value, match) function new_dm_particle(index_of_the_particle, mass, x, y, z, vx, vy, vz) implicit none integer :: index_of_the_particle - double precision :: mass, x, y, z, vx, vy, vz, radius + double precision :: mass, x, y, z, vx, vy, vz integer :: new_dm_particle call amuse_new_dm_particle(index_of_the_particle, mass, x, y, z, & - vx, vy, vz, radius) + vx, vy, vz) new_dm_particle=0 end function @@ -606,25 +643,26 @@ function get_potential_energy(potential_energy) get_potential_energy=0 end function -function get_state_sph(index_of_the_particle, mass, x, y, z, vx, vy, vz, u, & +function get_state_sph(index_of_the_particle, mass, x, y, z, vx, vy, vz, u, & h_smooth) implicit none integer :: index_of_the_particle, n double precision :: mass, x, y, z, vx, vy, vz, u, h_smooth integer :: get_state_sph get_state_sph = -1 - call amuse_get_number_of_sph_particles(n) + !call amuse_get_number_of_sph_particles(n) + call amuse_get_norig(n) if (index_of_the_particle < 1) then get_state_sph = -2 elseif (index_of_the_particle > n) then !write(*,*) index_of_the_particle, n, "error?" - call amuse_get_state_gas(index_of_the_particle, mass, x, y, z, vx, vy, vz, u, h_smooth) + !call amuse_get_state_gas(index_of_the_particle, mass, x, y, z, vx, vy, vz, u, h_smooth, phantom_index) !write(*,*) mass, x, y, z - !get_state_sph = -3 - get_state_sph = 0 + get_state_sph = -3 + !get_state_sph = 0 else call amuse_get_state_gas(index_of_the_particle, mass, x, y, z, vx, vy, vz, u, h_smooth) - get_state_sph=0 + get_state_sph = 0 endif end function @@ -1238,4 +1276,4 @@ function get_constant_planckh(planckh) ! create_wind = 0 !end function -!end module +end module From eaa2c5e03a53e785175f181454d7dd40cd655bda Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Fri, 18 Oct 2024 12:14:37 +0200 Subject: [PATCH 07/10] Add 64 bits indices to interface --- .gitignore | 2 - src/amuse/community/interface/gd/__init__.py | 23 + .../{gd.py => gd/gravitational_dynamics.py} | 1037 +++++------ .../interface/gd/gravitational_dynamics_64.py | 1512 +++++++++++++++++ .../community/interface/gd/gravity_field.py | 113 ++ src/amuse/community/phantom/interface.f90 | 109 +- src/amuse/community/phantom/interface.py | 78 +- 7 files changed, 2274 insertions(+), 600 deletions(-) create mode 100644 src/amuse/community/interface/gd/__init__.py rename src/amuse/community/interface/{gd.py => gd/gravitational_dynamics.py} (62%) create mode 100644 src/amuse/community/interface/gd/gravitational_dynamics_64.py create mode 100644 src/amuse/community/interface/gd/gravity_field.py diff --git a/.gitignore b/.gitignore index 86281d0f39..6009623c25 100644 --- a/.gitignore +++ b/.gitignore @@ -405,8 +405,6 @@ test_results/code-sockets.f90 test_results/code.c test_results/code.cc test_results/code.f90 -gas_interface -interface test_results/interface-sockets.f90 test_results/interface.cc test_results/interface.f90 diff --git a/src/amuse/community/interface/gd/__init__.py b/src/amuse/community/interface/gd/__init__.py new file mode 100644 index 0000000000..98cdbabcec --- /dev/null +++ b/src/amuse/community/interface/gd/__init__.py @@ -0,0 +1,23 @@ +from amuse.support.interface import InCodeComponentImplementation +from amuse.units import nbody_system +from amuse.units import generic_unit_converter +from amuse.community.interface import common + +from amuse.rfi.core import legacy_function +from amuse.rfi.core import LegacyFunctionSpecification + +from .gravitational_dynamics import ( + GravitationalDynamicsInterface, + GravitationalDynamicsDocumentation, + GravitationalDynamics, +) +from .gravity_field import ( + SinglePointGravityFieldInterface, + GravityFieldInterface, + GravityFieldCode, +) +from .gravitational_dynamics_64 import ( + GravitationalDynamics64Interface, + GravitationalDynamics64Documentation, + GravitationalDynamics64, +) diff --git a/src/amuse/community/interface/gd.py b/src/amuse/community/interface/gd/gravitational_dynamics.py similarity index 62% rename from src/amuse/community/interface/gd.py rename to src/amuse/community/interface/gd/gravitational_dynamics.py index 8de5a44b2d..c3bc929379 100644 --- a/src/amuse/community/interface/gd.py +++ b/src/amuse/community/interface/gd/gravitational_dynamics.py @@ -1,10 +1,8 @@ """ -Stellar Dynamics Interface Defintion +Stellar Dynamics Interface Definition """ -from amuse.support.interface import InCodeComponentImplementation from amuse.units import nbody_system -from amuse.units import generic_unit_converter from amuse.community.interface import common from amuse.rfi.core import legacy_function @@ -23,40 +21,67 @@ def new_particle(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.OUT, + "index_of_the_particle", + dtype="int32", + direction=function.OUT, description=( "An index assigned to the newly created particle. " "This index is supposed to be a local index for the code " "(and not valid in other instances of the code or in other " "codes)" - ) + ), ) function.addParameter( - 'mass', dtype='float64', direction=function.IN, - description="The mass of the particle") + "mass", + dtype="float64", + direction=function.IN, + description="The mass of the particle", + ) function.addParameter( - 'x', dtype='float64', direction=function.IN, - description="The initial position vector of the particle") + "x", + dtype="float64", + direction=function.IN, + description="The initial position vector of the particle", + ) function.addParameter( - 'y', dtype='float64', direction=function.IN, - description="The initial position vector of the particle") + "y", + dtype="float64", + direction=function.IN, + description="The initial position vector of the particle", + ) function.addParameter( - 'z', dtype='float64', direction=function.IN, - description="The initial position vector of the particle") + "z", + dtype="float64", + direction=function.IN, + description="The initial position vector of the particle", + ) function.addParameter( - 'vx', dtype='float64', direction=function.IN, - description="The initial velocity vector of the particle") + "vx", + dtype="float64", + direction=function.IN, + description="The initial velocity vector of the particle", + ) function.addParameter( - 'vy', dtype='float64', direction=function.IN, - description="The initial velocity vector of the particle") + "vy", + dtype="float64", + direction=function.IN, + description="The initial velocity vector of the particle", + ) function.addParameter( - 'vz', dtype='float64', direction=function.IN, - description="The initial velocity vector of the particle") + "vz", + dtype="float64", + direction=function.IN, + description="The initial velocity vector of the particle", + ) function.addParameter( - 'radius', dtype='float64', direction=function.IN, - description="The radius of the particle", default=0) - function.result_type = 'int32' + "radius", + dtype="float64", + direction=function.IN, + description="The radius of the particle", + default=0, + ) + function.result_type = "int32" function.result_doc = """ 0 - OK particle was created and added to the model -1 - ERROR @@ -74,13 +99,15 @@ def delete_particle(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle to be removed. This index must have " "been returned by an earlier call to :meth:`new_particle`" - ) + ), ) - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK particle was removed from the model @@ -101,37 +128,63 @@ def get_state(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle to get the state from. This index must " "have been returned by an earlier call to :meth:`new_particle`" - ) + ), ) function.addParameter( - 'mass', dtype='float64', direction=function.OUT, - description="The current mass of the particle") + "mass", + dtype="float64", + direction=function.OUT, + description="The current mass of the particle", + ) function.addParameter( - 'x', dtype='float64', direction=function.OUT, - description="The current position vector of the particle") + "x", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) function.addParameter( - 'y', dtype='float64', direction=function.OUT, - description="The current position vector of the particle") + "y", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) function.addParameter( - 'z', dtype='float64', direction=function.OUT, - description="The current position vector of the particle") + "z", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) function.addParameter( - 'vx', dtype='float64', direction=function.OUT, - description="The current velocity vector of the particle") + "vx", + dtype="float64", + direction=function.OUT, + description="The current velocity vector of the particle", + ) function.addParameter( - 'vy', dtype='float64', direction=function.OUT, - description="The current velocity vector of the particle") + "vy", + dtype="float64", + direction=function.OUT, + description="The current velocity vector of the particle", + ) function.addParameter( - 'vz', dtype='float64', direction=function.OUT, - description="The current velocity vector of the particle") + "vz", + dtype="float64", + direction=function.OUT, + description="The current velocity vector of the particle", + ) function.addParameter( - 'radius', dtype='float64', direction=function.OUT, - description="The current radius of the particle") - function.result_type = 'int32' + "radius", + dtype="float64", + direction=function.OUT, + description="The current radius of the particle", + ) + function.result_type = "int32" function.result_doc = """ 0 - OK particle was removed from the model @@ -150,38 +203,65 @@ def set_state(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle for which the state is to be updated. " "This index must have been returned by an earlier call to " ":meth:`new_particle`" - ) + ), ) function.addParameter( - 'mass', dtype='float64', direction=function.IN, - description="The new mass of the particle") + "mass", + dtype="float64", + direction=function.IN, + description="The new mass of the particle", + ) function.addParameter( - 'x', dtype='float64', direction=function.IN, - description="The new position vector of the particle") + "x", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) function.addParameter( - 'y', dtype='float64', direction=function.IN, - description="The new position vector of the particle") + "y", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) function.addParameter( - 'z', dtype='float64', direction=function.IN, - description="The new position vector of the particle") + "z", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) function.addParameter( - 'vx', dtype='float64', direction=function.IN, - description="The new velocity vector of the particle") + "vx", + dtype="float64", + direction=function.IN, + description="The new velocity vector of the particle", + ) function.addParameter( - 'vy', dtype='float64', direction=function.IN, - description="The new velocity vector of the particle") + "vy", + dtype="float64", + direction=function.IN, + description="The new velocity vector of the particle", + ) function.addParameter( - 'vz', dtype='float64', direction=function.IN, - description="The new velocity vector of the particle") + "vz", + dtype="float64", + direction=function.IN, + description="The new velocity vector of the particle", + ) function.addParameter( - 'radius', dtype='float64', direction=function.IN, - description="The new radius of the particle", default=0) - function.result_type = 'int32' + "radius", + dtype="float64", + direction=function.IN, + description="The new radius of the particle", + default=0, + ) + function.result_type = "int32" function.result_doc = """ 0 - OK particle was found in the model and the information was set @@ -202,17 +282,21 @@ def get_mass(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle to get the state from. This index must " "have been returned by an earlier call to :meth:`new_particle`" - ) + ), ) function.addParameter( - 'mass', dtype='float64', direction=function.OUT, - description="The current mass of the particle" + "mass", + dtype="float64", + direction=function.OUT, + description="The current mass of the particle", ) - function.result_type = 'int32' + function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK @@ -229,18 +313,22 @@ def set_mass(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle for which the state is to be updated. " "This index must have been returned by an earlier call to " ":meth:`new_particle`" - ) + ), ) function.addParameter( - 'mass', dtype='float64', direction=function.IN, - description="The new mass of the particle" + "mass", + dtype="float64", + direction=function.IN, + description="The new mass of the particle", ) - function.result_type = 'int32' + function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK @@ -260,17 +348,21 @@ def get_radius(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle to get the radius of. This index must " "have been returned by an earlier call to :meth:`new_particle`" - ) + ), ) function.addParameter( - 'radius', dtype='float64', direction=function.OUT, - description="The current radius of the particle" + "radius", + dtype="float64", + direction=function.OUT, + description="The current radius of the particle", ) - function.result_type = 'int32' + function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK @@ -288,17 +380,21 @@ def set_radius(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle to get the radius of. This index must " "have been returned by an earlier call to :meth:`new_particle`" - ) + ), ) function.addParameter( - 'radius', dtype='float64', direction=function.IN, - description="The new radius of the particle" + "radius", + dtype="float64", + direction=function.IN, + description="The new radius of the particle", ) - function.result_type = 'int32' + function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK @@ -316,22 +412,33 @@ def get_position(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle to get the state from. This index must " "have been returned by an earlier call to :meth:`new_particle`" - ) + ), ) function.addParameter( - 'x', dtype='float64', direction=function.OUT, - description="The current position vector of the particle") + "x", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) function.addParameter( - 'y', dtype='float64', direction=function.OUT, - description="The current position vector of the particle") + "y", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) function.addParameter( - 'z', dtype='float64', direction=function.OUT, - description="The current position vector of the particle") - function.result_type = 'int32' + "z", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) + function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK @@ -350,23 +457,34 @@ def set_position(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle for which the state is to be updated. " "This index must have been returned by an earlier call to " ":meth:`new_particle`" - ) + ), ) function.addParameter( - 'x', dtype='float64', direction=function.IN, - description="The new position vector of the particle") + "x", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) function.addParameter( - 'y', dtype='float64', direction=function.IN, - description="The new position vector of the particle") + "y", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) function.addParameter( - 'z', dtype='float64', direction=function.IN, - description="The new position vector of the particle") - function.result_type = 'int32' + "z", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) + function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK @@ -386,34 +504,40 @@ def get_velocity(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle to get the velocity from. This index " "must have been returned by an earlier call to " ":meth:`new_particle`" - ) + ), ) function.addParameter( - 'vx', dtype='float64', direction=function.OUT, + "vx", + dtype="float64", + direction=function.OUT, description=( - "The current x component of the position vector of the " - "particle" - ) + "The current x component of the position vector of the " "particle" + ), ) function.addParameter( - 'vy', dtype='float64', direction=function.OUT, + "vy", + dtype="float64", + direction=function.OUT, description=( - "The current y component of the position vector of the " - "particle" - ) + "The current y component of the position vector of the " "particle" + ), ) function.addParameter( - 'vz', dtype='float64', direction=function.OUT, + "vz", + dtype="float64", + direction=function.OUT, description=( - "The current z component of the position vector of the " - "particle") + "The current z component of the position vector of the " "particle" + ), ) - function.result_type = 'int32' + function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK @@ -432,34 +556,39 @@ def set_velocity(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle to get the state from. This index must " "have been returned by an earlier call to :meth:`new_particle`" - ) + ), ) function.addParameter( - 'vx', dtype='float64', direction=function.IN, + "vx", + dtype="float64", + direction=function.IN, description=( - "The current x component of the velocity vector of the " - "particle" - ) + "The current x component of the velocity vector of the " "particle" + ), ) function.addParameter( - 'vy', dtype='float64', direction=function.IN, + "vy", + dtype="float64", + direction=function.IN, description=( - "The current y component of the velocity vector of the " - "particle" - ) + "The current y component of the velocity vector of the " "particle" + ), ) function.addParameter( - 'vz', dtype='float64', direction=function.IN, + "vz", + dtype="float64", + direction=function.IN, description=( - "The current z component of the velocity vector of the " - "particle" - ) + "The current z component of the velocity vector of the " "particle" + ), ) - function.result_type = 'int32' + function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK @@ -479,22 +608,33 @@ def get_acceleration(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle to get the state from. This index must " "have been returned by an earlier call to :meth:`new_particle`" - ) + ), ) function.addParameter( - 'ax', dtype='float64', direction=function.OUT, - description="The current acceleration vector of the particle") + "ax", + dtype="float64", + direction=function.OUT, + description="The current acceleration vector of the particle", + ) function.addParameter( - 'ay', dtype='float64', direction=function.OUT, - description="The current acceleration vector of the particle") + "ay", + dtype="float64", + direction=function.OUT, + description="The current acceleration vector of the particle", + ) function.addParameter( - 'az', dtype='float64', direction=function.OUT, - description="The current acceleration vector of the particle") - function.result_type = 'int32' + "az", + dtype="float64", + direction=function.OUT, + description="The current acceleration vector of the particle", + ) + function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK @@ -515,17 +655,22 @@ def get_potential(): function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + "index_of_the_particle", + dtype="int32", + direction=function.IN, description=( "Index of the particle for which the state is to be updated. " "This index must have been returned by an earlier call to " ":meth:`new_particle`" - ) + ), ) function.addParameter( - 'potential', dtype='float64', direction=function.OUT, - description="The current scalar potential...") - function.result_type = 'int32' + "potential", + dtype="float64", + direction=function.OUT, + description="The current scalar potential...", + ) + function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK @@ -545,13 +690,15 @@ def evolve_model(): """ function = LegacyFunctionSpecification() function.addParameter( - 'time', dtype='float64', direction=function.IN, + "time", + dtype="float64", + direction=function.IN, description=( "Model time to evolve the code to. The model will be " "evolved until this time is reached exactly or just after." - ) + ), ) - function.result_type = 'int32' + function.result_type = "int32" return function @legacy_function @@ -563,7 +710,7 @@ def commit_particles(): after the last new_particle call. """ function = LegacyFunctionSpecification() - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK Model is initialized and evolution can start @@ -581,7 +728,7 @@ def synchronize_model(): at the current simulation time """ function = LegacyFunctionSpecification() - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK @@ -597,7 +744,7 @@ def recommit_particles(): the script. """ function = LegacyFunctionSpecification() - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK Model is initialized and evolution can start @@ -614,10 +761,12 @@ def get_eps2(): """ function = LegacyFunctionSpecification() function.addParameter( - 'epsilon_squared', dtype='float64', direction=function.OUT, - description="The current value of the smooting parameter, squared." + "epsilon_squared", + dtype="float64", + direction=function.OUT, + description="The current value of the smooting parameter, squared.", ) - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the smoothing parameter was set @@ -633,10 +782,12 @@ def set_eps2(): """ function = LegacyFunctionSpecification() function.addParameter( - 'epsilon_squared', dtype='float64', direction=function.IN, - description="The new value of the smooting parameter, squared." + "epsilon_squared", + dtype="float64", + direction=function.IN, + description="The new value of the smooting parameter, squared.", ) - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the smoothing parameter was set @@ -652,10 +803,12 @@ def get_kinetic_energy(): """ function = LegacyFunctionSpecification() function.addParameter( - 'kinetic_energy', dtype='float64', direction=function.OUT, - description="The kinetic energy of the model" + "kinetic_energy", + dtype="float64", + direction=function.OUT, + description="The kinetic energy of the model", ) - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the kinetic energy was set @@ -671,10 +824,12 @@ def get_potential_energy(): """ function = LegacyFunctionSpecification() function.addParameter( - 'potential_energy', dtype='float64', direction=function.OUT, - description="The potential energy of the model" + "potential_energy", + dtype="float64", + direction=function.OUT, + description="The potential energy of the model", ) - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the potential energy was set @@ -692,9 +847,12 @@ def get_time(): """ function = LegacyFunctionSpecification() function.addParameter( - 'time', dtype='float64', direction=function.OUT, - description="The current model time") - function.result_type = 'int32' + "time", + dtype="float64", + direction=function.OUT, + description="The current model time", + ) + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the time was retrieved @@ -710,9 +868,13 @@ def get_begin_time(): """ function = LegacyFunctionSpecification() function.addParameter( - 'time', dtype='float64', direction=function.OUT, - description="The begin time", unit=nbody_system.time) - function.result_type = 'int32' + "time", + dtype="float64", + direction=function.OUT, + description="The begin time", + unit=nbody_system.time, + ) + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the time was retrieved @@ -729,10 +891,13 @@ def set_begin_time(): """ function = LegacyFunctionSpecification() function.addParameter( - 'time', dtype='float64', direction=function.IN, + "time", + dtype="float64", + direction=function.IN, description="The model time to start at", - unit=nbody_system.time) - function.result_type = 'int32' + unit=nbody_system.time, + ) + function.result_type = "int32" function.result_doc = """ 0 - OK Time value was changed @@ -748,10 +913,12 @@ def get_time_step(): """ function = LegacyFunctionSpecification() function.addParameter( - 'time_step', dtype='float64', direction=function.OUT, - description="The current model timestep" + "time_step", + dtype="float64", + direction=function.OUT, + description="The current model timestep", ) - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the time step was retrieved @@ -767,10 +934,12 @@ def get_total_mass(): """ function = LegacyFunctionSpecification() function.addParameter( - 'mass', dtype='float64', direction=function.OUT, - description="The total mass of the model" + "mass", + dtype="float64", + direction=function.OUT, + description="The total mass of the model", ) - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the kinetic mass was retrieved @@ -789,15 +958,24 @@ def get_center_of_mass_position(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'x', dtype='float64', direction=function.OUT, - description="The center of mass of the model") + "x", + dtype="float64", + direction=function.OUT, + description="The center of mass of the model", + ) function.addParameter( - 'y', dtype='float64', direction=function.OUT, - description="The center of mass of the model") + "y", + dtype="float64", + direction=function.OUT, + description="The center of mass of the model", + ) function.addParameter( - 'z', dtype='float64', direction=function.OUT, - description="The center of mass of the model") - function.result_type = 'int32' + "z", + dtype="float64", + direction=function.OUT, + description="The center of mass of the model", + ) + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the center was retrieved @@ -817,15 +995,24 @@ def get_center_of_mass_velocity(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'vx', dtype='float64', direction=function.OUT, - description="The mean velocity of the model") + "vx", + dtype="float64", + direction=function.OUT, + description="The mean velocity of the model", + ) function.addParameter( - 'vy', dtype='float64', direction=function.OUT, - description="The mean velocity of the model") + "vy", + dtype="float64", + direction=function.OUT, + description="The mean velocity of the model", + ) function.addParameter( - 'vz', dtype='float64', direction=function.OUT, - description="The mean velocity of the model") - function.result_type = 'int32' + "vz", + dtype="float64", + direction=function.OUT, + description="The mean velocity of the model", + ) + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the center of mass velocity was retrieved @@ -842,13 +1029,14 @@ def get_total_radius(): """ function = LegacyFunctionSpecification() function.addParameter( - 'radius', dtype='float64', direction=function.OUT, + "radius", + dtype="float64", + direction=function.OUT, description=( - "The maximum distance from a star to the center of mass of " - "the model" - ) + "The maximum distance from a star to the center of mass of the model" + ), ) - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK Current value of the radius was retrieved @@ -860,13 +1048,16 @@ def get_total_radius(): @legacy_function def get_number_of_particles(): """ - Retrieve the total number of particles define d in the code + Retrieve the total number of particles defined in the code """ function = LegacyFunctionSpecification() function.addParameter( - 'number_of_particles', dtype='int32', direction=function.OUT, - description="Count of the particles in the code") - function.result_type = 'int32' + "number_of_particles", + dtype="int32", + direction=function.OUT, + description="Count of the particles in the code", + ) + function.result_type = "int32" function.result_doc = """ 0 - OK Count could be determined @@ -891,9 +1082,12 @@ def get_index_of_first_particle(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.OUT, - description="Index of the first particle") - function.result_type = 'int32' + "index_of_the_particle", + dtype="int32", + direction=function.OUT, + description="Index of the first particle", + ) + function.result_type = "int32" function.result_doc = """ 0 - OK Index was set @@ -910,17 +1104,21 @@ def get_index_of_next_particle(): """ function = LegacyFunctionSpecification() function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, - description="Index of the particle") + "index_of_the_particle", + dtype="int32", + direction=function.IN, + description="Index of the particle", + ) function.addParameter( - 'index_of_the_next_particle', dtype='int32', + "index_of_the_next_particle", + dtype="int32", direction=function.OUT, description=( "Index of the particle following the particle with the " "provided index" - ) + ), ) - function.result_type = 'int32' + function.result_type = "int32" function.result_doc = """ 0 - OK Index was set @@ -932,123 +1130,7 @@ def get_index_of_next_particle(): return function -class GravityFieldInterface(object): - """ - Codes implementing the gravity field interface provide functions to - calculate the force and potential energy fields at any point. - """ - - @legacy_function - def get_gravity_at_point(): - """ - Get the gravitational acceleration at the given points. To calculate - the force on bodies at those points, multiply with the mass of the - bodies - """ - function = LegacyFunctionSpecification() - for x in ['eps', 'x', 'y', 'z']: - function.addParameter( - x, - dtype='float64', - direction=function.IN, - unit=nbody_system.length - ) - for x in ['ax', 'ay', 'az']: - function.addParameter( - x, - dtype='float64', - direction=function.OUT, - unit=nbody_system.acceleration - ) - function.addParameter('npoints', dtype='i', direction=function.LENGTH) - function.result_type = 'int32' - function.must_handle_array = True - return function - - @legacy_function - def get_potential_at_point(): - """ - Determine the gravitational potential on any given point - """ - function = LegacyFunctionSpecification() - for x in ['eps', 'x', 'y', 'z']: - function.addParameter( - x, - dtype='float64', - direction=function.IN, - unit=nbody_system.length - ) - for x in ['phi']: - function.addParameter( - x, - dtype='float64', - direction=function.OUT, - unit=nbody_system.potential - ) - function.addParameter('npoints', dtype='i', direction=function.LENGTH) - function.result_type = 'int32' - function.must_handle_array = True - return function - - -class SinglePointGravityFieldInterface(object): - """ - Codes implementing the gravity field interface provide functions to - calculate the force and potential energy fields at any point. - """ - - @legacy_function - def get_gravity_at_point(): - """ - Get the gravitational acceleration at the given points. To calculate - the force on bodies at those points, multiply with the mass of the - bodies - """ - function = LegacyFunctionSpecification() - for x in ['eps', 'x', 'y', 'z']: - function.addParameter( - x, - dtype='float64', - direction=function.IN, - unit=nbody_system.length - ) - for x in ['ax', 'ay', 'az']: - function.addParameter( - x, - dtype='float64', - direction=function.OUT, - unit=nbody_system.acceleration - ) - function.result_type = 'int32' - function.can_handle_array = True - return function - - @legacy_function - def get_potential_at_point(): - """ - Determine the gravitational potential on any given point - """ - function = LegacyFunctionSpecification() - for x in ['eps', 'x', 'y', 'z']: - function.addParameter( - x, - dtype='float64', - direction=function.IN, - unit=nbody_system.length - ) - for x in ['phi']: - function.addParameter( - x, - dtype='float64', - direction=function.OUT, - unit=nbody_system.potential - ) - function.result_type = 'int32' - function.can_handle_array = True - return function - - -class GravitationalDynamicsDocumentation(object): +class GravitationalDynamicsDocumentation: def __get__(self, instance, owner): @@ -1063,7 +1145,7 @@ class GravitationalDynamics(common.CommonCode): __doc__ = GravitationalDynamicsDocumentation() - def __init__(self, legacy_interface, unit_converter=None, **options): + def __init__(self, legacy_interface, unit_converter=None, **options): self.unit_converter = unit_converter common.CommonCode.__init__(self, legacy_interface, **options) @@ -1075,66 +1157,67 @@ def define_properties(self, handler): handler.add_property("get_center_of_mass_position") handler.add_property("get_center_of_mass_velocity") handler.add_property("get_total_mass") - handler.add_property('get_time', public_name="model_time") + handler.add_property("get_time", public_name="model_time") def define_state(self, handler): common.CommonCode.define_state(self, handler) - handler.add_transition('END', 'INITIALIZED', 'initialize_code', False) + handler.add_transition("END", "INITIALIZED", "initialize_code", False) - handler.add_transition('INITIALIZED', 'EDIT', 'commit_parameters') - handler.add_transition( - 'RUN', 'CHANGE_PARAMETERS_RUN', 'before_set_parameter', False) + handler.add_transition("INITIALIZED", "EDIT", "commit_parameters") handler.add_transition( - 'EDIT', 'CHANGE_PARAMETERS_EDIT', 'before_set_parameter', False) - handler.add_transition( - 'UPDATE', 'CHANGE_PARAMETERS_UPDATE', 'before_set_parameter', False + "RUN", "CHANGE_PARAMETERS_RUN", "before_set_parameter", False ) handler.add_transition( - 'CHANGE_PARAMETERS_RUN', 'RUN', 'recommit_parameters') + "EDIT", "CHANGE_PARAMETERS_EDIT", "before_set_parameter", False + ) handler.add_transition( - 'CHANGE_PARAMETERS_EDIT', 'EDIT', 'recommit_parameters') + "UPDATE", "CHANGE_PARAMETERS_UPDATE", "before_set_parameter", False + ) + handler.add_transition("CHANGE_PARAMETERS_RUN", "RUN", "recommit_parameters") + handler.add_transition("CHANGE_PARAMETERS_EDIT", "EDIT", "recommit_parameters") handler.add_transition( - 'CHANGE_PARAMETERS_UPDATE', 'UPDATE', 'recommit_parameters') - - handler.add_method('CHANGE_PARAMETERS_RUN', 'before_set_parameter') - handler.add_method('CHANGE_PARAMETERS_EDIT', 'before_set_parameter') - handler.add_method('CHANGE_PARAMETERS_UPDATE', 'before_set_parameter') - - handler.add_method('CHANGE_PARAMETERS_RUN', 'before_get_parameter') - handler.add_method('CHANGE_PARAMETERS_EDIT', 'before_get_parameter') - handler.add_method('CHANGE_PARAMETERS_UPDATE', 'before_get_parameter') - handler.add_method('RUN', 'before_get_parameter') - handler.add_method('EDIT', 'before_get_parameter') - handler.add_method('UPDATE', 'before_get_parameter') - handler.add_method('EVOLVED', 'before_get_parameter') - - handler.add_method('EDIT', 'new_particle') - handler.add_method('EDIT', 'delete_particle') - handler.add_method('UPDATE', 'new_particle') - handler.add_method('UPDATE', 'delete_particle') - handler.add_transition('EDIT', 'RUN', 'commit_particles') - handler.add_transition('RUN', 'UPDATE', 'new_particle', False) - handler.add_transition('RUN', 'UPDATE', 'delete_particle', False) - handler.add_transition('UPDATE', 'RUN', 'recommit_particles') - handler.add_transition('RUN', 'EVOLVED', 'evolve_model', False) - handler.add_method('EVOLVED', 'evolve_model') - handler.add_transition('EVOLVED', 'RUN', 'synchronize_model') - handler.add_method('RUN', 'synchronize_model') - handler.add_method('RUN', 'get_state') - handler.add_method('RUN', 'get_mass') - handler.add_method('RUN', 'get_position') - handler.add_method('RUN', 'get_velocity') - handler.add_method('RUN', 'get_potential') - handler.add_method('RUN', 'get_potential_energy') - handler.add_method('RUN', 'get_kinetic_energy') - handler.add_transition('RUN', 'UPDATE', 'set_mass', False) - handler.add_transition('RUN', 'UPDATE', 'set_position', False) - handler.add_transition('RUN', 'UPDATE', 'set_velocity', False) - handler.add_transition('RUN', 'UPDATE', 'set_radius', False) - handler.add_method('UPDATE', 'set_mass') - handler.add_method('UPDATE', 'set_position') - handler.add_method('UPDATE', 'set_velocity') - handler.add_method('UPDATE', 'set_radius') + "CHANGE_PARAMETERS_UPDATE", "UPDATE", "recommit_parameters" + ) + + handler.add_method("CHANGE_PARAMETERS_RUN", "before_set_parameter") + handler.add_method("CHANGE_PARAMETERS_EDIT", "before_set_parameter") + handler.add_method("CHANGE_PARAMETERS_UPDATE", "before_set_parameter") + + handler.add_method("CHANGE_PARAMETERS_RUN", "before_get_parameter") + handler.add_method("CHANGE_PARAMETERS_EDIT", "before_get_parameter") + handler.add_method("CHANGE_PARAMETERS_UPDATE", "before_get_parameter") + handler.add_method("RUN", "before_get_parameter") + handler.add_method("EDIT", "before_get_parameter") + handler.add_method("UPDATE", "before_get_parameter") + handler.add_method("EVOLVED", "before_get_parameter") + + handler.add_method("EDIT", "new_particle") + handler.add_method("EDIT", "delete_particle") + handler.add_method("UPDATE", "new_particle") + handler.add_method("UPDATE", "delete_particle") + handler.add_transition("EDIT", "RUN", "commit_particles") + handler.add_transition("RUN", "UPDATE", "new_particle", False) + handler.add_transition("RUN", "UPDATE", "delete_particle", False) + handler.add_transition("UPDATE", "RUN", "recommit_particles") + handler.add_transition("RUN", "EVOLVED", "evolve_model", False) + handler.add_method("EVOLVED", "evolve_model") + handler.add_transition("EVOLVED", "RUN", "synchronize_model") + handler.add_method("RUN", "synchronize_model") + handler.add_method("RUN", "get_state") + handler.add_method("RUN", "get_mass") + handler.add_method("RUN", "get_position") + handler.add_method("RUN", "get_velocity") + handler.add_method("RUN", "get_potential") + handler.add_method("RUN", "get_potential_energy") + handler.add_method("RUN", "get_kinetic_energy") + handler.add_transition("RUN", "UPDATE", "set_mass", False) + handler.add_transition("RUN", "UPDATE", "set_position", False) + handler.add_transition("RUN", "UPDATE", "set_velocity", False) + handler.add_transition("RUN", "UPDATE", "set_radius", False) + handler.add_method("UPDATE", "set_mass") + handler.add_method("UPDATE", "set_position") + handler.add_method("UPDATE", "set_velocity") + handler.add_method("UPDATE", "set_radius") def define_parameters(self, handler): handler.add_method_parameter( @@ -1142,7 +1225,7 @@ def define_parameters(self, handler): None, "timestep", "constant timestep for iteration", - default_value=0.7 | nbody_system.time + default_value=0.7 | nbody_system.time, ) handler.add_method_parameter( @@ -1150,20 +1233,12 @@ def define_parameters(self, handler): "set_begin_time", "begin_time", "model time to start the simulation at", - default_value=0.0 | nbody_system.time + default_value=0.0 | nbody_system.time, ) def define_methods(self, handler): common.CommonCode.define_methods(self, handler) - handler.add_method( - 'evolve_model', - ( - nbody_system.time, - ), - ( - handler.ERROR_CODE, - ) - ) + handler.add_method("evolve_model", (nbody_system.time,), (handler.ERROR_CODE,)) handler.add_method( "new_particle", @@ -1180,22 +1255,12 @@ def define_methods(self, handler): ( handler.INDEX, handler.ERROR_CODE, - ) - ) - handler.add_method( - "delete_particle", - ( - handler.NO_UNIT, ), - ( - handler.ERROR_CODE, - ) ) + handler.add_method("delete_particle", (handler.NO_UNIT,), (handler.ERROR_CODE,)) handler.add_method( "get_state", - ( - handler.NO_UNIT, - ), + (handler.NO_UNIT,), ( nbody_system.mass, nbody_system.length, @@ -1205,8 +1270,8 @@ def define_methods(self, handler): nbody_system.speed, nbody_system.speed, nbody_system.length, - handler.ERROR_CODE - ) + handler.ERROR_CODE, + ), ) handler.add_method( "set_state", @@ -1221,9 +1286,7 @@ def define_methods(self, handler): nbody_system.speed, nbody_system.length, ), - ( - handler.ERROR_CODE - ) + (handler.ERROR_CODE), ) handler.add_method( "set_mass", @@ -1231,19 +1294,10 @@ def define_methods(self, handler): handler.NO_UNIT, nbody_system.mass, ), - ( - handler.ERROR_CODE - ) + (handler.ERROR_CODE), ) handler.add_method( - "get_mass", - ( - handler.NO_UNIT, - ), - ( - nbody_system.mass, - handler.ERROR_CODE - ) + "get_mass", (handler.NO_UNIT,), (nbody_system.mass, handler.ERROR_CODE) ) handler.add_method( "set_radius", @@ -1251,19 +1305,10 @@ def define_methods(self, handler): handler.NO_UNIT, nbody_system.length, ), - ( - handler.ERROR_CODE - ) + (handler.ERROR_CODE), ) handler.add_method( - "get_radius", - ( - handler.NO_UNIT, - ), - ( - nbody_system.length, - handler.ERROR_CODE - ) + "get_radius", (handler.NO_UNIT,), (nbody_system.length, handler.ERROR_CODE) ) handler.add_method( "set_position", @@ -1273,21 +1318,17 @@ def define_methods(self, handler): nbody_system.length, nbody_system.length, ), - ( - handler.ERROR_CODE - ) + (handler.ERROR_CODE), ) handler.add_method( "get_position", - ( - handler.NO_UNIT, - ), + (handler.NO_UNIT,), ( nbody_system.length, nbody_system.length, nbody_system.length, - handler.ERROR_CODE - ) + handler.ERROR_CODE, + ), ) handler.add_method( "set_velocity", @@ -1297,79 +1338,62 @@ def define_methods(self, handler): nbody_system.speed, nbody_system.speed, ), - ( - handler.ERROR_CODE - ) + (handler.ERROR_CODE), ) handler.add_method( "get_velocity", - ( - handler.NO_UNIT, - ), + (handler.NO_UNIT,), ( nbody_system.speed, nbody_system.speed, nbody_system.speed, - handler.ERROR_CODE - ) + handler.ERROR_CODE, + ), ) handler.add_method( "get_acceleration", - ( - handler.NO_UNIT, - ), + (handler.NO_UNIT,), ( nbody_system.acceleration, nbody_system.acceleration, nbody_system.acceleration, - handler.ERROR_CODE - ) + handler.ERROR_CODE, + ), ) handler.add_method( "get_potential", - ( - handler.NO_UNIT, - ), + (handler.NO_UNIT,), ( nbody_system.length**2 * nbody_system.time**-2, handler.ERROR_CODE, - ) + ), ) handler.add_method( - 'get_indices_of_colliding_particles', + "get_indices_of_colliding_particles", (), ( handler.NO_UNIT, handler.NO_UNIT, handler.ERROR_CODE, - ) + ), ) - handler.add_method( - 'commit_particles', - (), - (handler.ERROR_CODE) - ) + handler.add_method("commit_particles", (), (handler.ERROR_CODE)) - handler.add_method( - 'recommit_particles', - (), - (handler.ERROR_CODE) - ) + handler.add_method("recommit_particles", (), (handler.ERROR_CODE)) - handler.add_method( - 'synchronize_model', - (), - (handler.ERROR_CODE) - ) + handler.add_method("synchronize_model", (), (handler.ERROR_CODE)) handler.add_method( "get_time_step", (), - (nbody_system.time, handler.ERROR_CODE,) + ( + nbody_system.time, + handler.ERROR_CODE, + ), ) handler.add_method( @@ -1378,7 +1402,7 @@ def define_methods(self, handler): ( nbody_system.mass * nbody_system.length**2 * nbody_system.time**-2, handler.ERROR_CODE, - ) + ), ) handler.add_method( @@ -1387,22 +1411,27 @@ def define_methods(self, handler): ( nbody_system.mass * nbody_system.length**2 * nbody_system.time**-2, handler.ERROR_CODE, - ) + ), ) handler.add_method( "get_total_radius", (), - (nbody_system.length, handler.ERROR_CODE,) + ( + nbody_system.length, + handler.ERROR_CODE, + ), ) handler.add_method( "get_center_of_mass_position", (), ( - nbody_system.length, nbody_system.length, nbody_system.length, + nbody_system.length, + nbody_system.length, + nbody_system.length, handler.ERROR_CODE, - ) + ), ) handler.add_method( @@ -1413,38 +1442,45 @@ def define_methods(self, handler): nbody_system.length / nbody_system.time, nbody_system.length / nbody_system.time, handler.ERROR_CODE, - ) + ), ) handler.add_method( "get_total_mass", (), - (nbody_system.mass, handler.ERROR_CODE,) + ( + nbody_system.mass, + handler.ERROR_CODE, + ), ) handler.add_method( - 'get_time', + "get_time", (), - (nbody_system.time, handler.ERROR_CODE,) + ( + nbody_system.time, + handler.ERROR_CODE, + ), ) def define_particle_sets(self, handler): - handler.define_set('particles', 'index_of_the_particle') - handler.set_new('particles', 'new_particle') - handler.set_delete('particles', 'delete_particle') - handler.add_setter('particles', 'set_state') - handler.add_getter('particles', 'get_state') - handler.add_setter('particles', 'set_mass') - handler.add_getter('particles', 'get_mass', names=('mass',)) - handler.add_setter('particles', 'set_position') - handler.add_getter('particles', 'get_position') - handler.add_setter('particles', 'set_velocity') - handler.add_getter('particles', 'get_velocity') - handler.add_setter('particles', 'set_radius') - handler.add_getter('particles', 'get_radius') + handler.define_set("particles", "index_of_the_particle") + handler.set_new("particles", "new_particle") + handler.set_delete("particles", "delete_particle") + handler.add_setter("particles", "set_state") + handler.add_getter("particles", "get_state") + handler.add_setter("particles", "set_mass") + handler.add_getter("particles", "get_mass", names=("mass",)) + handler.add_setter("particles", "set_position") + handler.add_getter("particles", "get_position") + handler.add_setter("particles", "set_velocity") + handler.add_getter("particles", "get_velocity") + handler.add_setter("particles", "set_radius") + handler.add_getter("particles", "get_radius") handler.add_query( - 'particles', 'get_indices_of_colliding_particles', - public_name='select_colliding_particles' + "particles", + "get_indices_of_colliding_particles", + public_name="select_colliding_particles", ) def get_colliding_particles(self): @@ -1453,9 +1489,7 @@ def get_colliding_particles(self): def define_converter(self, handler): if self.unit_converter is not None: - handler.set_converter( - self.unit_converter.as_converter_from_si_to_generic() - ) + handler.set_converter(self.unit_converter.as_converter_from_si_to_generic()) def commit_parameters(self): self.parameters.send_not_set_parameters_to_code() @@ -1465,7 +1499,7 @@ def commit_parameters(self): def cleanup_code(self): self.overridden().cleanup_code() - handler = self.get_handler('PARTICLES') + handler = self.get_handler("PARTICLES") handler._cleanup_instances() def reset(self): @@ -1476,10 +1510,3 @@ def reset(self): def get_total_energy(self): return self.get_potential_energy() + self.get_kinetic_energy() - - -class GravityFieldCode(object): - - def define_state(self, handler): - handler.add_method('RUN', 'get_gravity_at_point') - handler.add_method('RUN', 'get_potential_at_point') diff --git a/src/amuse/community/interface/gd/gravitational_dynamics_64.py b/src/amuse/community/interface/gd/gravitational_dynamics_64.py new file mode 100644 index 0000000000..26f32db684 --- /dev/null +++ b/src/amuse/community/interface/gd/gravitational_dynamics_64.py @@ -0,0 +1,1512 @@ +""" +Stellar Dynamics Interface Definition +""" + +from amuse.units import nbody_system +from amuse.community.interface import common + +from amuse.rfi.core import legacy_function +from amuse.rfi.core import LegacyFunctionSpecification + + +class GravitationalDynamics64Interface(common.CommonCodeInterface): + + @legacy_function + def new_particle(): + """ + Define a new particle in the stellar dynamics code. The particle is + initialized with the provided mass, radius, position and velocity. This + function returns an index that can be used to refer to this particle. + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.OUT, + description=( + "An index assigned to the newly created particle. " + "This index is supposed to be a local index for the code " + "(and not valid in other instances of the code or in other " + "codes)" + ), + ) + + function.addParameter( + "mass", + dtype="float64", + direction=function.IN, + description="The mass of the particle", + ) + function.addParameter( + "x", + dtype="float64", + direction=function.IN, + description="The initial position vector of the particle", + ) + function.addParameter( + "y", + dtype="float64", + direction=function.IN, + description="The initial position vector of the particle", + ) + function.addParameter( + "z", + dtype="float64", + direction=function.IN, + description="The initial position vector of the particle", + ) + function.addParameter( + "vx", + dtype="float64", + direction=function.IN, + description="The initial velocity vector of the particle", + ) + function.addParameter( + "vy", + dtype="float64", + direction=function.IN, + description="The initial velocity vector of the particle", + ) + function.addParameter( + "vz", + dtype="float64", + direction=function.IN, + description="The initial velocity vector of the particle", + ) + function.addParameter( + "radius", + dtype="float64", + direction=function.IN, + description="The radius of the particle", + default=0, + ) + function.result_type = "int32" + function.result_doc = """ 0 - OK + particle was created and added to the model + -1 - ERROR + particle could not be created""" + return function + + @legacy_function + def delete_particle(): + """ + Remove the definition of particle from the code. After calling this + function the particle is no longer part of the model evolution. It is + up to the code if the index will be reused. + This function is optional. + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle to be removed. This index must have " + "been returned by an earlier call to :meth:`new_particle`" + ), + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + particle was removed from the model + -1 - ERROR + particle could not be removed + -2 - ERROR + not yet implemented + """ + return function + + @legacy_function + def get_state(): + """ + Retrieve the current state of a particle. The *minimal* information of + a stellar dynamics particle (mass, radius, position and velocity) is + returned. + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle to get the state from. This index must " + "have been returned by an earlier call to :meth:`new_particle`" + ), + ) + function.addParameter( + "mass", + dtype="float64", + direction=function.OUT, + description="The current mass of the particle", + ) + function.addParameter( + "x", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) + function.addParameter( + "y", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) + function.addParameter( + "z", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) + function.addParameter( + "vx", + dtype="float64", + direction=function.OUT, + description="The current velocity vector of the particle", + ) + function.addParameter( + "vy", + dtype="float64", + direction=function.OUT, + description="The current velocity vector of the particle", + ) + function.addParameter( + "vz", + dtype="float64", + direction=function.OUT, + description="The current velocity vector of the particle", + ) + function.addParameter( + "radius", + dtype="float64", + direction=function.OUT, + description="The current radius of the particle", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + particle was removed from the model + -1 - ERROR + particle could not be found + """ + return function + + @legacy_function + def set_state(): + """ + Update the current state of a particle. The *minimal* information of a + stellar dynamics particle (mass, radius, position and velocity) is + updated. + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle for which the state is to be updated. " + "This index must have been returned by an earlier call to " + ":meth:`new_particle`" + ), + ) + function.addParameter( + "mass", + dtype="float64", + direction=function.IN, + description="The new mass of the particle", + ) + function.addParameter( + "x", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) + function.addParameter( + "y", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) + function.addParameter( + "z", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) + function.addParameter( + "vx", + dtype="float64", + direction=function.IN, + description="The new velocity vector of the particle", + ) + function.addParameter( + "vy", + dtype="float64", + direction=function.IN, + description="The new velocity vector of the particle", + ) + function.addParameter( + "vz", + dtype="float64", + direction=function.IN, + description="The new velocity vector of the particle", + ) + function.addParameter( + "radius", + dtype="float64", + direction=function.IN, + description="The new radius of the particle", + default=0, + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + particle was found in the model and the information was set + -1 - ERROR + particle could not be found + -2 - ERROR + code does not support updating of a particle + -3 - ERROR + not yet implemented + """ + return function + + @legacy_function + def get_mass(): + """ + Retrieve the mass of a particle. Mass is a scalar property of a + particle, this function has one OUT argument. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle to get the state from. This index must " + "have been returned by an earlier call to :meth:`new_particle`" + ), + ) + function.addParameter( + "mass", + dtype="float64", + direction=function.OUT, + description="The current mass of the particle", + ) + function.result_type = "int32" + function.can_handle_array = True + function.result_doc = """ + 0 - OK + particle was removed from the model + -1 - ERROR + particle could not be found + """ + return function + + @legacy_function + def set_mass(): + """ + Update the mass of a particle. Mass is a scalar property of a particle. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle for which the state is to be updated. " + "This index must have been returned by an earlier call to " + ":meth:`new_particle`" + ), + ) + function.addParameter( + "mass", + dtype="float64", + direction=function.IN, + description="The new mass of the particle", + ) + function.result_type = "int32" + function.can_handle_array = True + function.result_doc = """ + 0 - OK + particle was found in the model and the information was set + -1 - ERROR + particle could not be found + -2 - ERROR + code does not support updating of a particle + """ + return function + + @legacy_function + def get_radius(): + """ + Retrieve the radius of a particle. Radius is a scalar property of a + particle, this function has one OUT argument. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle to get the radius of. This index must " + "have been returned by an earlier call to :meth:`new_particle`" + ), + ) + function.addParameter( + "radius", + dtype="float64", + direction=function.OUT, + description="The current radius of the particle", + ) + function.result_type = "int32" + function.can_handle_array = True + function.result_doc = """ + 0 - OK + particle was found in the model and the information was retreived + -1 - ERROR + particle could not be found + """ + return function + + @legacy_function + def set_radius(): + """ + Set the radius of a particle. Radius is a scalar property of a + particle. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle to get the radius of. This index must " + "have been returned by an earlier call to :meth:`new_particle`" + ), + ) + function.addParameter( + "radius", + dtype="float64", + direction=function.IN, + description="The new radius of the particle", + ) + function.result_type = "int32" + function.can_handle_array = True + function.result_doc = """ + 0 - OK + particle was found in the model and the information was retreived + -1 - ERROR + particle could not be found + """ + return function + + @legacy_function + def get_position(): + """ + Retrieve the position vector of a particle. Position is a vector + property, this function has 3 OUT arguments. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle to get the state from. This index must " + "have been returned by an earlier call to :meth:`new_particle`" + ), + ) + function.addParameter( + "x", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) + function.addParameter( + "y", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) + function.addParameter( + "z", + dtype="float64", + direction=function.OUT, + description="The current position vector of the particle", + ) + function.result_type = "int32" + function.can_handle_array = True + function.result_doc = """ + 0 - OK + current value was retrieved + -1 - ERROR + particle could not be found + -2 - ERROR + not yet implemented + """ + return function + + @legacy_function + def set_position(): + """ + Update the position of a particle. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle for which the state is to be updated. " + "This index must have been returned by an earlier call to " + ":meth:`new_particle`" + ), + ) + function.addParameter( + "x", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) + function.addParameter( + "y", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) + function.addParameter( + "z", + dtype="float64", + direction=function.IN, + description="The new position vector of the particle", + ) + function.result_type = "int32" + function.can_handle_array = True + function.result_doc = """ + 0 - OK + particle was found in the model and the information was set + -1 - ERROR + particle could not be found + -2 - ERROR + code does not support updating of a particle + """ + return function + + @legacy_function + def get_velocity(): + """ + Retrieve the velocity vector of a particle. Position is a vector + property, this function has 3 OUT arguments. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle to get the velocity from. This index " + "must have been returned by an earlier call to " + ":meth:`new_particle`" + ), + ) + function.addParameter( + "vx", + dtype="float64", + direction=function.OUT, + description=( + "The current x component of the position vector of the " "particle" + ), + ) + function.addParameter( + "vy", + dtype="float64", + direction=function.OUT, + description=( + "The current y component of the position vector of the " "particle" + ), + ) + function.addParameter( + "vz", + dtype="float64", + direction=function.OUT, + description=( + "The current z component of the position vector of the " "particle" + ), + ) + function.result_type = "int32" + function.can_handle_array = True + function.result_doc = """ + 0 - OK + current value was retrieved + -1 - ERROR + particle could not be found + -2 - ERROR + not yet implemented + """ + return function + + @legacy_function + def set_velocity(): + """ + Set the velocity vector of a particle. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle to get the state from. This index must " + "have been returned by an earlier call to :meth:`new_particle`" + ), + ) + function.addParameter( + "vx", + dtype="float64", + direction=function.IN, + description=( + "The current x component of the velocity vector of the " "particle" + ), + ) + function.addParameter( + "vy", + dtype="float64", + direction=function.IN, + description=( + "The current y component of the velocity vector of the " "particle" + ), + ) + function.addParameter( + "vz", + dtype="float64", + direction=function.IN, + description=( + "The current z component of the velocity vector of the " "particle" + ), + ) + function.result_type = "int32" + function.can_handle_array = True + function.result_doc = """ + 0 - OK + current value was retrieved + -1 - ERROR + particle could not be found + -2 - ERROR + not yet implemented + """ + return function + + @legacy_function + def get_acceleration(): + """ + Retrieve the acceleration vector of a particle. Second time derivative + of the position. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle to get the state from. This index must " + "have been returned by an earlier call to :meth:`new_particle`" + ), + ) + function.addParameter( + "ax", + dtype="float64", + direction=function.OUT, + description="The current acceleration vector of the particle", + ) + function.addParameter( + "ay", + dtype="float64", + direction=function.OUT, + description="The current acceleration vector of the particle", + ) + function.addParameter( + "az", + dtype="float64", + direction=function.OUT, + description="The current acceleration vector of the particle", + ) + function.result_type = "int32" + function.can_handle_array = True + function.result_doc = """ + 0 - OK + current value was retrieved + -1 - ERROR + particle could not be found + -2 - ERROR + not yet implemented + """ + return function + + @legacy_function + def get_potential(): + """ + Retrieve the potential at a particle position, for retrieving the + potential anywhere in the field use get_potential_at_point. + """ + function = LegacyFunctionSpecification() + + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description=( + "Index of the particle for which the state is to be updated. " + "This index must have been returned by an earlier call to " + ":meth:`new_particle`" + ), + ) + function.addParameter( + "potential", + dtype="float64", + direction=function.OUT, + description="The current scalar potential...", + ) + function.result_type = "int32" + function.can_handle_array = True + function.result_doc = """ + 0 - OK + current value was retrieved + -1 - ERROR + particle could not be found + -2 - ERROR + not yet implemented + """ + return function + + @legacy_function + def evolve_model(): + """ + Evolve the model until the given time, or until a stopping condition is + set. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "time", + dtype="float64", + direction=function.IN, + description=( + "Model time to evolve the code to. The model will be " + "evolved until this time is reached exactly or just after." + ), + ) + function.result_type = "int32" + return function + + @legacy_function + def commit_particles(): + """ + Let the code perform initialization actions + after all particles have been loaded in the model. + Should be called before the first evolve call and + after the last new_particle call. + """ + function = LegacyFunctionSpecification() + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Model is initialized and evolution can start + -1 - ERROR + Error happened during initialization, this error needs to be + further specified by every code implemention + """ + return function + + @legacy_function + def synchronize_model(): + """ + After an evolve the particles may be at different simulation + times. Synchronize the particles to a consistent stat + at the current simulation time + """ + function = LegacyFunctionSpecification() + function.result_type = "int32" + function.result_doc = """ + 0 - OK + + """ + return function + + @legacy_function + def recommit_particles(): + """ + Let the code perform initialization actions + after the number of particles have been updated + or particle attributes have been updated from + the script. + """ + function = LegacyFunctionSpecification() + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Model is initialized and evolution can start + -1 - ERROR + Error happened during initialization, this error needs to be + further specified by every code implemention + """ + return function + + @legacy_function + def get_eps2(): + """ + Retrieve the current value of the squared smoothing parameter. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "epsilon_squared", + dtype="float64", + direction=function.OUT, + description="The current value of the smooting parameter, squared.", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the smoothing parameter was set + -1 - ERROR + The code does not have support for a smoothing parameter + """ + return function + + @legacy_function + def set_eps2(): + """ + Update the value of the squared smoothing parameter. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "epsilon_squared", + dtype="float64", + direction=function.IN, + description="The new value of the smooting parameter, squared.", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the smoothing parameter was set + -1 - ERROR + The code does not have support for a smoothing parameter + """ + return function + + @legacy_function + def get_kinetic_energy(): + """ + Retrieve the current kinetic energy of the model + """ + function = LegacyFunctionSpecification() + function.addParameter( + "kinetic_energy", + dtype="float64", + direction=function.OUT, + description="The kinetic energy of the model", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the kinetic energy was set + -1 - ERROR + Kinetic energy could not be provided + """ + return function + + @legacy_function + def get_potential_energy(): + """ + Retrieve the current potential energy of the model + """ + function = LegacyFunctionSpecification() + function.addParameter( + "potential_energy", + dtype="float64", + direction=function.OUT, + description="The potential energy of the model", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the potential energy was set + -1 - ERROR + Kinetic potential could not be provided + """ + return function + + @legacy_function + def get_time(): + """ + Retrieve the model time. This time should be close to the end time + specified in the evolve code. Or, when a collision was detected, it + will be the model time of the collision. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "time", + dtype="float64", + direction=function.OUT, + description="The current model time", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the time was retrieved + -1 - ERROR + The code does not have support for querying the time + """ + return function + + @legacy_function + def get_begin_time(): + """ + Retrieve the model time to start the evolution at. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "time", + dtype="float64", + direction=function.OUT, + description="The begin time", + unit=nbody_system.time, + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the time was retrieved + -2 - ERROR + The code does not have support for querying the begin time + """ + return function + + @legacy_function + def set_begin_time(): + """ + Set the model time to start the evolution at. This is an offset for + all further calculations in the code. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "time", + dtype="float64", + direction=function.IN, + description="The model time to start at", + unit=nbody_system.time, + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Time value was changed + -2 - ERROR + The code does not support setting the begin time + """ + return function + + @legacy_function + def get_time_step(): + """ + Retrieve the model timestep. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "time_step", + dtype="float64", + direction=function.OUT, + description="The current model timestep", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the time step was retrieved + -1 - ERROR + The code does not have support for querying the time + """ + return function + + @legacy_function + def get_total_mass(): + """ + Retrieve the sum of the masses of all particles. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "mass", + dtype="float64", + direction=function.OUT, + description="The total mass of the model", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the kinetic mass was retrieved + -1 - ERROR + Total mass could not be provided + -2 - ERROR + not yet implemented + """ + return function + + @legacy_function + def get_center_of_mass_position(): + """ + Retrieve the center of mass (a point in space) of all particles. + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + "x", + dtype="float64", + direction=function.OUT, + description="The center of mass of the model", + ) + function.addParameter( + "y", + dtype="float64", + direction=function.OUT, + description="The center of mass of the model", + ) + function.addParameter( + "z", + dtype="float64", + direction=function.OUT, + description="The center of mass of the model", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the center was retrieved + -1 - ERROR + The mass could not be provided + -2 - ERROR + not yet implemented + """ + return function + + @legacy_function + def get_center_of_mass_velocity(): + """ + Retrieve the velocity of the center of mass of all particles. This + velocity is mass weighted mean of the velocity of all particles. + """ + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + "vx", + dtype="float64", + direction=function.OUT, + description="The mean velocity of the model", + ) + function.addParameter( + "vy", + dtype="float64", + direction=function.OUT, + description="The mean velocity of the model", + ) + function.addParameter( + "vz", + dtype="float64", + direction=function.OUT, + description="The mean velocity of the model", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the center of mass velocity was retrieved + -1 - ERROR + The value could not be provided + """ + return function + + @legacy_function + def get_total_radius(): + """ + Return the radius of the sphere, centered on the center of mass that + contains all the particles. *get_size?* + """ + function = LegacyFunctionSpecification() + function.addParameter( + "radius", + dtype="float64", + direction=function.OUT, + description=( + "The maximum distance from a star to the center of mass of the model" + ), + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Current value of the radius was retrieved + -1 - ERROR + The value could not be provided + """ + return function + + @legacy_function + def get_number_of_particles(): + """ + Retrieve the total number of particles defined in the code + """ + function = LegacyFunctionSpecification() + function.addParameter( + "number_of_particles", + dtype="int32", + direction=function.OUT, + description="Count of the particles in the code", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Count could be determined + -1 - ERROR + Unable to determine the count + """ + return function + + @legacy_function + def get_index_of_first_particle(): + """ + Retrieve the index of first particle. When this index is used as the + starting index for the ``get_index_of_next_particle`` method, all + particles can be iterated over:: + + error, first_index = instance.get_index_of_first_particle() + current_index = first_index + while error == 0: + status, mass = instance.get_mass(current_index) + print mass + error, current_index = instance.get_index_of_next_particle(current_index) + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.OUT, + description="Index of the first particle", + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Index was set + 1 - ERROR + Code has no particles, or cannot set the index + """ + return function + + @legacy_function + def get_index_of_next_particle(): + """ + Retrieve the index of the particle following the provided index. The + iteration direction is determined by the code. + """ + function = LegacyFunctionSpecification() + function.addParameter( + "index_of_the_particle", + dtype="int64", + direction=function.IN, + description="Index of the particle", + ) + function.addParameter( + "index_of_the_next_particle", + dtype="int64", + direction=function.OUT, + description=( + "Index of the particle following the particle with the " + "provided index" + ), + ) + function.result_type = "int32" + function.result_doc = """ + 0 - OK + Index was set + 1 - STATE + Last index was reached + -1 - ERROR + Particle could not be found + """ + return function + + +class GravitationalDynamics64Documentation: + + def __get__(self, instance, owner): + + string = "" + + string = instance.parameters.__doc__ + return string + + +class GravitationalDynamics64(common.CommonCode): + NBODY = object() + + __doc__ = GravitationalDynamics64Documentation() + + def __init__(self, legacy_interface, unit_converter=None, **options): + self.unit_converter = unit_converter + + common.CommonCode.__init__(self, legacy_interface, **options) + + def define_properties(self, handler): + handler.add_property("get_kinetic_energy") + handler.add_property("get_potential_energy") + handler.add_property("get_total_radius") + handler.add_property("get_center_of_mass_position") + handler.add_property("get_center_of_mass_velocity") + handler.add_property("get_total_mass") + handler.add_property("get_time", public_name="model_time") + + def define_state(self, handler): + common.CommonCode.define_state(self, handler) + handler.add_transition("END", "INITIALIZED", "initialize_code", False) + + handler.add_transition("INITIALIZED", "EDIT", "commit_parameters") + handler.add_transition( + "RUN", "CHANGE_PARAMETERS_RUN", "before_set_parameter", False + ) + handler.add_transition( + "EDIT", "CHANGE_PARAMETERS_EDIT", "before_set_parameter", False + ) + handler.add_transition( + "UPDATE", "CHANGE_PARAMETERS_UPDATE", "before_set_parameter", False + ) + handler.add_transition("CHANGE_PARAMETERS_RUN", "RUN", "recommit_parameters") + handler.add_transition("CHANGE_PARAMETERS_EDIT", "EDIT", "recommit_parameters") + handler.add_transition( + "CHANGE_PARAMETERS_UPDATE", "UPDATE", "recommit_parameters" + ) + + handler.add_method("CHANGE_PARAMETERS_RUN", "before_set_parameter") + handler.add_method("CHANGE_PARAMETERS_EDIT", "before_set_parameter") + handler.add_method("CHANGE_PARAMETERS_UPDATE", "before_set_parameter") + + handler.add_method("CHANGE_PARAMETERS_RUN", "before_get_parameter") + handler.add_method("CHANGE_PARAMETERS_EDIT", "before_get_parameter") + handler.add_method("CHANGE_PARAMETERS_UPDATE", "before_get_parameter") + handler.add_method("RUN", "before_get_parameter") + handler.add_method("EDIT", "before_get_parameter") + handler.add_method("UPDATE", "before_get_parameter") + handler.add_method("EVOLVED", "before_get_parameter") + + handler.add_method("EDIT", "new_particle") + handler.add_method("EDIT", "delete_particle") + handler.add_method("UPDATE", "new_particle") + handler.add_method("UPDATE", "delete_particle") + handler.add_transition("EDIT", "RUN", "commit_particles") + handler.add_transition("RUN", "UPDATE", "new_particle", False) + handler.add_transition("RUN", "UPDATE", "delete_particle", False) + handler.add_transition("UPDATE", "RUN", "recommit_particles") + handler.add_transition("RUN", "EVOLVED", "evolve_model", False) + handler.add_method("EVOLVED", "evolve_model") + handler.add_transition("EVOLVED", "RUN", "synchronize_model") + handler.add_method("RUN", "synchronize_model") + handler.add_method("RUN", "get_state") + handler.add_method("RUN", "get_mass") + handler.add_method("RUN", "get_position") + handler.add_method("RUN", "get_velocity") + handler.add_method("RUN", "get_potential") + handler.add_method("RUN", "get_potential_energy") + handler.add_method("RUN", "get_kinetic_energy") + handler.add_transition("RUN", "UPDATE", "set_mass", False) + handler.add_transition("RUN", "UPDATE", "set_position", False) + handler.add_transition("RUN", "UPDATE", "set_velocity", False) + handler.add_transition("RUN", "UPDATE", "set_radius", False) + handler.add_method("UPDATE", "set_mass") + handler.add_method("UPDATE", "set_position") + handler.add_method("UPDATE", "set_velocity") + handler.add_method("UPDATE", "set_radius") + + def define_parameters(self, handler): + handler.add_method_parameter( + "get_time_step", + None, + "timestep", + "constant timestep for iteration", + default_value=0.7 | nbody_system.time, + ) + + handler.add_method_parameter( + "get_begin_time", + "set_begin_time", + "begin_time", + "model time to start the simulation at", + default_value=0.0 | nbody_system.time, + ) + + def define_methods(self, handler): + common.CommonCode.define_methods(self, handler) + handler.add_method("evolve_model", (nbody_system.time,), (handler.ERROR_CODE,)) + + handler.add_method( + "new_particle", + ( + nbody_system.mass, + nbody_system.length, + nbody_system.length, + nbody_system.length, + nbody_system.speed, + nbody_system.speed, + nbody_system.speed, + nbody_system.length, + ), + ( + handler.INDEX, + handler.ERROR_CODE, + ), + ) + handler.add_method("delete_particle", (handler.NO_UNIT,), (handler.ERROR_CODE,)) + handler.add_method( + "get_state", + (handler.NO_UNIT,), + ( + nbody_system.mass, + nbody_system.length, + nbody_system.length, + nbody_system.length, + nbody_system.speed, + nbody_system.speed, + nbody_system.speed, + nbody_system.length, + handler.ERROR_CODE, + ), + ) + handler.add_method( + "set_state", + ( + handler.NO_UNIT, + nbody_system.mass, + nbody_system.length, + nbody_system.length, + nbody_system.length, + nbody_system.speed, + nbody_system.speed, + nbody_system.speed, + nbody_system.length, + ), + (handler.ERROR_CODE), + ) + handler.add_method( + "set_mass", + ( + handler.NO_UNIT, + nbody_system.mass, + ), + (handler.ERROR_CODE), + ) + handler.add_method( + "get_mass", (handler.NO_UNIT,), (nbody_system.mass, handler.ERROR_CODE) + ) + handler.add_method( + "set_radius", + ( + handler.NO_UNIT, + nbody_system.length, + ), + (handler.ERROR_CODE), + ) + handler.add_method( + "get_radius", (handler.NO_UNIT,), (nbody_system.length, handler.ERROR_CODE) + ) + handler.add_method( + "set_position", + ( + handler.NO_UNIT, + nbody_system.length, + nbody_system.length, + nbody_system.length, + ), + (handler.ERROR_CODE), + ) + handler.add_method( + "get_position", + (handler.NO_UNIT,), + ( + nbody_system.length, + nbody_system.length, + nbody_system.length, + handler.ERROR_CODE, + ), + ) + handler.add_method( + "set_velocity", + ( + handler.NO_UNIT, + nbody_system.speed, + nbody_system.speed, + nbody_system.speed, + ), + (handler.ERROR_CODE), + ) + handler.add_method( + "get_velocity", + (handler.NO_UNIT,), + ( + nbody_system.speed, + nbody_system.speed, + nbody_system.speed, + handler.ERROR_CODE, + ), + ) + + handler.add_method( + "get_acceleration", + (handler.NO_UNIT,), + ( + nbody_system.acceleration, + nbody_system.acceleration, + nbody_system.acceleration, + handler.ERROR_CODE, + ), + ) + + handler.add_method( + "get_potential", + (handler.NO_UNIT,), + ( + nbody_system.length**2 * nbody_system.time**-2, + handler.ERROR_CODE, + ), + ) + + handler.add_method( + "get_indices_of_colliding_particles", + (), + ( + handler.NO_UNIT, + handler.NO_UNIT, + handler.ERROR_CODE, + ), + ) + + handler.add_method("commit_particles", (), (handler.ERROR_CODE)) + + handler.add_method("recommit_particles", (), (handler.ERROR_CODE)) + + handler.add_method("synchronize_model", (), (handler.ERROR_CODE)) + + handler.add_method( + "get_time_step", + (), + ( + nbody_system.time, + handler.ERROR_CODE, + ), + ) + + handler.add_method( + "get_kinetic_energy", + (), + ( + nbody_system.mass * nbody_system.length**2 * nbody_system.time**-2, + handler.ERROR_CODE, + ), + ) + + handler.add_method( + "get_potential_energy", + (), + ( + nbody_system.mass * nbody_system.length**2 * nbody_system.time**-2, + handler.ERROR_CODE, + ), + ) + + handler.add_method( + "get_total_radius", + (), + ( + nbody_system.length, + handler.ERROR_CODE, + ), + ) + + handler.add_method( + "get_center_of_mass_position", + (), + ( + nbody_system.length, + nbody_system.length, + nbody_system.length, + handler.ERROR_CODE, + ), + ) + + handler.add_method( + "get_center_of_mass_velocity", + (), + ( + nbody_system.length / nbody_system.time, + nbody_system.length / nbody_system.time, + nbody_system.length / nbody_system.time, + handler.ERROR_CODE, + ), + ) + + handler.add_method( + "get_total_mass", + (), + ( + nbody_system.mass, + handler.ERROR_CODE, + ), + ) + + handler.add_method( + "get_time", + (), + ( + nbody_system.time, + handler.ERROR_CODE, + ), + ) + + def define_particle_sets(self, handler): + handler.define_set("particles", "index_of_the_particle") + handler.set_new("particles", "new_particle") + handler.set_delete("particles", "delete_particle") + handler.add_setter("particles", "set_state") + handler.add_getter("particles", "get_state") + handler.add_setter("particles", "set_mass") + handler.add_getter("particles", "get_mass", names=("mass",)) + handler.add_setter("particles", "set_position") + handler.add_getter("particles", "get_position") + handler.add_setter("particles", "set_velocity") + handler.add_getter("particles", "get_velocity") + handler.add_setter("particles", "set_radius") + handler.add_getter("particles", "get_radius") + handler.add_query( + "particles", + "get_indices_of_colliding_particles", + public_name="select_colliding_particles", + ) + + def get_colliding_particles(self): + subset = self.colliding_particles_method._run(self, self.particles) + return subset + + def define_converter(self, handler): + if self.unit_converter is not None: + handler.set_converter(self.unit_converter.as_converter_from_si_to_generic()) + + def commit_parameters(self): + self.parameters.send_not_set_parameters_to_code() + self.parameters.send_cached_parameters_to_code() + self.overridden().commit_parameters() + + def cleanup_code(self): + self.overridden().cleanup_code() + + handler = self.get_handler("PARTICLES") + handler._cleanup_instances() + + def reset(self): + parameters = self.parameters.copy() + self.cleanup_code() + self.initialize_code() + self.parameters.reset_from_memento(parameters) + + def get_total_energy(self): + return self.get_potential_energy() + self.get_kinetic_energy() diff --git a/src/amuse/community/interface/gd/gravity_field.py b/src/amuse/community/interface/gd/gravity_field.py new file mode 100644 index 0000000000..9838148284 --- /dev/null +++ b/src/amuse/community/interface/gd/gravity_field.py @@ -0,0 +1,113 @@ +""" +Stellar Dynamics Interface Definition +""" + +from amuse.units import nbody_system + +from amuse.rfi.core import legacy_function +from amuse.rfi.core import LegacyFunctionSpecification + + +class SinglePointGravityFieldInterface: + """ + Codes implementing the gravity field interface provide functions to + calculate the force and potential energy fields at any point. + """ + + @legacy_function + def get_gravity_at_point(): + """ + Get the gravitational acceleration at the given points. To calculate + the force on bodies at those points, multiply with the mass of the + bodies + """ + function = LegacyFunctionSpecification() + for x in ["eps", "x", "y", "z"]: + function.addParameter( + x, dtype="float64", direction=function.IN, unit=nbody_system.length + ) + for x in ["ax", "ay", "az"]: + function.addParameter( + x, + dtype="float64", + direction=function.OUT, + unit=nbody_system.acceleration, + ) + function.result_type = "int32" + function.can_handle_array = True + return function + + @legacy_function + def get_potential_at_point(): + """ + Determine the gravitational potential on any given point + """ + function = LegacyFunctionSpecification() + for x in ["eps", "x", "y", "z"]: + function.addParameter( + x, dtype="float64", direction=function.IN, unit=nbody_system.length + ) + for x in ["phi"]: + function.addParameter( + x, dtype="float64", direction=function.OUT, unit=nbody_system.potential + ) + function.result_type = "int32" + function.can_handle_array = True + return function + + +class GravityFieldInterface: + """ + Codes implementing the gravity field interface provide functions to + calculate the force and potential energy fields at any point. + """ + + @legacy_function + def get_gravity_at_point(): + """ + Get the gravitational acceleration at the given points. To calculate + the force on bodies at those points, multiply with the mass of the + bodies + """ + function = LegacyFunctionSpecification() + for x in ["eps", "x", "y", "z"]: + function.addParameter( + x, dtype="float64", direction=function.IN, unit=nbody_system.length + ) + for x in ["ax", "ay", "az"]: + function.addParameter( + x, + dtype="float64", + direction=function.OUT, + unit=nbody_system.acceleration, + ) + function.addParameter("npoints", dtype="i", direction=function.LENGTH) + function.result_type = "int32" + function.must_handle_array = True + return function + + @legacy_function + def get_potential_at_point(): + """ + Determine the gravitational potential on any given point + """ + function = LegacyFunctionSpecification() + for x in ["eps", "x", "y", "z"]: + function.addParameter( + x, dtype="float64", direction=function.IN, unit=nbody_system.length + ) + for x in ["phi"]: + function.addParameter( + x, dtype="float64", direction=function.OUT, unit=nbody_system.potential + ) + function.addParameter("npoints", dtype="i", direction=function.LENGTH) + function.result_type = "int32" + function.must_handle_array = True + return function + + +class GravityFieldCode: + + def define_state(self, handler): + handler.add_method("RUN", "get_gravity_at_point") + handler.add_method("RUN", "get_potential_at_point") diff --git a/src/amuse/community/phantom/interface.f90 b/src/amuse/community/phantom/interface.f90 index dc94961ad5..9767f234e7 100644 --- a/src/amuse/community/phantom/interface.f90 +++ b/src/amuse/community/phantom/interface.f90 @@ -58,7 +58,7 @@ function get_time(time) function get_mass(index_of_the_particle, mass) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass integer :: get_mass call amuse_get_mass(index_of_the_particle, mass) @@ -67,7 +67,7 @@ function get_mass(index_of_the_particle, mass) function set_mass(index_of_the_particle, mass) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass integer :: set_mass call amuse_set_mass(index_of_the_particle, mass) @@ -76,7 +76,7 @@ function set_mass(index_of_the_particle, mass) function set_smoothing_length(index_of_the_particle, h_smooth) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: h_smooth integer :: set_smoothing_length call amuse_set_smoothing_length(index_of_the_particle, h_smooth) @@ -85,7 +85,7 @@ function set_smoothing_length(index_of_the_particle, h_smooth) function set_internal_energy(index_of_the_particle, u) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: u integer :: set_internal_energy call amuse_set_internal_energy(index_of_the_particle, u) @@ -94,7 +94,7 @@ function set_internal_energy(index_of_the_particle, u) function set_h2ratio(index_of_the_particle, h2ratio) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: h2ratio integer :: set_h2ratio call amuse_set_h2ratio(index_of_the_particle, h2ratio) @@ -103,7 +103,7 @@ function set_h2ratio(index_of_the_particle, h2ratio) function set_hi_abundance(index_of_the_particle, hi_abundance) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: hi_abundance integer :: set_hi_abundance call amuse_set_hi_abundance(index_of_the_particle, hi_abundance) @@ -112,7 +112,7 @@ function set_hi_abundance(index_of_the_particle, hi_abundance) function set_proton_abundance(index_of_the_particle, proton_abundance) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: proton_abundance integer :: set_proton_abundance call amuse_set_proton_abundance(index_of_the_particle, proton_abundance) @@ -121,7 +121,7 @@ function set_proton_abundance(index_of_the_particle, proton_abundance) function set_electron_abundance(index_of_the_particle, electron_abundance) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: electron_abundance integer :: set_electron_abundance call amuse_set_electron_abundance(index_of_the_particle, electron_abundance) @@ -130,7 +130,7 @@ function set_electron_abundance(index_of_the_particle, electron_abundance) function set_co_abundance(index_of_the_particle, co_abundance) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: co_abundance integer :: set_co_abundance call amuse_set_co_abundance(index_of_the_particle, co_abundance) @@ -139,7 +139,7 @@ function set_co_abundance(index_of_the_particle, co_abundance) function get_h2ratio(index_of_the_particle, h2ratio) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: h2ratio integer :: get_h2ratio call amuse_get_h2ratio(index_of_the_particle, h2ratio) @@ -148,7 +148,7 @@ function get_h2ratio(index_of_the_particle, h2ratio) function get_hi_abundance(index_of_the_particle, hi_abundance) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: hi_abundance integer :: get_hi_abundance call amuse_get_hi_abundance(index_of_the_particle, hi_abundance) @@ -157,7 +157,7 @@ function get_hi_abundance(index_of_the_particle, hi_abundance) function get_proton_abundance(index_of_the_particle, proton_abundance) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: proton_abundance integer :: get_proton_abundance call amuse_get_proton_abundance(index_of_the_particle, proton_abundance) @@ -166,7 +166,7 @@ function get_proton_abundance(index_of_the_particle, proton_abundance) function get_electron_abundance(index_of_the_particle, electron_abundance) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: electron_abundance integer :: get_electron_abundance call amuse_get_electron_abundance(index_of_the_particle, electron_abundance) @@ -175,7 +175,7 @@ function get_electron_abundance(index_of_the_particle, electron_abundance) function get_co_abundance(index_of_the_particle, co_abundance) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: co_abundance integer :: get_co_abundance call amuse_get_co_abundance(index_of_the_particle, co_abundance) @@ -184,7 +184,7 @@ function get_co_abundance(index_of_the_particle, co_abundance) function get_pressure(index_of_the_particle, p) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: p integer :: get_pressure call amuse_get_pressure(index_of_the_particle, p) @@ -193,7 +193,7 @@ function get_pressure(index_of_the_particle, p) function get_density(index_of_the_particle, rho) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: rho integer :: get_density call amuse_get_density(index_of_the_particle, rho) @@ -202,7 +202,7 @@ function get_density(index_of_the_particle, rho) function get_smoothing_length(index_of_the_particle, h_smooth) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: h_smooth integer :: get_smoothing_length call amuse_get_smoothing_length(index_of_the_particle, h_smooth) @@ -211,7 +211,7 @@ function get_smoothing_length(index_of_the_particle, h_smooth) function get_internal_energy(index_of_the_particle, u) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: u integer :: get_internal_energy call amuse_get_internal_energy(index_of_the_particle, u) @@ -220,7 +220,7 @@ function get_internal_energy(index_of_the_particle, u) function get_particle_type(index_of_the_particle, itype) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle integer :: itype integer :: get_particle_type !call amuse_get_particle_type(index_of_the_particle, itype) @@ -229,7 +229,7 @@ function get_particle_type(index_of_the_particle, itype) function get_index_of_first_particle(index_of_the_particle) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle integer :: get_index_of_first_particle get_index_of_first_particle=0 end function @@ -262,7 +262,8 @@ function evolve_model(tmax) double precision :: tmax integer :: evolve_model integer :: sc - integer :: i, nmax + integer(kind=8) :: i + integer :: nmax integer :: is_density_limit_detection_enabled, stopping_index integer :: error double precision :: minimum_density_parameter, maximum_density_parameter, rho, radius @@ -289,7 +290,7 @@ function evolve_model(tmax) stopping_index = next_index_for_stopping_condition() if (stopping_index > 0) then error = set_stopping_condition_info(stopping_index, DENSITY_LIMIT_DETECTION) - error = set_stopping_condition_particle_index(stopping_index, 0, i) + error = set_stopping_condition_particle_index(stopping_index, 0, int(i, 4)) endif endif ! else @@ -303,7 +304,7 @@ function evolve_model(tmax) function set_state_sph(index_of_the_particle, mass, x, y, z, & vx, vy, vz, u, h_smooth) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz, u, h_smooth integer :: set_state_sph call amuse_set_state_gas(index_of_the_particle, mass, x, y, z, & @@ -321,7 +322,7 @@ function set_eps2(epsilon_squared) function set_state_star(index_of_the_particle, mass, x, y, z, & vx, vy, vz, tform, radius) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz, tform, radius integer :: set_state_star set_state_star=-1 @@ -344,7 +345,7 @@ function get_eps2(epsilon_squared) function get_index_of_next_particle(index_of_the_particle, & index_of_the_next_particle) implicit none - integer :: index_of_the_particle, index_of_the_next_particle + integer(kind=8) :: index_of_the_particle, index_of_the_next_particle integer :: get_index_of_next_particle get_index_of_next_particle=-1 end function @@ -369,17 +370,17 @@ function reset_number_of_new_gas_particles() function new_sph_particle(index_of_the_particle, mass, x, y, z, & vx, vy, vz, u, h_smooth) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz, u, h_smooth integer :: new_sph_particle - call amuse_new_sph_particle(index_of_the_particle, mass, x, y, z, & + call amuse_new_sph_particle(int(index_of_the_particle, 4), mass, x, y, z, & vx, vy, vz, u, h_smooth) new_sph_particle=0 end function function delete_particle(index_of_the_particle) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle integer :: delete_particle call amuse_delete_particle(index_of_the_particle) delete_particle=0 @@ -387,7 +388,7 @@ function delete_particle(index_of_the_particle) function get_potential(index_of_the_particle, potential) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: potential integer :: get_potential get_potential=0 @@ -402,7 +403,7 @@ function synchronize_model() function set_state_sink(index_of_the_particle, mass, x, y, z, & vx, vy, vz, radius, accretion_radius, h_smooth) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius, h_smooth integer :: set_state_sink call amuse_set_state_sink(index_of_the_particle, mass, x, y, z, & @@ -413,7 +414,7 @@ function set_state_sink(index_of_the_particle, mass, x, y, z, & function get_state_sink(index_of_the_particle, mass, x, y, z, & vx, vy, vz, radius, accretion_radius, h_smooth) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius, h_smooth integer :: get_state_sink call amuse_get_state_sink(index_of_the_particle, mass, x, y, z, & @@ -424,7 +425,7 @@ function get_state_sink(index_of_the_particle, mass, x, y, z, & function set_state_dm(index_of_the_particle, mass, x, y, z, & vx, vy, vz) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz integer :: set_state_dm call amuse_set_state_dm(index_of_the_particle, mass, x, y, z, & @@ -435,7 +436,7 @@ function set_state_dm(index_of_the_particle, mass, x, y, z, & function get_state_dm(index_of_the_particle, mass, x, y, z, & vx, vy, vz) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz integer :: get_state_dm call amuse_get_state_dm(index_of_the_particle, mass, x, y, z, & @@ -484,7 +485,7 @@ function get_thermal_energy(thermal_energy) function get_maximum_particle_index(i) implicit none - integer :: i + integer(kind=8) :: i integer :: get_maximum_particle_index call amuse_get_norig(i) get_maximum_particle_index=0 @@ -514,7 +515,7 @@ function get_center_of_mass_velocity(vx, vy, vz) function get_radius(index_of_the_particle, radius) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: radius integer :: get_radius call amuse_get_radius(index_of_the_particle, radius) @@ -523,7 +524,7 @@ function get_radius(index_of_the_particle, radius) function get_sink_accretion_radius(index_of_the_particle, radius) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: radius integer :: get_sink_accretion_radius call amuse_get_sink_accretion_radius(index_of_the_particle, radius) @@ -532,7 +533,7 @@ function get_sink_accretion_radius(index_of_the_particle, radius) function get_sink_temperature(index_of_the_particle, temperature) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: temperature integer :: get_sink_temperature call amuse_get_sink_temperature(index_of_the_particle, temperature) @@ -541,7 +542,7 @@ function get_sink_temperature(index_of_the_particle, temperature) function get_sink_luminosity(index_of_the_particle, luminosity) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: luminosity integer :: get_sink_luminosity call amuse_get_sink_luminosity(index_of_the_particle, luminosity) @@ -551,7 +552,7 @@ function get_sink_luminosity(index_of_the_particle, luminosity) function get_state_star(index_of_the_particle, mass, x, y, z, vx, vy, vz, & tform, radius) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz, tform, radius integer :: get_state_star get_state_star=-1 @@ -574,10 +575,10 @@ integer function set_phantom_option(name, value, match) function new_dm_particle(index_of_the_particle, mass, x, y, z, vx, vy, vz) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz integer :: new_dm_particle - call amuse_new_dm_particle(index_of_the_particle, mass, x, y, z, & + call amuse_new_dm_particle(int(index_of_the_particle, 4), mass, x, y, z, & vx, vy, vz) new_dm_particle=0 end function @@ -585,17 +586,17 @@ function new_dm_particle(index_of_the_particle, mass, x, y, z, vx, vy, vz) function new_sink_particle(index_of_the_particle, mass, x, y, z, vx, vy, vz, & radius, accretion_radius, h_smooth) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius, h_smooth integer :: new_sink_particle - call amuse_new_sink_particle(index_of_the_particle, mass, x, y, z, & + call amuse_new_sink_particle(int(index_of_the_particle, 4), mass, x, y, z, & vx, vy, vz, radius, accretion_radius, h_smooth) new_sink_particle=0 end function function set_radius(index_of_the_particle, radius) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: radius integer :: set_radius call amuse_set_radius(index_of_the_particle, radius) @@ -604,7 +605,7 @@ function set_radius(index_of_the_particle, radius) function set_sink_accretion_radius(index_of_the_particle, radius) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: radius integer :: set_sink_accretion_radius call amuse_set_sink_accretion_radius(index_of_the_particle, radius) @@ -613,7 +614,7 @@ function set_sink_accretion_radius(index_of_the_particle, radius) function set_sink_temperature(index_of_the_particle, temperature) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: temperature integer :: set_sink_temperature call amuse_set_sink_temperature(index_of_the_particle, temperature) @@ -622,7 +623,7 @@ function set_sink_temperature(index_of_the_particle, temperature) function set_sink_luminosity(index_of_the_particle, luminosity) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: luminosity integer :: set_sink_luminosity call amuse_set_sink_luminosity(index_of_the_particle, luminosity) @@ -646,7 +647,7 @@ function get_potential_energy(potential_energy) function get_state_sph(index_of_the_particle, mass, x, y, z, vx, vy, vz, u, & h_smooth) implicit none - integer :: index_of_the_particle, n + integer(kind=8) :: index_of_the_particle, n double precision :: mass, x, y, z, vx, vy, vz, u, h_smooth integer :: get_state_sph get_state_sph = -1 @@ -676,7 +677,7 @@ function get_gravity_at_point(eps, x, y, z, ax, ay, az, npoints) function get_velocity(index_of_the_particle, vx, vy, vz) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: vx, vy, vz integer :: get_velocity call amuse_get_velocity(index_of_the_particle, vx, vy, vz) @@ -685,7 +686,7 @@ function get_velocity(index_of_the_particle, vx, vy, vz) function get_acceleration(index_of_the_particle, ax, ay, az) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: ax, ay, az integer :: get_acceleration call amuse_get_acceleration(index_of_the_particle, ax, ay, az) @@ -695,7 +696,7 @@ function get_acceleration(index_of_the_particle, ax, ay, az) function new_star_particle(index_of_the_particle, mass, x, y, z, vx, vy, & vz, tform, radius) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: mass, x, y, z, vx, vy, vz, tform, radius integer :: new_star_particle new_star_particle=-1 @@ -703,7 +704,7 @@ function new_star_particle(index_of_the_particle, mass, x, y, z, vx, vy, & function get_position(index_of_the_particle, x, y, z) implicit none - integer :: index_of_the_particle, i + integer(kind=8) :: index_of_the_particle, i double precision :: x, y, z integer :: get_position call amuse_get_position(index_of_the_particle, x, y, z) @@ -712,7 +713,7 @@ function get_position(index_of_the_particle, x, y, z) function set_position(index_of_the_particle, x, y, z) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: x, y, z integer :: set_position call amuse_set_position(index_of_the_particle, x, y, z) @@ -728,7 +729,7 @@ function commit_parameters() function set_velocity(index_of_the_particle, vx, vy, vz) implicit none - integer :: index_of_the_particle + integer(kind=8) :: index_of_the_particle double precision :: vx, vy, vz integer :: set_velocity call amuse_set_velocity(index_of_the_particle, vx, vy, vz) diff --git a/src/amuse/community/phantom/interface.py b/src/amuse/community/phantom/interface.py index a4f9486f90..165be07197 100644 --- a/src/amuse/community/phantom/interface.py +++ b/src/amuse/community/phantom/interface.py @@ -9,8 +9,8 @@ ) from amuse.community.interface.gd import ( - GravitationalDynamicsInterface, - GravitationalDynamics, + GravitationalDynamics64Interface, + GravitationalDynamics64, # GravityFieldInterface, GravityFieldCode, ) @@ -24,7 +24,7 @@ class PhantomInterface( CodeInterface, - GravitationalDynamicsInterface, + GravitationalDynamics64Interface, LiteratureReferencesMixIn, StoppingConditionInterface, # GravityFieldInterface, @@ -69,7 +69,7 @@ def new_dm_particle(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.OUT, + 'index_of_the_particle', dtype='int64', direction=function.OUT, ) for x in ['mass', 'x', 'y', 'z', 'vx', 'vy', 'vz']: function.addParameter(x, dtype='float64', direction=function.IN) @@ -80,7 +80,7 @@ def new_dm_particle(): def get_maximum_particle_index(): function = LegacyFunctionSpecification() function.addParameter( - 'norig', dtype='int32', direction=function.OUT) + 'norig', dtype='int64', direction=function.OUT) function.result_type = 'int32' return function @@ -92,7 +92,7 @@ def new_sph_particle(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.OUT, + 'index_of_the_particle', dtype='int64', direction=function.OUT, ) for x in ['mass', 'x', 'y', 'z', 'vx', 'vy', 'vz', 'u']: function.addParameter(x, dtype='float64', direction=function.IN) @@ -107,7 +107,7 @@ def new_sink_particle(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.OUT, + 'index_of_the_particle', dtype='int64', direction=function.OUT, ) for x in ['mass', 'x', 'y', 'z', 'vx', 'vy', 'vz']: function.addParameter(x, dtype='float64', direction=function.IN) @@ -134,7 +134,7 @@ def get_state_dm(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description="""Index of the particle to get the state from. This index must have been returned by an earlier call to :meth:`new_particle`""") @@ -177,7 +177,7 @@ def get_state_sink(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description="""Index of the particle to get the state from. This index must have been returned by an earlier call to :meth:`new_particle`""") @@ -232,7 +232,7 @@ def get_state_sph(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description="""Index of the particle to get the state from. This index must have been returned by an earlier call to :meth:`new_particle`""") @@ -285,7 +285,7 @@ def set_state_sph(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description="""Index of the particle for which the state is to be updated. This index must have been returned by an earlier call to :meth:`new_particle`""") @@ -344,7 +344,7 @@ def set_state_dm(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description="""Index of the particle for which the state is to be updated. This index must have been returned by an earlier call to :meth:`new_particle`""") @@ -394,7 +394,7 @@ def set_state_sink(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description="""Index of the particle for which the state is to be updated. This index must have been returned by an earlier call to :meth:`new_particle`""") @@ -446,7 +446,7 @@ def set_sink_accretion_radius(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -465,7 +465,7 @@ def set_sink_temperature(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -484,7 +484,7 @@ def set_sink_luminosity(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -503,7 +503,7 @@ def set_internal_energy(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -522,7 +522,7 @@ def set_h2ratio(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -541,7 +541,7 @@ def set_hi_abundance(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -560,7 +560,7 @@ def set_proton_abundance(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -579,7 +579,7 @@ def set_electron_abundance(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -598,7 +598,7 @@ def set_co_abundance(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -617,7 +617,7 @@ def set_smoothing_length(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -636,7 +636,7 @@ def get_density(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='' ) function.addParameter( @@ -655,7 +655,7 @@ def get_pressure(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -674,7 +674,7 @@ def get_h2ratio(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -693,7 +693,7 @@ def get_hi_abundance(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -712,7 +712,7 @@ def get_proton_abundance(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -731,7 +731,7 @@ def get_electron_abundance(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -750,7 +750,7 @@ def get_co_abundance(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='', ) function.addParameter( @@ -769,7 +769,7 @@ def get_sink_accretion_radius(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='' ) function.addParameter( @@ -788,7 +788,7 @@ def get_sink_temperature(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='' ) function.addParameter( @@ -807,7 +807,7 @@ def get_sink_luminosity(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='' ) function.addParameter( @@ -826,7 +826,7 @@ def get_internal_energy(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='' ) function.addParameter( @@ -845,7 +845,7 @@ def get_smoothing_length(): function = LegacyFunctionSpecification() function.can_handle_array = True function.addParameter( - 'index_of_the_particle', dtype='int32', direction=function.IN, + 'index_of_the_particle', dtype='int64', direction=function.IN, description='' ) function.addParameter( @@ -1701,7 +1701,7 @@ def reset_number_of_new_gas_particles(): return function -class Phantom(GravitationalDynamics, GravityFieldCode): +class Phantom(GravitationalDynamics64, GravityFieldCode): __interface__ = PhantomInterface def __init__( @@ -1737,7 +1737,7 @@ def __init__( self.stopping_conditions = StoppingConditions(self) - GravitationalDynamics.__init__( + GravitationalDynamics64.__init__( self, PhantomInterface(**options), convert_nbody, @@ -1774,7 +1774,7 @@ def update_gas_particle_set(self): self.gas_particles.remove_particles(deleted_particles) def define_state(self, handler): - GravitationalDynamics.define_state(self, handler) + GravitationalDynamics64.define_state(self, handler) GravityFieldCode.define_state(self, handler) # self.stopping_conditions.define_state(handler) @@ -2137,7 +2137,7 @@ def define_particle_sets(self, handler): self.stopping_conditions.define_particle_set(handler, 'particles') def define_methods(self, handler): - GravitationalDynamics.define_methods(self, handler) + GravitationalDynamics64.define_methods(self, handler) handler.add_method( "new_dm_particle", From d7a0507d4c591c43717b0f7fed6c674437beed84 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Fri, 18 Oct 2024 12:20:49 +0200 Subject: [PATCH 08/10] update phantom version used --- src/amuse/community/phantom/download.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/amuse/community/phantom/download.py b/src/amuse/community/phantom/download.py index 5f5a3d63f7..224b3ad90c 100755 --- a/src/amuse/community/phantom/download.py +++ b/src/amuse/community/phantom/download.py @@ -70,8 +70,7 @@ def new_option_parser(): result = OptionParser() result.add_option( "--version", - # default="ca6907f5f9a95f866b1d7e520cc356ab8cec8dd0", - default="1e7bd070b79f2408e333797f544c3b1daa047c8a", + default="d39f1d6cc8218800187ccd104573ce26142f898e", dest="version", help="version number to download", type="string" From 707f3f8b7a356dbcf5569fb2198ccccc3847c6cc Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Fri, 18 Oct 2024 12:21:41 +0200 Subject: [PATCH 09/10] Download Phantom version --- src/amuse/community/phantom/Makefile | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/amuse/community/phantom/Makefile b/src/amuse/community/phantom/Makefile index d115a1c67c..22c08cf035 100644 --- a/src/amuse/community/phantom/Makefile +++ b/src/amuse/community/phantom/Makefile @@ -32,10 +32,10 @@ $(CODEDIR)/Makefile: download #make -C . download download: - echo "not downloading" - #$(RM) -Rf src - #mkdir src - #$(DOWNLOAD_FROM_WEB) + #echo "not downloading" + $(RM) -Rf src + mkdir src + $(DOWNLOAD_FROM_WEB) clean: $(RM) -f *.so *.o *.pyc worker_code.cc worker_code.h From 43dafdfc54f915a9fd34e50888fb207818dd63eb Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Fri, 18 Oct 2024 12:28:10 +0200 Subject: [PATCH 10/10] plotting updates --- src/amuse/plot/hydro.py | 53 +++++++++++++++++++++++++++++++++----- src/amuse/plot/mapper.py | 55 ++++++++++++++++++++++++++++++++-------- 2 files changed, 91 insertions(+), 17 deletions(-) diff --git a/src/amuse/plot/hydro.py b/src/amuse/plot/hydro.py index 364de588ca..bc07a56046 100755 --- a/src/amuse/plot/hydro.py +++ b/src/amuse/plot/hydro.py @@ -44,16 +44,18 @@ def plot_column_density( plot_axes, maps, unit_col_density=units.MSun * units.pc**-2, - vmin=-1, - vmax=4, + vmin=None, + vmax=None, cmap='viridis', ): + print(f"{vmin} {vmax} VMINMAX") cmap = copy.copy(matplotlib.colormaps[cmap]) cmap.set_bad('k', alpha=1.0) column_density = maps.column_density logscale = numpy.log10(column_density.value_in(unit_col_density)) return plot_axes.imshow( + # column_density.number, logscale, extent=maps.extent, vmin=vmin, @@ -80,6 +82,38 @@ def plot_temperature( return plot_axes.imshow( logscale_temperature_map, extent=maps.extent, + vmin=2, + vmax=6, + cmap=cmap, + origin="lower", + ) + + +def plot_attribute( + plot_axes, + maps, + attribute_name=None, + vmin=0, + vmax=5, + cmap="inferno", + logscale=True, + unit=None, +): + if attribute_name is None: + raise ValueError("need to name attribute") + cmap = copy.copy(matplotlib.colormaps[cmap]) + cmap.set_bad('k', alpha=1.0) + attribute_map = getattr(maps, attribute_name) + if unit is not None: + attribute_map = attribute_map.value_in(unit) + if logscale: + attribute_map = numpy.log10( + attribute_map + ) + + return plot_axes.imshow( + attribute_map, + extent=maps.extent, vmin=vmin, vmax=vmax, cmap=cmap, @@ -216,6 +250,8 @@ def plot_hydro_and_stars( plot="column density", length_unit=units.parsec, colorbar=True, + vmin=None, + vmax=None, **kwargs ): "Plot gas and stars" @@ -236,12 +272,12 @@ def plot_hydro_and_stars( else: fig, ax = new_figure() cax = False - if "stars" in plot: - ax.set_facecolor('black') + # if "stars" in plot: + # ax.set_facecolor('black') gasplot = False if plot == "column density": - gasplot = plot_column_density(ax, maps) + gasplot = plot_column_density(ax, maps, vmin=vmin, vmax=vmax, **kwargs) gasplot_unit = units.MSun * units.pc**-2 elif plot == "temperature": gasplot = plot_temperature(ax, maps) @@ -276,6 +312,9 @@ def plot_hydro_and_stars( extent=maps.extent, origin="lower", ) + else: + gasplot = plot_attribute(ax, maps, attribute_name=plot, vmin=vmin, vmax=vmax, **kwargs) + gasplot_unit = None cmap = copy.copy(matplotlib.colormaps["coolwarm"]) @@ -283,10 +322,10 @@ def plot_hydro_and_stars( (maps.stars is not None) and not maps.stars.is_empty() ): - s = 1 * ( + s = 2 * ( (maps.stars.mass / (7 | units.MSun))**(3.5 / 2) ) - # s = 0.5 + # s = 2 x = getattr(maps.stars, 'x').value_in(length_unit) y = getattr(maps.stars, 'y').value_in(length_unit) z = getattr(maps.stars, 'z').value_in(length_unit) diff --git a/src/amuse/plot/mapper.py b/src/amuse/plot/mapper.py index 6bf21bcc4b..9681812ff5 100644 --- a/src/amuse/plot/mapper.py +++ b/src/amuse/plot/mapper.py @@ -206,6 +206,8 @@ def axes(self, axes): self.__axis_x = axes[0] self.__axis_y = axes[1] self.__axis_z = axes[2] + # TODO: call rotate + # self.__state = "EDIT" @property def phi(self): @@ -361,22 +363,25 @@ def temperature(self): self.__new_gas_mapper() gas = self.__gas - if hasattr(gas, "mu"): - print(f"h2ratio: {gas.mu.mean()}") - meanmwt = gas.mu - elif hasattr(gas, "h2ratio"): - print(f"h2ratio: {gas.h2ratio.mean()}") - meanmwt = gas_mean_molecular_weight(gas.h2ratio) + if hasattr(gas, "temperature"): + temperature = gas.temperature else: - meanmwt = gas_mean_molecular_weight() - temperature = u_to_temperature(gas.u, meanmwt=meanmwt) + if hasattr(gas, "mu"): + print(f"h2ratio: {gas.mu.mean()}") + meanmwt = gas.mu + elif hasattr(gas, "h2ratio"): + print(f"h2ratio: {gas.h2ratio.mean()}") + meanmwt = gas_mean_molecular_weight(gas.h2ratio) + else: + meanmwt = gas_mean_molecular_weight() + temperature = u_to_temperature(gas.u, meanmwt=meanmwt) print(f"temperature range: {min(temperature)} - {max(temperature)}") self.__mapper.particles.weight = temperature.value_in( self.__unit_temperature) # counts = self.counts self.__maps.temperature = numpy.nan_to_num( - self.__mapper.image.pixel_value.transpose(), - # / counts, + self.__mapper.image.pixel_value.transpose() + / self.counts, nan=0, ) | self.__unit_temperature print( @@ -385,6 +390,36 @@ def temperature(self): ) return self.__maps.temperature + @property + def h2o_abundance(self): + "Return a H2O map" + if self.__state != "RUN": + self.__new_gas_mapper() + + self.__mapper.particles.weight = self.__gas.h2o_abundance + self.__maps.h2o = numpy.nan_to_num( + self.__mapper.image.pixel_value.transpose(), + # / self.counts, + nan=0, + ) + + return self.__maps.h2o + + @property + def co_abundance(self): + "Return a CO map" + if self.__state != "RUN": + self.__new_gas_mapper() + + self.__mapper.particles.weight = self.__gas.co_abundance + self.__maps.co = numpy.nan_to_num( + self.__mapper.image.pixel_value.transpose(), + # / self.counts, + nan=0, + ) + + return self.__maps.co + @property def vx(self): "Return a vx map"