diff --git a/.gitignore b/.gitignore index d9b5dcd68b..6009623c25 100644 --- a/.gitignore +++ b/.gitignore @@ -405,7 +405,6 @@ test_results/code-sockets.f90 test_results/code.c test_results/code.cc test_results/code.f90 -gas_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 92% rename from src/amuse/community/interface/gd.py rename to src/amuse/community/interface/gd/gravitational_dynamics.py index 2a0d9059ac..c3bc929379 100644 --- a/src/amuse/community/interface/gd.py +++ b/src/amuse/community/interface/gd/gravitational_dynamics.py @@ -2,9 +2,7 @@ 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 @@ -1132,104 +1130,6 @@ def get_index_of_next_particle(): 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 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 GravitationalDynamicsDocumentation: def __get__(self, instance, owner): @@ -1610,10 +1510,3 @@ def reset(self): def get_total_energy(self): return self.get_potential_energy() + self.get_kinetic_energy() - - -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/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/Makefile b/src/amuse/community/phantom/Makefile index dbcd5da12c..22c08cf035 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,10 +28,11 @@ DOWNLOAD_FROM_WEB = $(PYTHON) ./download.py all: phantom_worker -$(CODEDIR)/Makefile: - make -C . download +$(CODEDIR)/Makefile: download + #make -C . download download: + #echo "not downloading" $(RM) -Rf src mkdir src $(DOWNLOAD_FROM_WEB) @@ -38,17 +40,20 @@ download: 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/download.py b/src/amuse/community/phantom/download.py index 959efcc131..224b3ad90c 100755 --- a/src/amuse/community/phantom/download.py +++ b/src/amuse/community/phantom/download.py @@ -70,7 +70,7 @@ def new_option_parser(): result = OptionParser() result.add_option( "--version", - default="ca6907f5f9a95f866b1d7e520cc356ab8cec8dd0", + default="d39f1d6cc8218800187ccd104573ce26142f898e", dest="version", help="version number to download", type="string" diff --git a/src/amuse/community/phantom/interface.f90 b/src/amuse/community/phantom/interface.f90 index cec2cd2419..9767f234e7 100644 --- a/src/amuse/community/phantom/interface.f90 +++ b/src/amuse/community/phantom/interface.f90 @@ -1,11 +1,19 @@ -!MODULE PhantomInterface -! -!CONTAINS - -function initialize_code() +module AmuseInterface + use AmusePhantom use StoppingConditions + +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 + +integer function initialize_code() result(ret) implicit none - integer :: initialize_code + !integer :: initialize_code integer :: error double precision :: polyk call amuse_initialize_code() @@ -15,7 +23,7 @@ function initialize_code() !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() @@ -25,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 @@ -42,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) @@ -51,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) @@ -60,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) @@ -69,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) @@ -78,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) @@ -87,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) @@ -96,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) @@ -105,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) @@ -114,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) @@ -123,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) @@ -132,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) @@ -141,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) @@ -150,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) @@ -159,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) @@ -168,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) @@ -177,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) @@ -186,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) @@ -195,16 +211,25 @@ 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) get_internal_energy=0 end function +function get_particle_type(index_of_the_particle, itype) + implicit none + integer(kind=8) :: 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 + integer(kind=8) :: index_of_the_particle integer :: get_index_of_first_particle get_index_of_first_particle=0 end function @@ -237,10 +262,12 @@ 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 + double precision :: time error = reset_stopping_conditions() error = is_stopping_condition_enabled(& @@ -248,11 +275,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 (& @@ -262,9 +290,11 @@ 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 + ! nmax += 1 endif enddo endif @@ -274,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, & @@ -292,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 @@ -315,25 +345,42 @@ 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 +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 - 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 @@ -341,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 @@ -354,31 +401,31 @@ 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 + 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, & - vx, vy, vz, radius, h_smooth) + vx, vy, vz, radius, accretion_radius) 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 + 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, & - vx, vy, vz, radius, h_smooth) + vx, vy, vz, radius, accretion_radius) get_state_sink=0 end function 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, & @@ -389,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, & @@ -436,6 +483,14 @@ function get_thermal_energy(thermal_energy) get_thermal_energy=0 end function +function get_maximum_particle_index(i) + implicit none + integer(kind=8) :: 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 @@ -460,17 +515,44 @@ 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) get_radius=0 end function +function get_sink_accretion_radius(index_of_the_particle, radius) + implicit none + 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) + get_sink_accretion_radius=0 +end function + +function get_sink_temperature(index_of_the_particle, temperature) + implicit none + integer(kind=8) :: 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(kind=8) :: 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 - 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 @@ -483,36 +565,71 @@ 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 - double precision :: mass, x, y, z, vx, vy, vz, radius + 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, & - vx, vy, vz, radius) + call amuse_new_dm_particle(int(index_of_the_particle, 4), mass, x, y, z, & + vx, vy, vz) new_dm_particle=0 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 + 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, & - vx, vy, vz, radius, h_smooth) + 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) set_radius=0 end function +function set_sink_accretion_radius(index_of_the_particle, radius) + implicit none + 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) + set_sink_accretion_radius=0 +end function + +function set_sink_temperature(index_of_the_particle, temperature) + implicit none + integer(kind=8) :: 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(kind=8) :: 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 @@ -527,14 +644,27 @@ 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 + integer(kind=8) :: 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) + 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, phantom_index) + !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) @@ -547,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) @@ -556,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) @@ -566,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 @@ -574,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) @@ -583,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) @@ -593,12 +723,13 @@ 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 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) @@ -1125,4 +1256,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..165be07197 100644 --- a/src/amuse/community/phantom/interface.py +++ b/src/amuse/community/phantom/interface.py @@ -5,11 +5,12 @@ LegacyFunctionSpecification, legacy_function, LiteratureReferencesMixIn, + CodeWithDataDirectories, ) from amuse.community.interface.gd import ( - GravitationalDynamicsInterface, - GravitationalDynamics, + GravitationalDynamics64Interface, + GravitationalDynamics64, # GravityFieldInterface, GravityFieldCode, ) @@ -22,11 +23,12 @@ class PhantomInterface( - CodeInterface, - LiteratureReferencesMixIn, - GravitationalDynamicsInterface, - StoppingConditionInterface, - # SinglePointGravityFieldInterface, + CodeInterface, + GravitationalDynamics64Interface, + 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__( @@ -42,18 +45,45 @@ 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() 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) function.result_type = 'int32' return function + @legacy_function + def get_maximum_particle_index(): + function = LegacyFunctionSpecification() + function.addParameter( + 'norig', dtype='int64', 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) @@ -62,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) @@ -77,12 +107,15 @@ 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) 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( @@ -101,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`""") @@ -144,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`""") @@ -171,6 +204,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, @@ -196,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`""") @@ -232,7 +268,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 @@ -245,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`""") @@ -304,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`""") @@ -354,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`""") @@ -381,6 +421,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, @@ -398,12 +441,69 @@ 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='int64', 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() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int64', 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='int64', 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() 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( @@ -422,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( @@ -441,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( @@ -460,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( @@ -479,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( @@ -498,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( @@ -517,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( @@ -536,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( @@ -555,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( @@ -574,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( @@ -593,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( @@ -612,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( @@ -631,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( @@ -650,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( @@ -664,12 +764,69 @@ 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='int64', 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() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int64', 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='int64', 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() 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( @@ -688,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( @@ -1521,8 +1678,30 @@ 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 -class Phantom(GravitationalDynamics, GravityFieldCode): + @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(GravitationalDynamics64, GravityFieldCode): __interface__ = PhantomInterface def __init__( @@ -1534,6 +1713,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 @@ -1541,7 +1721,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 * phantom_au phantom_time = (phantom_length**3 / (phantom_gg*phantom_mass))**0.5 unit_converter = ConvertBetweenGenericAndSiUnits( @@ -1556,7 +1737,7 @@ def __init__( self.stopping_conditions = StoppingConditions(self) - GravitationalDynamics.__init__( + GravitationalDynamics64.__init__( self, PhantomInterface(**options), convert_nbody, @@ -1576,14 +1757,32 @@ 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) + GravitationalDynamics64.define_state(self, handler) GravityFieldCode.define_state(self, handler) # self.stopping_conditions.define_state(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') @@ -1926,13 +2125,19 @@ 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') + handler.add_setter('sink_particles', 'set_sink_luminosity') handler.add_getter('sink_particles', 'get_smoothing_length') handler.add_setter('sink_particles', 'set_smoothing_length') 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", @@ -1982,6 +2187,7 @@ def define_methods(self, handler): generic_unit_system.speed, generic_unit_system.length, generic_unit_system.length, + generic_unit_system.length, ), ( handler.INDEX, @@ -2077,6 +2283,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, ) ) @@ -2094,8 +2301,75 @@ 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", + ( + 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, + generic_unit_system.energy / generic_unit_system.time, + ), + ( + handler.ERROR_CODE, + ) + ) + + handler.add_method( + "get_sink_luminosity", + ( + handler.INDEX, ), ( + generic_unit_system.energy / generic_unit_system.time, handler.ERROR_CODE, ) ) 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"