diff --git a/src/amuse/community/mizuki/Makefile b/src/amuse/community/mizuki/Makefile new file mode 100644 index 0000000000..01ceb4fc21 --- /dev/null +++ b/src/amuse/community/mizuki/Makefile @@ -0,0 +1,75 @@ +# standard amuse configuration include +# config.mk will be made after ./configure has run +ifeq ($(origin AMUSE_DIR), undefined) +AMUSE_DIR := $(shell amusifier --get-amuse-dir) +endif +-include $(AMUSE_DIR)/config.mk + +MPICXX ?= mpicxx + +use_phantom_grape_x86 = no # speedup from ~1M particles or so? + +FDPS_PATH = src/FDPS/src +SUI_PATH = src/mizuki + +INCLUDE += -I $(FDPS_PATH) -I $(SUI_PATH) + +MPIFLAGS += -DPARTICLE_SIMULATOR_MPI_PARALLEL + +CXXFLAGS = $(INCLUDE) -fPIC -O3 -Wall -std=c++20 -fopenmp + + +CXXFLAGS += -DPARTICLE_SIMULATOR_THREAD_PARALLEL +CXXFLAGS += -DENABLE_HYDRO_INTERACT +# CXXFLAGS += -DPARTICLE_SIMULATOR_MPI_PARALLEL +CXXFLAGS += -DENABLE_VARIABLE_SMOOTHING_LENGTH +# CXXFLAGS += -DUSE_ENTROPY +# CXXFLAGS += -DUSE_BALSARA_SWITCH +# CXXFLAGS += -DUSE_PRESCR_OF_THOMAS_COUCHMAN_1992 +# CXXFLAGS += -DISOTHERMAL_EOS + +ifeq ($(use_phantom_grape_x86),yes) + PG_ROOT = $(FDPS_PATH)/phantom_grape_x86/G5/newton/libpg5 + CXXFLAGS += -I $(PG_ROOT) + CXXFLAGS += -DENABLE_PHANTOM_GRAPE_X86 + CXXFLAGS += -L$(PG_ROOT) -lpg5 + PG_BUILD = make -C $(PG_ROOT) distclean libpg5.a + PG_CLEAN = make -C $(PG_ROOT) distclean +else + PG_BUILD = + PG_CLEAN = +endif + +LDFLAGS += -lm $(MUSE_LD_FLAGS) + +OBJS = interface.o + + +all: mizuki_worker $(CODELIB) + +$(CODELIB): + $(PG_BUILD) + #amusi $(MPICXX) $(CXXFLAGS) + #$(PG_BUILD) + #$(AR) $@ + +clean: + $(RM) -rf __pycache__ + $(RM) -f *.so *.o *.pyc worker_code.cc worker_code.h + $(RM) *~ mizuki_worker worker_code.cc + make -C src clean + +distclean: clean + make -C src distclean + +worker_code.cc: interface.py + $(CODE_GENERATOR) --type=c interface.py MizukiInterface -o $@ + +worker_code.h: interface.py + $(CODE_GENERATOR) --type=H -i amuse.support.codes.stopping_conditions.StoppingConditionInterface interface.py MizukiInterface -o $@ + +mizuki_worker: worker_code.cc worker_code.h $(CODELIB) $(OBJS) interface.hpp + $(MPICXX) $(CXXFLAGS) $(SC_FLAGS) $(LDFLAGS) $< $(OBJS) -o $@ + +interface.o: interface.cc + $(MPICXX) $(CXXFLAGS) $(SC_FLAGS) -c -o $@ $< diff --git a/src/amuse/community/mizuki/__init__.py b/src/amuse/community/mizuki/__init__.py new file mode 100644 index 0000000000..a6ecf098d2 --- /dev/null +++ b/src/amuse/community/mizuki/__init__.py @@ -0,0 +1 @@ +from .interface import Mizuki diff --git a/src/amuse/community/mizuki/download.py b/src/amuse/community/mizuki/download.py new file mode 100644 index 0000000000..f27a057f53 --- /dev/null +++ b/src/amuse/community/mizuki/download.py @@ -0,0 +1,108 @@ +#!/usr/bin/env python + +import subprocess +import os +import urllib.request +import urllib.parse +import urllib.error +from shutil import which +from optparse import OptionParser + + +class GetCodeFromHttp: + filename_template = "{version}.tar.gz" + name = ["FDPS"] + url_template = [ + "https://github.com/FDPS/FDPS/archive/{version}.tar.gz", + ] + version = [ + "", + ] + + def directory(self): + return os.path.abspath(os.path.dirname(__file__)) + + def src_directory(self): + return os.path.join(self.directory(), 'src') + + def unpack_downloaded_file(self, filename, name, version): + print("unpacking", filename) + arguments = ['tar', '-xf'] + arguments.append(filename) + subprocess.call( + arguments, + cwd=os.path.join(self.src_directory()) + ) + subprocess.call( + [ + 'mv', '{name}-{version}'.format(name=name, version=version), + name + ], + cwd=os.path.join(self.src_directory()) + ) + print("done") + + def start(self): + # if os.path.exists('src'): + # counter = 0 + # while os.path.exists('src.{0}'.format(counter)): + # counter += 1 + # if counter > 100: + # print("too many backup directories") + # break + # os.rename('src', 'src.{0}'.format(counter)) + + # os.mkdir('src') + + for i, url_template in enumerate(self.url_template): + url = url_template.format(version=self.version[i]) + filename = self.filename_template.format(version=self.version[i]) + filepath = os.path.join(self.src_directory(), filename) + print( + "downloading version", self.version[i], + "from", url, "to", filename + ) + if which('wget') is not None: + arguments = ['wget', url] + subprocess.call( + arguments, + cwd=os.path.join(self.src_directory()) + ) + elif which('curl') is not None: + arguments = ['curl', '-L', '-O', url] + subprocess.call( + arguments, + cwd=os.path.join(self.src_directory()) + ) + else: + urllib.request.urlretrieve(url, filepath) + print("downloading finished") + self.unpack_downloaded_file( + filename, self.name[i], self.version[i] + ) + + +def main(fdps_version=''): + version = [ + fdps_version, + ] + instance = GetCodeFromHttp() + instance.version = version + instance.start() + + +def new_option_parser(): + result = OptionParser() + result.add_option( + "--fdps-version", + default='57c73ed2213bfbd9dbd05700fd60c8922be22103', + dest="fdps_version", + help="FDPS commit hash to download", + type="string" + ) + return result + + +if __name__ == "__main__": + options, arguments = new_option_parser().parse_args() + main(**options.__dict__) diff --git a/src/amuse/community/mizuki/interface.cc b/src/amuse/community/mizuki/interface.cc new file mode 100644 index 0000000000..60ab7525d1 --- /dev/null +++ b/src/amuse/community/mizuki/interface.cc @@ -0,0 +1,407 @@ +#include "worker_code.h" +#include "stopcond.h" +#include "src/mizuki/user_defined.hpp" +#include "src/mizuki/mizuki.hpp" + +static Mizuki* mizuki = new Mizuki; + +int initialize_code(){ + mizuki->initialize(); + return 0; +} + +int commit_particles(){ + mizuki->commit_particles(); + return 0; +} + +int new_sph_particle(int * index_of_the_particle, double mass, double x, + double y, double z, double vx, double vy, double vz, double eng){ + *index_of_the_particle = mizuki->add_sph_particle( + mass, x, y, z, vx, vy, vz, eng + ); + return 0; +} + +int new_particle(int * index_of_the_particle, double mass, double x, + double y, double z, double vx, double vy, double vz, double radius){ + return 0; +} + +int get_mass(int index_of_the_particle, double * mass){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *mass = p->mass; + return 0; +} + +int set_mass(int index_of_the_particle, double mass){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + p->mass = mass; + return 0; +} + +int get_position(int index_of_the_particle, double * x, double * y, + double * z){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *x = p->pos.x; + *y = p->pos.y; + *z = p->pos.z; + return 0; +} + +int set_position(int index_of_the_particle, double x, double y, + double z){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + p->pos.x = x; + p->pos.y = y; + p->pos.z = z; + return 0; +} + +int get_velocity(int index_of_the_particle, double * vx, double * vy, + double * vz){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *vx = p->vel.x; + *vy = p->vel.y; + *vz = p->vel.z; + return 0; +} + +int set_velocity(int index_of_the_particle, double vx, double vy, + double vz){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + p->vel.x = vx; + p->vel.y = vy; + p->vel.z = vz; + return 0; +} + +int get_internal_energy(int index_of_the_particle, double * eng){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *eng = p->eng; + return 0; +} + +int set_internal_energy(int index_of_the_particle, double eng){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + p->eng = eng; + return 0; +} + +int get_smoothing_length(int index_of_the_particle, double * smth){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *smth = p->smth/2; + return 0; +} + +int set_smoothing_length(int index_of_the_particle, double smth){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + p->smth = smth*2; + return 0; +} + +int get_radius(int index_of_the_particle, double * radius){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *radius = p->smth/2; + return 0; +} + +int set_radius(int index_of_the_particle, double radius){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + p->smth = radius*2; + return 0; +} + +int get_density(int index_of_the_particle, double * dens){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *dens = p->dens; + return 0; +} + +int get_pressure(int index_of_the_particle, double * pres){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *pres = p->pres; + return 0; +} + +int get_state_sph(int index_of_the_particle, double * mass, double * x, + double * y, double * z, double * vx, double * vy, double * vz, + double * eng, double * smth){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *mass = p->mass; + *x = p->pos.x; + *y = p->pos.y; + *z = p->pos.z; + *vx = p->vel.x; + *vy = p->vel.y; + *vz = p->vel.z; + *eng = p->eng; + *smth = p->smth/2; + return 0; +} + +int set_state_sph(int index_of_the_particle, double mass, double x, + double y, double z, double vx, double vy, double vz, + double eng){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + p->mass = mass; + p->pos.x = x; + p->pos.y = y; + p->pos.z = z; + p->vel.x = vx; + p->vel.y = vy; + p->vel.z = vz; + p->eng = eng; + return 0; +} + +int get_state(int index_of_the_particle, double * mass, double * x, + double * y, double * z, double * vx, double * vy, double * vz, + double * radius){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *mass = p->mass; + *x = p->pos.x; + *y = p->pos.y; + *z = p->pos.z; + *vx = p->vel.x; + *vy = p->vel.y; + *vz = p->vel.z; + *radius = 0; + return 0; +} + +int set_stopping_condition_timeout_parameter(double value){ + return 0; +} + +int get_time(double * time){ + *time = mizuki->time; + return 0; +} + +int get_stopping_condition_maximum_density_parameter(double * value){ + return 0; +} + +int set_stopping_condition_out_of_box_use_center_of_mass_parameter( + _Bool value){ + return 0; +} + +int get_index_of_first_particle(int * index_of_the_particle){ + return 0; +} + +int get_stopping_condition_out_of_box_use_center_of_mass_parameter( + _Bool * value){ + return 0; +} + +int enable_stopping_condition(int type){ + return 0; +} + +int get_total_radius(double * radius){ + return 0; +} + +int get_stopping_condition_maximum_internal_energy_parameter( + double * value){ + return 0; +} + +int get_potential_at_point(double eps, double x, double y, double z, + double * phi, int npoints){ + return 0; +} + +int get_number_of_stopping_conditions_set(int * result){ + return 0; +} + +int is_stopping_condition_set(int type, int * result){ + return 0; +} + +int get_total_mass(double * mass){ + return 0; +} + +int evolve_model(double time){ + mizuki->evolve_model(time); + return 0; +} + +int set_stopping_condition_out_of_box_parameter(double value){ + return 0; +} + +int set_eps2(double epsilon_squared){ + return 0; +} + +int set_stopping_condition_number_of_steps_parameter(int value){ + return 0; +} + +int get_stopping_condition_timeout_parameter(double * value){ + return 0; +} + +int get_begin_time(double * time){ + return 0; +} + +int get_eps2(double * epsilon_squared){ + *epsilon_squared = mizuki->epsilon_gravity * mizuki->epsilon_gravity; + return 0; +} + +int get_stopping_condition_minimum_internal_energy_parameter( + double * value){ + return 0; +} + +int get_index_of_next_particle(int index_of_the_particle, + int * index_of_the_next_particle){ + return 0; +} + +int delete_particle(int index_of_the_particle){ + return 0; +} + +int is_stopping_condition_enabled(int type, int * result){ + return 0; +} + +int get_potential(int index_of_the_particle, double * potential){ + return 0; +} + +int synchronize_model(){ + return 0; +} + +int set_state(int index_of_the_particle, double mass, double x, double y, + double z, double vx, double vy, double vz, double radius){ + return 0; +} + +int get_stopping_condition_minimum_density_parameter(double * value){ + return 0; +} + +int get_time_step(double * time_step){ + *time_step = mizuki->dt_max; + return 0; +} + +int set_time_step(double time_step){ + mizuki->dt_max = time_step; + mizuki->getTimeStep(); + return 0; +} + +int recommit_particles(){ + return 0; +} + +int get_kinetic_energy(double * kinetic_energy){ + return 0; +} + +int get_number_of_particles(int * number_of_particles){ + return 0; +} + +int get_stopping_condition_number_of_steps_parameter(int * value){ + return 0; +} + +int disable_stopping_condition(int type){ + return 0; +} + +int set_acceleration(int index_of_the_particle, double ax, double ay, + double az){ + return 0; +} + +int get_center_of_mass_position(double * x, double * y, double * z){ + return 0; +} + +int get_center_of_mass_velocity(double * vx, double * vy, double * vz){ + return 0; +} + +int set_stopping_condition_minimum_internal_energy_parameter( + double value){ + return 0; +} + +int set_begin_time(double time){ + return 0; +} + +int set_stopping_condition_minimum_density_parameter(double value){ + return 0; +} + +int has_stopping_condition(int type, int * result){ + return 0; +} + +int cleanup_code(){ + return 0; +} + +int set_stopping_condition_maximum_density_parameter(double value){ + return 0; +} + +int recommit_parameters(){ + return 0; +} + +int get_potential_energy(double * potential_energy){ + return 0; +} + +int get_gravity_at_point(double eps, double x, double y, double z, + double * ax, double * ay, double * az, int npoints){ + return 0; +} + +int get_stopping_condition_out_of_box_parameter(double * value){ + return 0; +} + +int set_stopping_condition_maximum_internal_energy_parameter( + double value){ + return 0; +} + +int get_stopping_condition_info(int index, int * type, + int * number_of_particles){ + return 0; +} + +int get_acceleration(int index_of_the_particle, double * ax, double * ay, + double * az){ + FP_sph* p = &(mizuki->psys_sph[index_of_the_particle]); + *ax = p->acc_hydro.x + p->acc_grav.x; + *ay = p->acc_hydro.y + p->acc_grav.y; + *az = p->acc_hydro.z + p->acc_grav.z; + return 0; +} + +int commit_parameters(){ + return 0; +} + +int get_stopping_condition_particle_index(int index, + int index_of_the_column, int * index_of_particle){ + return 0; +} diff --git a/src/amuse/community/mizuki/interface.py b/src/amuse/community/mizuki/interface.py new file mode 100644 index 0000000000..fb28da7b12 --- /dev/null +++ b/src/amuse/community/mizuki/interface.py @@ -0,0 +1,869 @@ +from amuse.community import ( + CodeInterface, + LegacyFunctionSpecification, + legacy_function, + LiteratureReferencesMixIn, +) + +from amuse.community.interface.gd import ( + GravitationalDynamicsInterface, + GravitationalDynamics, + # GravityFieldInterface, + # GravityFieldCode, +) +from amuse.community.interface.stopping_conditions import ( + StoppingConditionInterface, + StoppingConditions, +) +from amuse.units import nbody_system + + +class MizukiInterface( + CodeInterface, + LiteratureReferencesMixIn, + GravitationalDynamicsInterface, + StoppingConditionInterface, + # GravityFieldInterface, +): + """ + Mizuki: based on c++ nbody+sph example of FDPS + + References: + .. [#] Iwasawa et al. (FDPS) + """ + + include_headers = ['worker_code.h', 'stopcond.h'] + + def __init__(self, **keyword_arguments): + CodeInterface.__init__( + self, + name_of_the_worker="mizuki_worker", + **keyword_arguments + ) + LiteratureReferencesMixIn.__init__(self) + + @legacy_function + def new_sph_particle(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.OUT, + ) + for x in ['mass', 'x', 'y', 'z', 'vx', 'vy', 'vz', 'u']: + function.addParameter(x, dtype='float64', direction=function.IN) + # NOTE: Do not set h_smooth manually - it will be calculated by Mizuki! + # function.addParameter( + # 'h_smooth', dtype='float64', direction=function.IN, default=0.01, + # ) + function.result_type = 'int32' + return function + + + def get_state(self, index_of_the_particle): + return self.get_state_sph(index_of_the_particle) + + @legacy_function + def get_state_sph(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN, + description="""Index of the particle to get the state from. This + index must have been returned by an earlier call to + :meth:`new_particle`""") + function.addParameter( + 'mass', dtype='float64', direction=function.OUT, + description="The current mass of the particle") + 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( + 'u', dtype='float64', direction=function.OUT, + description="The current specific energy of the particle") + function.addParameter( + 'h_smooth', dtype='float64', direction=function.OUT, + description="The current smoothing length of the particle") + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + particle was found + -1 - ERROR + particle could not be found""" + return function + + @legacy_function + def set_state_sph(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN, + description="Identifier of the 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( + 'u', dtype='float64', direction=function.IN, + description="The new specific energy of the particle") + # function.addParameter( + # 'h_smooth', dtype='float64', direction=function.IN, + # description="The new smoothing length of the particle") + 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 + """ + 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.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN, + description=( + "Index of the particle to get the state from." + " This index must have been returned by an earlier call to" + " :meth:`new_particle`" + ) + ) + function.addParameter( + 'mass', dtype='float64', direction=function.OUT, + description="The current mass of the particle" + ) + 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 get_position(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN) + function.addParameter( + 'x', dtype='float64', direction=function.OUT, + unit=nbody_system.length, + ) + function.addParameter( + 'y', dtype='float64', direction=function.OUT, + unit=nbody_system.length, + ) + function.addParameter( + 'z', dtype='float64', direction=function.OUT, + unit=nbody_system.length, + ) + function.result_type = 'int32' + return function + + @legacy_function + def get_velocity(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN) + function.addParameter( + 'vx', dtype='float64', direction=function.OUT, + unit=nbody_system.speed, + ) + function.addParameter( + 'vy', dtype='float64', direction=function.OUT, + unit=nbody_system.speed, + ) + function.addParameter( + 'vz', dtype='float64', direction=function.OUT, + unit=nbody_system.speed, + ) + function.result_type = 'int32' + return function + + @legacy_function + def get_density(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN) + function.addParameter( + 'density', dtype='float64', direction=function.OUT, + description="The density of the particle", + unit=nbody_system.density, + ) + function.result_type = 'int32' + return function + + @legacy_function + def get_pressure(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN) + function.addParameter( + 'pressure', dtype='float64', direction=function.OUT, + description="The pressure of the particle", + unit=nbody_system.pressure, + ) + function.result_type = 'int32' + return function + + @legacy_function + def get_smoothing_length(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN) + function.addParameter( + 'h_smooth', dtype='float64', direction=function.OUT, + description="Smoothing length", + unit=nbody_system.length, + ) + function.result_type = 'int32' + return function + + # @legacy_function + # def set_smoothing_length(): + # function = LegacyFunctionSpecification() + # function.can_handle_array = True + # function.addParameter( + # 'index_of_the_particle', dtype='int32', direction=function.IN) + # function.addParameter( + # 'h_smooth', dtype='float64', direction=function.IN, + # description="Smoothing length", + # unit=nbody_system.length, + # ) + # function.result_type = 'int32' + # 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) + function.addParameter( + 'u', dtype='float64', direction=function.OUT, + description='Specific energy', + ) + function.result_type = 'int32' + 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, + description='', + ) + function.addParameter( + 'u', dtype='float64', direction=function.IN, + description='', + ) + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + -1 - ERROR + """ + return function + + @legacy_function + def get_radius(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN) + function.addParameter( + 'radius', dtype='float64', direction=function.OUT, + description="Smoothing length", + unit=nbody_system.length, + ) + function.result_type = 'int32' + return function + + @legacy_function + def set_radius(): + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + 'index_of_the_particle', dtype='int32', direction=function.IN) + function.addParameter( + 'radius', dtype='float64', direction=function.IN, + description="Smoothing length", + unit=nbody_system.length, + ) + function.result_type = 'int32' + return function + + @legacy_function + def get_time_step(): + function = LegacyFunctionSpecification() + function.addParameter( + 'time_step', dtype='float64', direction=function.OUT, + unit=nbody_system.time, + ) + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + -1 - ERROR + """ + return function + + @legacy_function + def set_time_step(): + function = LegacyFunctionSpecification() + function.addParameter( + 'time_step', dtype='float64', direction=function.IN, + unit=nbody_system.time, + ) + function.result_type = 'int32' + function.result_doc = """ + 0 - OK + -1 - ERROR + """ + return function + +class Mizuki(GravitationalDynamics,): + __interface__ = MizukiInterface + + def __init__( + self, + convert_nbody=None, + **options + ): + self.stopping_conditions = StoppingConditions(self) + + GravitationalDynamics.__init__( + self, + MizukiInterface(**options), + convert_nbody, + **options + ) + + def define_state(self, handler): + GravitationalDynamics.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('EDIT', 'new_sph_particle') + handler.add_method('UPDATE', 'new_sph_particle') + handler.add_transition('RUN', 'UPDATE', 'new_sph_particle', False) + handler.add_method('RUN', 'get_velocity') + handler.add_method('RUN', 'get_acceleration') + handler.add_method('RUN', 'get_internal_energy') + # handler.add_method('RUN', 'get_dinternal_energy_dt') + handler.add_method('RUN', 'get_smoothing_length') + handler.add_method('RUN', 'get_density') + handler.add_method('RUN', 'get_pressure') + handler.add_method('RUN', 'get_state_sph') + + self.stopping_conditions.define_state(handler) + + def define_parameters(self, handler): + handler.add_method_parameter( + "get_time_step", + "set_time_step", + "time_step", + "Maximum internal time step", + default_value=0.01 | nbody_system.time + ) + # handler.add_method_parameter( + # "get_eps2", + # "set_eps2", + # "epsilon_squared", + # "smoothing parameter for gravity calculations", + # default_value=0.0 | nbody_system.length * nbody_system.length + # ) + + # handler.add_method_parameter( + # "get_dtime", + # "set_dtime", + # "timestep", + # "timestep for system", + # default_value=1.0/8 | nbody_system.time + # ) + + # handler.add_boolean_parameter( + # "get_use_hydro", + # "set_use_hydro", + # "use_hydro_flag", + # "Hydrodynamics flag. True means: SPH hydro included," + # " False means: gravity only.", + # True + # ) + pass + + def define_particle_sets(self, handler): + handler.define_super_set( + 'particles', + ['gas_particles', ], + index_to_default_set=0 + ) + + # handler.define_set('dm_particles', 'index_of_the_particle') + # handler.set_new('dm_particles', 'new_dm_particle') + # handler.set_delete('dm_particles', 'delete_particle') + # handler.add_setter('dm_particles', 'set_state') + # handler.add_getter('dm_particles', 'get_state') + # handler.add_setter('dm_particles', 'set_mass') + # handler.add_getter('dm_particles', 'get_mass', names = ('mass',)) + # handler.add_setter('dm_particles', 'set_position') + # handler.add_getter('dm_particles', 'get_position') + # handler.add_setter('dm_particles', 'set_radius') + # handler.add_getter('dm_particles', 'get_radius') + # handler.add_setter('dm_particles', 'set_velocity') + # handler.add_getter('dm_particles', 'get_velocity') + + handler.define_set('gas_particles', 'index_of_the_particle') + handler.set_new('gas_particles', 'new_sph_particle') + handler.set_delete('gas_particles', 'delete_particle') + handler.add_getter('gas_particles', 'get_state_sph') + handler.add_setter('gas_particles', 'set_state_sph') + handler.add_getter('gas_particles', 'get_mass') + handler.add_setter('gas_particles', 'set_mass') + # handler.add_getter('gas_particles', 'get_radius') + # handler.add_setter('gas_particles', 'set_radius') + handler.add_getter('gas_particles', 'get_position') + handler.add_setter('gas_particles', 'set_position') + handler.add_getter('gas_particles', 'get_velocity') + handler.add_setter('gas_particles', 'set_velocity') + handler.add_getter('gas_particles', 'get_acceleration') + handler.add_getter('gas_particles', 'get_internal_energy') + handler.add_setter('gas_particles', 'set_internal_energy') + # handler.add_getter('gas_particles', 'get_dinternal_energy_dt') + handler.add_getter('gas_particles', 'get_smoothing_length') + # handler.add_setter('gas_particles', 'set_smoothing_length') + handler.add_getter('gas_particles', 'get_density', names=('rho',)) + handler.add_getter('gas_particles', 'get_density', names=('density',)) + handler.add_getter('gas_particles', 'get_pressure') + + def define_methods(self, handler): + GravitationalDynamics.define_methods(self, handler) + + # handler.add_method( + # "new_dm_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( + "new_sph_particle", + ( + nbody_system.mass, + nbody_system.length, + nbody_system.length, + nbody_system.length, + nbody_system.speed, + nbody_system.speed, + nbody_system.speed, + nbody_system.specific_energy, + # nbody_system.length, + ), + ( + handler.INDEX, + handler.ERROR_CODE, + ) + ) + + handler.add_method( + "get_state_sph", + ( + handler.INDEX, + ), + ( + nbody_system.mass, + nbody_system.length, + nbody_system.length, + nbody_system.length, + nbody_system.speed, + nbody_system.speed, + nbody_system.speed, + nbody_system.specific_energy, + nbody_system.length, + handler.ERROR_CODE + ) + ) + + handler.add_method( + "set_state_sph", + ( + handler.INDEX, + nbody_system.mass, + nbody_system.length, + nbody_system.length, + nbody_system.length, + nbody_system.speed, + nbody_system.speed, + nbody_system.speed, + nbody_system.specific_energy, + ), + ( + handler.ERROR_CODE, + ) + ) + + handler.add_method( + "get_density", + ( + handler.INDEX, + ), + ( + nbody_system.density, + handler.ERROR_CODE, + ) + ) + handler.add_method( + "get_smoothing_length", + ( + handler.INDEX, + ), + ( + nbody_system.length, + handler.ERROR_CODE, + ) + ) + + # handler.add_method( + # "set_smoothing_length", + # ( + # handler.INDEX, + # nbody_system.length, + # ), + # ( + # handler.ERROR_CODE, + # ) + # ) + + handler.add_method( + "get_pressure", + ( + handler.INDEX, + ), + ( + nbody_system.pressure, + handler.ERROR_CODE, + ) + ) + + # handler.add_method( + # "set_velocity", + # ( + # handler.INDEX, + # nbody_system.speed, + # nbody_system.speed, + # nbody_system.speed, + # ), + # ( + # handler.ERROR_CODE + # ) + # ) + + # handler.add_method( + # "get_velocity", + # ( + # handler.INDEX, + # ), + # ( + # nbody_system.speed, + # nbody_system.speed, + # nbody_system.speed, + # handler.ERROR_CODE + # ) + # ) + + # handler.add_method( + # "get_acceleration", + # ( + # handler.INDEX, + # ), + # ( + # nbody_system.speed / nbody_system.time, + # nbody_system.speed / nbody_system.time, + # nbody_system.speed / nbody_system.time, + # handler.ERROR_CODE + # ) + # ) + + handler.add_method( + "get_internal_energy", + ( + handler.INDEX, + ), + ( + nbody_system.specific_energy, + handler.ERROR_CODE + ) + ) + + handler.add_method( + "set_internal_energy", + ( + handler.INDEX, + nbody_system.specific_energy, + ), + ( + handler.ERROR_CODE, + ) + ) + # handler.add_method( + # "get_dinternal_energy_dt", + # ( + # handler.INDEX, + # ), + # ( + # nbody_system.specific_energy/nbody_system.time, + # handler.ERROR_CODE + # ) + # ) + + # handler.add_method( + # "new_star_particle", + # ( + # nbody_system.mass, + # nbody_system.length, + # nbody_system.length, + # nbody_system.length, + # nbody_system.speed, + # nbody_system.speed, + # nbody_system.speed, + # nbody_system.time, + # nbody_system.length, + # ), + # ( + # handler.INDEX, + # handler.ERROR_CODE, + # ) + # ) + # handler.add_method( + # "get_state_star", + # ( + # handler.INDEX, + # ), + # ( + # nbody_system.mass, + # nbody_system.length, + # nbody_system.length, + # nbody_system.length, + # nbody_system.speed, + # nbody_system.speed, + # nbody_system.speed, + # nbody_system.time, + # nbody_system.length, + # handler.ERROR_CODE + # ) + # ) + # handler.add_method( + # "set_state_star", + # ( + # handler.INDEX, + # nbody_system.mass, + # nbody_system.length, + # nbody_system.length, + # nbody_system.length, + # nbody_system.speed, + # nbody_system.speed, + # nbody_system.speed, + # nbody_system.time, + # nbody_system.length, + # ), + # ( + # handler.ERROR_CODE, + # ) + # ) + # handler.add_method( + # "set_star_tform", + # ( + # handler.INDEX, + # nbody_system.time, + # ), + # ( + # handler.ERROR_CODE, + # ) + # ) + # handler.add_method( + # "get_star_tform", + # ( + # handler.INDEX, + # ), + # ( + # nbody_system.time, + # handler.ERROR_CODE + # ) + # ) + + # handler.add_method( + # 'get_hydro_state_at_point', + # ( + # nbody_system.length, nbody_system.length, + # nbody_system.length, nbody_system.speed, nbody_system.speed, + # nbody_system.speed + # ), + # ( + # nbody_system.density, nbody_system.momentum_density, + # nbody_system.momentum_density, nbody_system.momentum_density, + # nbody_system.energy_density, handler.ERROR_CODE + # ) + # ) + + handler.add_method( + "get_eps2", + ( + ), + ( + nbody_system.length * nbody_system.length, + handler.ERROR_CODE, + ) + ) + + handler.add_method( + "set_eps2", + ( + nbody_system.length * nbody_system.length, + ), + ( + handler.ERROR_CODE, + ) + ) + + # handler.add_method( + # "get_dtime", + # (), + # (nbody_system.time, handler.ERROR_CODE,) + # ) + + # handler.add_method( + # "set_dtime", + # (nbody_system.time, ), + # (handler.ERROR_CODE,) + # ) + + # handler.add_method( + # "get_verbosity", + # (), + # (handler.NO_UNIT, handler.ERROR_CODE,) + # ) + # + # handler.add_method( + # "set_verbosity", + # (handler.NO_UNIT, ), + # (handler.ERROR_CODE,) + # ) + + # handler.add_method( + # "get_pboxsize", + # (), + # (nbody_system.length, handler.ERROR_CODE,) + # ) + # + # handler.add_method( + # "set_pboxsize", + # (nbody_system.length, ), + # (handler.ERROR_CODE,) + # ) + + # handler.add_method( + # "get_gamma", + # (), + # (handler.NO_UNIT, handler.ERROR_CODE,) + # ) + # + # handler.add_method( + # "set_gamma", + # (handler.NO_UNIT, ), + # (handler.ERROR_CODE,) + # ) + + # handler.add_method( + # "get_alpha", + # (), + # (handler.NO_UNIT, handler.ERROR_CODE,) + # ) + # + # handler.add_method( + # "set_alpha", + # (handler.NO_UNIT, ), + # (handler.ERROR_CODE,) + # ) + + # handler.add_method( + # "get_beta", + # (), + # (handler.NO_UNIT, handler.ERROR_CODE,) + # ) + + # handler.add_method( + # "set_beta", + # (handler.NO_UNIT, ), + # (handler.ERROR_CODE,) + # ) + # handler.add_method( + # "get_thermal_energy", + # (), + # (nbody_system.energy, handler.ERROR_CODE,) + # ) + # + # handler.add_method( + # "get_total_energy", + # (), + # (nbody_system.energy, handler.ERROR_CODE,) + # ) + + # self.stopping_conditions.define_methods(handler) diff --git a/src/amuse/community/mizuki/src/mizuki/Makefile b/src/amuse/community/mizuki/src/mizuki/Makefile new file mode 100644 index 0000000000..cf8cce84a6 --- /dev/null +++ b/src/amuse/community/mizuki/src/mizuki/Makefile @@ -0,0 +1,62 @@ +PS_PATH = ../../../src/ +INC = -I$(PS_PATH) + +#CXX = g++ +CXX = mpicxx +CXXFLAGS = -std=c++17 -O3 -ffast-math -funroll-loops -Wall +CXXFLAGS += -DPARTICLE_SIMULATOR_THREAD_PARALLEL -fopenmp +# CXXFLAGS += -DPARTICLE_SIMULATOR_MPI_PARALLEL + +# fdps-autotest-set-vars (DO NOT CHANGE THIS LINE) + +# Simulation control macros +CXXFLAGS += -DENABLE_VARIABLE_SMOOTHING_LENGTH +CXXFLAGS += -DUSE_ENTROPY +CXXFLAGS += -DUSE_BALSARA_SWITCH +#CXXFLAGS += -DUSE_PRESCR_OF_THOMAS_COUCHMAN_1992 +CXXFLAGS += -DISOTHERMAL_EOS + +CXXFLAGS += -DINITIAL_CONDITION=0 +#CXXFLAGS += -DINITIAL_CONDITION=1 +#CXXFLAGS += -DINITIAL_CONDITION=2 +#CXXFLAGS += -DINITIAL_CONDITION=3 + +use_phantom_grape_x86 = yes + +# fdps-autotest-set-vars (DO NOT CHANGE THIS LINE) + +ifeq ($(use_phantom_grape_x86),yes) +PG_ROOT = $(PS_PATH)/phantom_grape_x86/G5/newton/libpg5 +INC += -I$(PG_ROOT) +CXXFLAGS += -DENABLE_PHANTOM_GRAPE_X86 +LIBS = -L$(PG_ROOT) -lpg5 +PG_BUILD = cd $(PG_ROOT) && $(MAKE) distclean libpg5.a +PG_CLEAN = cd $(PG_ROOT) && $(MAKE) distclean +else +PG_BUILD = +PG_CLEAN = +endif + +CPPOBJS = $(patsubst %.cpp, %.o, $(wildcard *.cpp)) +CPPHDRS = $(wildcard *.h) +PROGRAM = mizuki.out + +.PHONY: clean all + +all: $(CPPOBJS) $(CPPHDRS) + $(PG_BUILD) + $(CXX) $(CXXFLAGS) $(CPPOBJS) -o $(PROGRAM) $(LIBS) $(INC) + +%.o: %.cpp $(CPPHDRS) + $(CXX) -c $< $(CXXFLAGS) $(INC) + +clean: + rm -f $(CPPOBJS) + +distclean: clean + $(PG_CLEAN) + rm -f $(PROGRAM) + rm -rf result + +# fdps-autotest-run (DO NOT CHANGE THIS LINE) + diff --git a/src/amuse/community/mizuki/src/mizuki/leapfrog.hpp b/src/amuse/community/mizuki/src/mizuki/leapfrog.hpp new file mode 100644 index 0000000000..4371b42084 --- /dev/null +++ b/src/amuse/community/mizuki/src/mizuki/leapfrog.hpp @@ -0,0 +1,62 @@ + +/* Leapfrog integrators */ +void InitialKick(PS::ParticleSystem & psys, const PS::F64 dt){ + for (PS::S32 i = 0; i < psys.getNumberOfParticleLocal(); i++){ + psys[i].vel += 0.5 * dt * psys[i].acc; + } +} + +void InitialKick(PS::ParticleSystem & psys, const PS::F64 dt){ + for (PS::S32 i = 0; i < psys.getNumberOfParticleLocal(); i++){ + psys[i].vel_half = psys[i].vel + 0.5 * dt * (psys[i].acc_grav + psys[i].acc_hydro); +#if !defined(ISOTHERMAL_EOS) + psys[i].eng_half = psys[i].eng + 0.5 * dt * psys[i].eng_dot; + psys[i].ent_half = psys[i].ent + 0.5 * dt * psys[i].ent_dot; +#endif + } +} + + +void FullDrift(PS::ParticleSystem& psys, const PS::F64 dt){ + for (PS::S32 i = 0; i < psys.getNumberOfParticleLocal(); i++){ + psys[i].pos += dt * psys[i].vel; + } +} +void FullDrift(PS::ParticleSystem& psys, const PS::F64 dt){ + for (PS::S32 i = 0; i < psys.getNumberOfParticleLocal(); i++){ + psys[i].pos += dt * psys[i].vel_half; + } +} + + +void Predict(PS::ParticleSystem& psys, const PS::F64 dt){ + for(PS::S32 i = 0; i < psys.getNumberOfParticleLocal(); i++){ + psys[i].vel += dt * (psys[i].acc_grav + psys[i].acc_hydro); +#if !defined(ISOTHERMAL_EOS) + psys[i].eng += dt * psys[i].eng_dot; + psys[i].ent += dt * psys[i].ent_dot; +#endif + } +} + + +void FinalKick(PS::ParticleSystem& psys, const PS::F64 dt){ + for (PS::S32 i = 0; i < psys.getNumberOfParticleLocal(); i++){ + psys[i].vel += 0.5 * dt * psys[i].acc; + } +} +void FinalKick(PS::ParticleSystem& psys, const PS::F64 dt){ + for (PS::S32 i = 0; i < psys.getNumberOfParticleLocal(); i++){ + psys[i].vel = psys[i].vel_half + 0.5 * dt * (psys[i].acc_grav + psys[i].acc_hydro); +#if !defined(ISOTHERMAL_EOS) + psys[i].eng = psys[i].eng_half + 0.5 * dt * psys[i].eng_dot; + psys[i].ent = psys[i].ent_half + 0.5 * dt * psys[i].ent_dot; +#endif + } +#if defined(DUMP_VELOCITY_OF_SPH_PARTICLE) + const PS::F64 coeff = 0.1; + for (PS::S32 i = 0; i < psys.getNumberOfParticleLocal(); i++) { + psys[i].vel *= std::exp(- coeff * (CFL_hydro/0.1) * psys[i].snds * dt / psys[i].smth); + } +#endif +} diff --git a/src/amuse/community/mizuki/src/mizuki/main.cpp b/src/amuse/community/mizuki/src/mizuki/main.cpp new file mode 100644 index 0000000000..fcc6da7ba2 --- /dev/null +++ b/src/amuse/community/mizuki/src/mizuki/main.cpp @@ -0,0 +1,77 @@ +// Include the standard C++ headers +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +// Include the header file of FDPS +#include +// Include the header file of Phantom-GRAPE library +#if defined(ENABLE_PHANTOM_GRAPE_X86) +#include +#endif + +#include "mathematical_constants.h" +#include "physical_constants.h" +#include "user_defined.hpp" +#include "mizuki.hpp" + +static Mizuki* mizuki=NULL; + +int main(int argc, char* argv[]){ + // Initialize FDPS + mizuki = new Mizuki; + PS::S64 id; + PS::F64 mass; + PS::F64 x; + PS::F64 y; + PS::F64 z; + PS::F64 vx; + PS::F64 vy; + PS::F64 vz; + PS::F64 eng; + + mizuki->initialize(argc, argv); + mizuki->set_defaults(); + + // Populate SPH with AMUSE here + PS::S32 N_sph = 10000; + for (PS::S32 i = 0; i < N_sph; i++) { + mass = 1; + x = i; + y = 0.; + z = i / 2.; + vx = 0.; + vy = 0.; + vz = 0.; + eng = 100.; + mizuki->add_hydro_particle( + i, mass, x, y, z, vx, vy, vz, eng + ); + } + PS::F64vec pos = mizuki->get_position(50); + x = pos.x; + y = pos.y; + z = pos.z; + std::cout << id << " " << x << " " << y << " " << z << std::endl; + + // Commit particles + mizuki->commit_particles(); + + std::cout << "There are " << mizuki->get_number_of_particles_sph() + << " particles!" << std::endl; + + // Finalize FDPS + mizuki->finalize(); + delete mizuki; + mizuki = NULL; + return 0; +} diff --git a/src/amuse/community/mizuki/src/mizuki/mathematical_constants.h b/src/amuse/community/mizuki/src/mizuki/mathematical_constants.h new file mode 100644 index 0000000000..15e78ca015 --- /dev/null +++ b/src/amuse/community/mizuki/src/mizuki/mathematical_constants.h @@ -0,0 +1,8 @@ +#pragma once + +namespace mathematical_constants { + +extern const double pi; // circular constant + +} +namespace math_const = mathematical_constants; diff --git a/src/amuse/community/mizuki/src/mizuki/mizuki.hpp b/src/amuse/community/mizuki/src/mizuki/mizuki.hpp new file mode 100644 index 0000000000..885ed2818d --- /dev/null +++ b/src/amuse/community/mizuki/src/mizuki/mizuki.hpp @@ -0,0 +1,470 @@ +// Include the standard C++ headers +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +// Include the header file of FDPS +#include +// Include the header file of Phantom-GRAPE library +#if defined(ENABLE_PHANTOM_GRAPE_X86) +#include +#endif + +#include "leapfrog.hpp" +#include "mathematical_constants.h" +//#include "user_defined.hpp" + +void calcDensity(PS::ParticleSystem & psys, + PS::DomainInfo & dinfo, + PS::TreeForForceShort::Gather & tree) { +#if defined(ENABLE_VARIABLE_SMOOTHING_LENGTH) + const PS::S32 n_loc = psys.getNumberOfParticleLocal(); + const PS::S64 n_glb = psys.getNumberOfParticleGlobal(); + // Determine the density and the smoothing length so that Eq.(6) in Springel (2005) + // holds within a specified accuracy. + SCF_smth = 1.25; + PS::S32 iter = 0; + for (;;) { + iter++; + // Compute density, etc. + tree.calcForceAllAndWriteBack(CalcDensity(), psys, dinfo); + // Check convergence + PS::S32 n_compl_loc = 0; + for (PS::S32 i = 0; i < n_loc; i++) { + if (psys[i].flag == 1) n_compl_loc++; + } + const PS::S64 n_compl = PS::Comm::getSum(n_compl_loc); + if (n_compl == n_glb) break; + } + // Reset SCF_smth + SCF_smth = 1.0; +#else + SCF_smth = 1.0; + tree.calcForceAllAndWriteBack(CalcDensity(), psys, dinfo); +#endif +} + +void setEntropy(PS::ParticleSystem& psys){ + for (PS::S32 i = 0; i < psys.getNumberOfParticleLocal(); i++) { + psys[i].setEntropy(); + } +} + +void setPressure(PS::ParticleSystem& psys){ + for (PS::S32 i = 0; i < psys.getNumberOfParticleLocal(); i++) { + psys[i].setPressure(); + } +} + + +class Mizuki { +private: + PS::S64 nstep; + PS::S64 id_offset; + PS::S64 N_sph_max; + PS::F32 theta_gravity; + PS::BOUNDARY_CONDITION boundary_condition; + PS::F64ort pos_root_domain; + PS::S32 nprocs; + PS::S32 nthrds; + //PS::ParticleSystem psys_test; + PS::DomainInfo dinfo; + PS::S64 numPtclSPH; + PS::S64 numPtclAll; + PS::TreeForForceLong::Monopole tree_grav; + PS::TreeForForceShort::Gather tree_dens; + PS::TreeForForceShort::Symmetry tree_hydro; + PS::F64 dt; + bool enable_gravity_interact = true; + bool enable_hydro_interact = true; + bool use_entropy = false; + +public: + PS::ParticleSystem psys_nbody; + PS::ParticleSystem psys_sph; + PS::F64 time; + PS::F64 epsilon_gravity; + //std::str equation_of_state = "isothermal"; + PS::F64 dt_max; + PS::F64 dt_min; + Mizuki() { + this->boundary_condition = PS::BOUNDARY_CONDITION_OPEN; + this->nstep = 0; + this->time = 0.; + this->dt_max = 1./64.; // 2**6 + this->dt_min = 1./1048576.; // 2**20 + this->id_offset = 0; + this->epsilon_gravity = 1.0e-4; + this->theta_gravity = 0.5; + this->N_sph_max = 100000; + } + + void set_defaults(){ + //PS::F32 theta_grav = 0.5; + //this->epsilon_gravity = 1.0e-3; + } + + void initialize(){ + int argc = 0; + char **argv = NULL; + PS::Initialize(argc, argv); + this->psys_nbody.initialize(); + this->psys_sph.initialize(); + //psys_test.initialize(); + this->psys_nbody.setNumberOfParticleLocal(0); + this->psys_sph.setNumberOfParticleLocal(0); + //psys_test.setNumberOfParticleLocal(0); + this->dinfo.initialize(); + this->nprocs = PS::Comm::getNumberOfProc(); + this->nthrds = PS::Comm::getNumberOfThread(); + if (PS::Comm::getRank() == 0) { + std::cout << "=================================" << std::endl + << " This is Mizuki, " << std::endl + << " we are using FDPS " << std::endl + << " # of processes is " << this->nprocs << std::endl + << " # of thread is " << this->nthrds << std::endl + << "=================================" << std::endl; + } + } + + void checkConservativeVariables() { + PS::F64 ekin_loc = 0.0; + PS::F64 epot_loc = 0.0; + PS::F64 eth_loc = 0.0; + PS::F64vec mom_loc = 0.0; + for (PS::S32 i = 0; i < this->psys_nbody.getNumberOfParticleLocal(); i++) { + ekin_loc += 0.5 * this->psys_nbody[i].mass * this->psys_nbody[i].vel * this->psys_nbody[i].vel; + epot_loc += 0.5 * this->psys_nbody[i].mass * (this->psys_nbody[i].pot + this->psys_nbody[i].mass / this->epsilon_gravity); + mom_loc += this->psys_nbody[i].mass * this->psys_nbody[i].vel; + } + for (PS::S32 i = 0; i < this->psys_sph.getNumberOfParticleLocal(); i++) { + ekin_loc += 0.5 * this->psys_sph[i].mass * this->psys_sph[i].vel * this->psys_sph[i].vel; + epot_loc += 0.5 * this->psys_sph[i].mass * (this->psys_sph[i].pot_grav + this->psys_sph[i].mass / this->epsilon_gravity); + eth_loc += this->psys_sph[i].mass * this->psys_sph[i].eng; + mom_loc += this->psys_sph[i].mass * this->psys_sph[i].vel; + } + PS::F64 ekin = PS::Comm::getSum(ekin_loc); + PS::F64 epot = PS::Comm::getSum(epot_loc); + PS::F64 eth = PS::Comm::getSum(eth_loc); + PS::F64vec mom = PS::Comm::getSum(mom_loc); + + static bool is_initialized = false; + static PS::F64 emech_ini, etot_ini; + if (is_initialized == false) { + emech_ini = ekin + epot; + etot_ini = ekin + epot + eth; + is_initialized = true; + } + + if (PS::Comm::getRank() == 0){ + const PS::F64 emech = ekin + epot; + const PS::F64 etot = ekin + epot + eth; + const PS::F64 relerr_mech = std::fabs((emech - emech_ini)/emech_ini); + const PS::F64 relerr_tot = std::fabs((etot - etot_ini)/etot_ini); + std::cout << "-------------------------" << std::endl; + std::cout << "E_kin = " << ekin << std::endl; + std::cout << "E_pot = " << epot << std::endl; + std::cout << "E_th = " << eth << std::endl; + std::cout << "E_mech = " << emech << " (" << relerr_mech << ")" << std::endl; + std::cout << "E_tot = " << etot << " (" << relerr_tot << ")" << std::endl; + std::cout << "Mom (x) = " << mom.x << std::endl; + std::cout << "Mom (y) = " << mom.y << std::endl; + std::cout << "Mom (z) = " << mom.z << std::endl; + std::cout << "-------------------------" << std::endl; + } + } + + void getTimeStep(){ + this->dt = DBL_MAX; + if (this->dt_max > 0.0) this->dt = this->dt_max; + + // Timescale for N-body system + if (this->enable_gravity_interact) { + for (PS::S32 i = 0; i < this->psys_nbody.getNumberOfParticleLocal(); i++) { + const PS::F64 acc = std::sqrt(this->psys_nbody[i].acc * this->psys_nbody[i].acc); + if (acc > 0.0) + this->dt = std::min(this->dt, CFL_dyn * std::sqrt(this->epsilon_gravity / acc)); + } + } + + // Timescale for SPH system + for (PS::S32 i = 0; i < this->psys_sph.getNumberOfParticleLocal(); i++) { + if (this->enable_gravity_interact){ + const PS::F64 acc = std::sqrt((psys_sph[i].acc_grav + psys_sph[i].acc_hydro) + *(psys_sph[i].acc_grav + psys_sph[i].acc_hydro)); + if (acc > 0.0) + this->dt = std::min(this->dt, CFL_dyn * std::sqrt(epsilon_gravity / acc)); + } + if (this->enable_hydro_interact){ + this->dt = std::min(this->dt, psys_sph[i].dt); + } + } + this->dt = std::max(this->dt, this->dt_min); + this->dt = PS::Comm::getMinValue(this->dt); + } + + PS::S64 get_number_of_particles_sph(){ + return id_offset; + } + + PS::S64 add_sph_particle( + PS::F64 mass, + PS::F64 x, + PS::F64 y, + PS::F64 z, + PS::F64 vx, + PS::F64 vy, + PS::F64 vz, + PS::F64 eng + ) + { + const PS::S64 N_sph = psys_sph.getNumberOfParticleLocal(); + FP_sph p; + //psys_sph.setNumberOfParticleLocal(N_sph+1); + p.id = id_offset; + p.mass = mass; + p.pos.x = x; + p.pos.y = y; + p.pos.z = z; + //p.vel = 0.0; + p.vel.x = vx; + p.vel.y = vy; + p.vel.z = vz; + p.eng = eng; + p.acc_grav = 0.0; + //p.acc_grav.x = 0.; + //p.acc_grav.y = 0.; + //p.acc_grav.z = 0.; + p.pot_grav = 0.0; + p.acc_hydro = 0.0; + //p.acc_hydro.x = 0.; + //p.acc_hydro.y = 0.; + //p.acc_hydro.z = 0.; + p.flag = 0; + p.dens = 0.0; + p.ent = 0.0; + p.pres = 0.0; + p.smth = 0.1; + p.gradh = 0.0; + p.divv = 0.0; + p.rotv = 0.0; + //p.rotv = [0,0,0]; + p.BalSW = 0; + p.snds = 0.0; + p.eng_dot = 0.0; + p.ent_dot = 0.0; + p.dt = this->dt; + p.vel_half = 0.0; + //p.vel_half = [0,0,0]; + p.eng_half = 0; + p.ent_half = 0; + this->psys_sph.addOneParticle(p); + this->id_offset += 1; + //std::cout << id << " " << x << " " << y << " " << z << std::endl; + return p.id; + } + + PS::F64 get_mass( + PS::S64 id + ){ + return psys_sph[id].mass; + } + + PS::F64vec get_position( + PS::S64 id + ){ + return psys_sph[id].pos; + } + + PS::F64vec get_velocity( + PS::S64 id + ){ + return psys_sph[id].vel; + } + + void commit_particles(){ + // Broadcast + PS::Comm::broadcast(&this->boundary_condition, 1, 0); + PS::Comm::broadcast(&this->pos_root_domain, 1, 0); + PS::Comm::broadcast(&this->epsilon_gravity, 1, 0); + //PS::Comm::broadcast(&dt_dump, 1, 0); + //PS::Comm::broadcast(&time_dump, 1, 0); + //PS::Comm::broadcast(&time_end, 1, 0); + PS::Comm::broadcast(&dt_max, 1, 0); + + // Set the boundary condition and the size of the computational domain if needed. + this->dinfo.setBoundaryCondition(boundary_condition); + if (boundary_condition != PS::BOUNDARY_CONDITION_OPEN) { + this->dinfo.setPosRootDomain(pos_root_domain.low_, + pos_root_domain.high_); + } + + // Compute the average mass of SPH particles + PS::F64 m_sum_loc = 0.0; + for (PS::S32 i = 0; i < psys_sph.getNumberOfParticleLocal(); i++) + m_sum_loc += psys_sph[i].getCharge(); + mass_avg = PS::Comm::getSum(m_sum_loc) / psys_sph.getNumberOfParticleGlobal(); + + + // Perform domain decomposition + // TODO: what does 'false' mean here? clear? + this->dinfo.collectSampleParticle(this->psys_nbody,false); + this->dinfo.collectSampleParticle(this->psys_sph,false); + this->dinfo.decomposeDomain(); + + // Perform particle exchange + this->psys_nbody.exchangeParticle(this->dinfo); + this->psys_sph.exchangeParticle(this->dinfo); + + + // Make tree structures + numPtclSPH = std::max(psys_sph.getNumberOfParticleLocal(),1); + numPtclAll = psys_nbody.getNumberOfParticleLocal() + numPtclSPH; + //numPtclAll = std::max(psys_sph.getNumberOfParticleLocal(),1); + //std::cout << "N_sph: " << psys_sph.getNumberOfParticleLocal() << std::endl; + + //PS::TreeForForceLong::Monopole this->tree_grav; + this->tree_grav.initialize(3 * numPtclAll, theta_gravity); + + //PS::TreeForForceShort::Gather this->tree_dens; + this->tree_dens.initialize(3 * numPtclSPH); + + //PS::TreeForForceShort::Symmetry this->tree_hydro; + this->tree_hydro.initialize(3 * numPtclSPH); + +#if defined(ENABLE_PHANTOM_GRAPE_X86) + g5_open(); + g5_set_eps_to_all(epsilon_gravity); +#endif + + // Peform force calculations + //- Gravity calculations + if (this->enable_gravity_interact){ + tree_grav.setParticleLocalTree(psys_nbody); + tree_grav.setParticleLocalTree(psys_sph,false); + tree_grav.calcForceMakingTree(CalcGravity, + CalcGravity, + dinfo); + for (PS::S32 i = 0; i < psys_nbody.getNumberOfParticleLocal(); i++) { + psys_nbody[i].copyFromForce(tree_grav.getForce(i)); + } + const PS::S32 offset = psys_nbody.getNumberOfParticleLocal(); + for (PS::S32 i = 0; i < psys_sph.getNumberOfParticleLocal(); i++) { + psys_sph[i].copyFromForce(tree_grav.getForce(i+offset)); + //psys_sph[i].copyFromForce(tree_grav.getForce(i)); + } + } + + //- SPH calculations + if (this->enable_hydro_interact){ + calcDensity(psys_sph, dinfo, tree_dens); + if (this->use_entropy){ + setEntropy(psys_sph); + } + setPressure(this->psys_sph); + this->tree_hydro.calcForceAllAndWriteBack(CalcHydroForce(), this->psys_sph, this->dinfo); + } + + // Set the timestep + getTimeStep(); + //checkConservativeVariables(psys_nbody, psys_sph); + //std::cout << "finished commit particles" << std::endl; + } + + void evolve_model(PS::F64 time_end){ + for ( + this->time; + this->time < time_end; + this->time += this->dt, + this->nstep++ + ){ + if (PS::Comm::getRank() == 0) { + std::cout << "nstep = " << this->nstep + << " dt = " << this->dt + << " time = " << this->time + << " time_end = " << time_end + << std::endl; + } + //std::cout << "***** a: " << psys_sph[0].acc_grav << " " << psys_sph[0].acc_hydro << std::endl; + // Leap frog: Initial Kick & Full Drift + //InitialKick(psys_nbody, dt); + InitialKick(this->psys_sph, this->dt); + //FullDrift(psys_nbody, this->dt); + FullDrift(this->psys_sph, this->dt); + if (dinfo.getBoundaryCondition() != PS::BOUNDARY_CONDITION_OPEN) { + //psys_nbody.adjustPositionIntoRootDomain(dinfo); + psys_sph.adjustPositionIntoRootDomain(dinfo); + } + + // Leap frog: Predict + Predict(this->psys_sph, this->dt); + + // Perform domain decomposition again + //dinfo.collectSampleParticle(psys_nbody); + dinfo.collectSampleParticle(psys_sph,false); + dinfo.decomposeDomain(); + + // Exchange the particles between the (MPI) processes + //psys_nbody.exchangeParticle(dinfo); + psys_sph.exchangeParticle(dinfo); + + // Peform force calculations + PS::F64 t_start; + //- Gravity calculations + PS::Comm::barrier(); t_start = PS::GetWtime(); + if (this->enable_gravity_interact){ + //tree_grav.setParticleLocalTree(psys_nbody); + //tree_grav.setParticleLocalTree(psys_sph,false); + tree_grav.setParticleLocalTree(psys_sph,true); + tree_grav.calcForceMakingTree(CalcGravity, + CalcGravity, + dinfo); + //for (PS::S32 i = 0; i < psys_nbody.getNumberOfParticleLocal(); i++) { + // psys_nbody[i].copyFromForce(tree_grav.getForce(i)); + //} + //const PS::S32 offset = psys_nbody.getNumberOfParticleLocal(); + for (PS::S32 i = 0; i < psys_sph.getNumberOfParticleLocal(); i++) { + //psys_sph[i].copyFromForce(tree_grav.getForce(i+offset)); + psys_sph[i].copyFromForce(tree_grav.getForce(i)); + } + } + PS::Comm::barrier(); + //if (PS::Comm::getRank() == 0) std::cout << "t_grav = " << (PS::GetWtime() - t_start) << std::endl; + //- SPH calculations + PS::Comm::barrier(); t_start = PS::GetWtime(); + + if (this->enable_hydro_interact){ + calcDensity(psys_sph, dinfo, tree_dens); + setPressure(psys_sph); + tree_hydro.calcForceAllAndWriteBack(CalcHydroForce(), psys_sph, dinfo); + } + + PS::Comm::barrier(); + //if (PS::Comm::getRank() == 0) std::cout << "t_hydro = " << (PS::GetWtime() - t_start) << std::endl; + + // Get a new timestep + //this->dt = getTimeStep(psys_nbody, psys_sph); + getTimeStep(); + + // Leap frog: Final Kick + //FinalKick(psys_nbody, dt); + FinalKick(psys_sph, this->dt); + } + } + + void finalize(){ +#if defined(ENABLE_PHANTOM_GRAPE_X86) + g5_close(); +#endif + PS::Finalize(); + } +}; diff --git a/src/amuse/community/mizuki/src/mizuki/physical_constants.h b/src/amuse/community/mizuki/src/mizuki/physical_constants.h new file mode 100644 index 0000000000..2e4890cf9b --- /dev/null +++ b/src/amuse/community/mizuki/src/mizuki/physical_constants.h @@ -0,0 +1,23 @@ +#pragma once + +namespace physical_constants { + +// Physical constants +extern const double Ggrav; // Gravitational constants +extern const double kBoltz; // Boltzmann constant + +// Mass units +extern const double Msolar; // Solar mass +extern const double dalton; // Unified atomic mass unit +extern const double Mhydrogen; // Mass of hydrogen atom + +// Length units +extern const double km; // kilo-meters +extern const double AU; // Astronomical unit +extern const double pc, kpc, Mpc, Gpc; // parsec + +// Time units +extern const double yr, kyr, Myr, Gyr; // year + +} +namespace phys_const = physical_constants; diff --git a/src/amuse/community/mizuki/src/mizuki/user_defined.hpp b/src/amuse/community/mizuki/src/mizuki/user_defined.hpp new file mode 100644 index 0000000000..76bbf9a1ab --- /dev/null +++ b/src/amuse/community/mizuki/src/mizuki/user_defined.hpp @@ -0,0 +1,667 @@ +// Include the standard C++ headers +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +// Include the header file of FDPS +#include +// Include the header file of Phantom-GRAPE library +#if defined(ENABLE_PHANTOM_GRAPE_X86) +#include +#endif + +//#include "mathematical_constants.h" +//#include "physical_constants.h" + + +//const PS::F64 specific_heat_ratio = 5.0/3.0; // FIXME needs to be AMUSE settable +const PS::F64 specific_heat_ratio = 1.0; //5.0/3.0; // FIXME needs to be AMUSE settable +PS::F64 SCF_smth = 1.0; // scale factor for smoothing length +PS::F64 mass_avg; +const PS::F64 CFL_hydro = 0.3; // CFL(Courant-Friedrichs-Lewy) number for hydrodynamics +const PS::F64 alpha_AV = 1.0; // coefficient of artificial viscosity +const PS::S32 N_neighbor = 50; // number of neighbor particles +const PS::F64 CFL_dyn = 0.3; // coefficient used to limit a timestep in terms of dynamics +PS::F64 epsilon_gravity = 0.0001; + +/* Kernel Function */ +PS::F64 W(const PS::F64 r, const PS::F64 h){ + // M4 Cubic spline kernel + // (see Eq. (4) in Springel (2005)[MNRAS,364,1105]) + const PS::F64 u = r/h; + const PS::F64 cc=8.0/(std::numbers::pi*h*h*h); + if (u <= 0.5) { + const PS::F64 u2 = u*u; + return cc*(1.0+u2*6.0*(-1.0+u)); + } + else if ((0.5 < u) && (u <= 1.0)) { + const PS::F64 s = 1.0-u; + return cc*2.0*s*s*s; + } + else { + return 0.0; + } +} + +PS::F64vec gradW(const PS::F64vec dr, const PS::F64 h){ + // This subroutine gives \nabla W(r,h), i.e., + // \dfrac{\partial W(r,h)}{\partial r}\dfrac{dr}{r}. + const PS::F64 r = std::sqrt(dr * dr); + const PS::F64 u=r/h; + const PS::F64 cc = -48.0/(std::numbers::pi*h*h*h*h); +#if defined(USE_PRESCR_OF_THOMAS_COUCHMAN_1992) + if (u <= 1.0/3.0) { + return dr * cc*(1.0/3.0)/(r); + } + else if ((1.0/3.0 < u) && (u <= 0.5)) { + return dr * cc*u*(2.0-3.0*u)/(r); + } + else if ((0.5 < u) && (u < 1.0)) { + return dr * cc*(1.0-u)*(1.0-u)/(r); + } + else { + return PS::F64vec(0.0, 0.0, 0.0); + } +#else + if ((0.0 < u) && (u <= 0.5)) { + return dr * cc*u*(2.0-3.0*u)/(r); + } + else if ((0.5 < u) && (u < 1.0)) { + return dr * cc*(1.0-u)*(1.0-u)/(r); + } + else { + // r=0 case is included in this branch + return PS::F64vec(0.0, 0.0, 0.0); + } +#endif +} + + +PS::F64 dWdh(const PS::F64 r, const PS::F64 h){ + // This subroutine gives dW(r,h)/dh, i.e., + // \dfrac{\partial W(r,h)}{\partial h}. + const PS::F64 u=r/h; + const PS::F64 cc=-24.0/(std::numbers::pi*h*h*h*h); + if (u <= 0.5) { + const PS::F64 u2 = u*u; + return cc*(1.0 + +u2*(-10.0 + +12.0*u)); + } + else if ((0.5 < u) && (u < 1.0)) { + const PS::F64 s = 1.0-u; + return cc*2.0*s*s*(1.0-2.0*u); + } + else { + return 0.0; + } + +} + +//** Force class for gravity calculation +class Force_grav { +public: + PS::F64vec acc; + PS::F64 pot; + void clear() { + acc = 0.0; + pot = 0.0; + } +}; +//** Force classes for SPH calculation +class Force_dens{ +public: + PS::S32 flag; + PS::F64 dens; + PS::F64 smth; + PS::F64 gradh; + PS::F64 divv; + PS::F64vec rotv; + void clear(){ + flag = 0; + dens = 0.0; + gradh = 0.0; + divv = 0.0; + rotv = 0.0; + } +}; +class Force_hydro{ +public: + PS::F64vec acc; + PS::F64 eng_dot; + PS::F64 ent_dot; + PS::F64 dt; + void clear(){ + acc = 0.0; + eng_dot = 0.0; + ent_dot = 0.0; + } +}; + +class FP_test { +public: + PS::S64 id; + PS::F64vec pos; + PS::F64vec vel; + PS::F64vec acc_grav; // gravitational acceleration + PS::F64 pot_grav; // gravitational potential + PS::F64vec acc_hydro; // acceleration due to pressure-gradient + + PS::F64 dens; // mass density + PS::F64 eng; // specific internal energy + //PS::F64 pres; // pressure + PS::F64 smth; // smoothing length + PS::F64 gradh; // grad-h term + PS::F64 divv; // divergence of velocity + PS::F64vec rotv; // rotation of velocity + //PS::F64 snds; // sound speed + + PS::F64vec getPos() const { + return pos; + } + + void setPos(const PS::F64vec& pos){ + this->pos = pos; + } + + PS::F64 getCharge() const { + return 0; + } + + void copyFromForce(const Force_grav& f) { + this->acc_grav = f.acc; + this->pot_grav = f.pot; + } + + void copyFromForce(const Force_dens& f){ + //this->flag = f.flag; + this->dens = f.dens; + this->smth = f.smth; + this->gradh = f.gradh; + this->divv = f.divv; + this->rotv = f.rotv; + + } + + void copyFromForce(const Force_hydro& f){ + this->acc_hydro = f.acc; + //this->eng_dot = f.eng_dot; + //this->ent_dot = f.ent_dot; + //this->dt = f.dt; + } + +}; + +class FP_sph { +public: + PS::S64 id; + PS::F64 mass; + PS::F64vec pos; + PS::F64vec vel; + PS::F64vec acc_grav; // gravitational acceleration + PS::F64 pot_grav; // gravitational potential + PS::F64vec acc_hydro; // acceleration due to pressure-gradient + PS::S32 flag; + PS::F64 dens; // mass density + PS::F64 eng; // specific internal energy + PS::F64 ent; // entropy + PS::F64 pres; // pressure + PS::F64 smth; // smoothing length + PS::F64 gradh; // grad-h term + PS::F64 divv; // divergence of velocity + PS::F64vec rotv; // rotation of velocity + PS::F64 BalSW; // Balsara switch + PS::F64 snds; // sound speed + PS::F64 eng_dot; // time rate of change of `eng` + PS::F64 ent_dot; // time rate of change of `ent` + PS::F64 dt; // hydrodynamic time step for this particle + PS::F64vec vel_half; + PS::F64 eng_half; + PS::F64 ent_half; + + void copyFromForce(const Force_grav& f) { + this->acc_grav = f.acc; + this->pot_grav = f.pot; + } + void copyFromForce(const Force_dens& f){ + this->flag = f.flag; + this->dens = f.dens; + this->smth = f.smth; + this->gradh = f.gradh; + this->divv = f.divv; + this->rotv = f.rotv; + + } + void copyFromForce(const Force_hydro& f){ + this->acc_hydro = f.acc; + this->eng_dot = f.eng_dot; + this->ent_dot = f.ent_dot; + this->dt = f.dt; + } + PS::F64 getCharge() const{ + return this->mass; + } + PS::F64vec getPos() const{ + return this->pos; + } + PS::F64 getRSearch() const{ + return this->smth; + } + void setPos(const PS::F64vec& pos){ + this->pos = pos; + } + void writeAscii(FILE* fp) const{ + fprintf(fp, + "%lld\t%lf\t%lf\t%lf\t%lf\t%lf\t" + "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", + this->id, this->mass, + this->pos.x, this->pos.y, this->pos.z, + this->vel.x, this->vel.y, this->vel.z, + this->dens, this->eng, this->ent, this->pres); + } + void readAscii(FILE* fp){ + fscanf(fp, + "%lld\t%lf\t%lf\t%lf\t%lf\t%lf\t" + "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", + &this->id, &this->mass, + &this->pos.x, &this->pos.y, &this->pos.z, + &this->vel.x, &this->vel.y, &this->vel.z, + &this->dens, &this->eng, &this->ent, &this->pres); + } + void writeBinaryPos(FILE* fp) const { + fwrite(&this->pos, sizeof(this->pos), 1, fp); + } + void setEntropy(){ + ent = (specific_heat_ratio - 1.0) * eng / std::pow(dens, specific_heat_ratio - 1.0); + } + void setPressure(){ +#if defined(ISOTHERMAL_EOS) + // In this case, eng = const. + pres = (specific_heat_ratio - 1.0) * dens * eng; + ent = pres / std::pow(dens, specific_heat_ratio); +#else +#if defined(USE_ENTROPY) + pres = ent * std::pow(dens, specific_heat_ratio); + eng = pres / ((specific_heat_ratio - 1.0) * dens); +#else + pres = (specific_heat_ratio - 1.0) * dens * eng; + ent = pres / std::pow(dens, specific_heat_ratio); +#endif +#endif + snds = std::sqrt(specific_heat_ratio * pres / dens); +#if defined(USE_BALSARA_SWITCH) + BalSW = std::fabs(divv) / (std::fabs(divv) + std::sqrt(rotv * rotv) + 1.0e-4 * snds / smth); +#else + BalSW = 1.0; +#endif + } + +}; + +//** Full Particle Classes +class FP_nbody{ +public: + PS::S64 id; + PS::F64 mass; + PS::F64vec pos; + PS::F64vec vel; + PS::F64vec acc; + PS::F64 pot; + + PS::F64vec getPos() const { + return pos; + } + PS::F64 getCharge() const { + return mass; + } + void setPos(const PS::F64vec& pos){ + this->pos = pos; + } + void copyFromForce(const Force_grav & f) { + this->acc = f.acc; + this->pot = f.pot; + } + void writeAscii(FILE* fp) const { + fprintf(fp, "%lld\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", + this->id, this->mass, + this->pos.x, this->pos.y, this->pos.z, + this->vel.x, this->vel.y, this->vel.z); + } + void readAscii(FILE* fp) { + fscanf(fp, "%lld\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", + &this->id, &this->mass, + &this->pos.x, &this->pos.y, &this->pos.z, + &this->vel.x, &this->vel.y, &this->vel.z); + } +}; + +//** Essential Particle Class +class EP_grav { +public: + PS::S64 id; + PS::F64 mass; + PS::F64vec pos; + PS::F64 smth; + + void copyFromFP(const FP_nbody& fp) { + this->id = fp.id; + this->mass = fp.mass; + this->pos = fp.pos; + //this->smth = fp.smth; -> fixed epsilon_gravity in this case? + } + void copyFromFP(const FP_sph& fp) { + this->id = fp.id; + this->mass = fp.mass; + this->pos = fp.pos; + this->smth = fp.smth; + } + PS::F64 getCharge() const { + return this->mass; + } + PS::F64vec getPos() const { + return this->pos; + } +}; + +class EP_hydro { +public: + PS::S64 id; + PS::F64vec pos; + PS::F64vec vel; + PS::F64 mass; + PS::F64 smth; + PS::F64 dens; + PS::F64 pres; + PS::F64 gradh; + PS::F64 snds; + PS::F64 BalSW; + + void copyFromFP(const FP_sph& fp){ + this->id = fp.id; + this->pos = fp.pos; + this->vel = fp.vel; + this->mass = fp.mass; + this->smth = fp.smth; + this->dens = fp.dens; + this->pres = fp.pres; + this->gradh = fp.gradh; + this->snds = fp.snds; + this->BalSW = fp.BalSW; + } + PS::F64vec getPos() const{ + return this->pos; + } + PS::F64 getRSearch() const{ + return SCF_smth * this->smth; + } + void setPos(const PS::F64vec& pos){ + this->pos = pos; + } +}; + +class CalcDensity{ +public: + void operator () (const EP_hydro * ep_i, + const PS::S32 n_ip, + const EP_hydro * ep_j, + const PS::S32 n_jp, + Force_dens * force){ +#if defined(ENABLE_VARIABLE_SMOOTHING_LENGTH) + const PS::F64 eps = 1.0e-6; + const PS::F64 M_trgt = mass_avg * N_neighbor; + for (PS::S32 i = 0; i < n_ip; i++) { + PS::F64 dens = 0.0; + PS::F64 h = ep_i[i].smth; + const PS::F64 h_max_alw = SCF_smth * h; // maximum allowance + PS::F64 h_L = 0.0; + PS::F64 h_U = h_max_alw; + PS::F64 dh_prev = 0.0; + PS::S32 n_unchanged = 0; + // Software caches + PS::F64 * mj = (PS::F64 *)malloc(sizeof(PS::F64) * n_jp); + PS::F64 * rij = (PS::F64 *)malloc(sizeof(PS::F64) * n_jp); + for (PS::S32 j = 0; j < n_jp; j++) { + mj[j] = ep_j[j].mass; + const PS::F64vec dr = ep_j[j].pos - ep_i[i].pos; + rij[j] = std::sqrt(dr * dr); + } + for (;;) { + // Calculate density + dens = 0.0; + for (PS::S32 j = 0; j < n_jp; j++) { + dens += mj[j] * W(rij[j], h); + } + // Check if the current value of the smoohting length satisfies + // Eq.(5) in Springel (2005). + const PS::F64 M = 4.0 * std::numbers::pi * h * h * h * dens / 3.0; + if ((h < h_max_alw) && (std::abs(M/M_trgt - 1.0) < eps)) { + // In this case, Eq.(5) holds within a specified accuracy. + force[i].flag = 1; + force[i].dens = dens; + force[i].smth = h; + break; + } + if (((h == h_max_alw) && (M < M_trgt)) || (n_unchanged == 4)) { + // In this case, we skip this particle forcibly. + // In order to determine consistently the density + // and the smoohting length for this particle, + // we must re-perform calcForceAllAndWriteBack(). + force[i].flag = 0; + force[i].dens = dens; + force[i].smth = h_max_alw; + break; + } + // Update h_L & h_U + if (M < M_trgt) { + if (h_L < h) h_L = h; + } + else if (M_trgt < M) { + if (h < h_U) h_U = h; + } + const PS::F64 dh = h_U - h_L; + if (dh == dh_prev) { + n_unchanged++; + } + else { + dh_prev = dh; + n_unchanged = 0; + } + // Update smoothing length + h = std::pow((3.0 * M_trgt)/(4.0 * std::numbers::pi * dens), 1.0/3.0); + if ((h <= h_L) || (h == h_U)) { + // In this case, we switch to the bisection search. + // The inclusion of '=' in the if statement is very + // important to escape a limit cycle. + h = 0.5 * (h_L + h_U); + } + else if (h_U < h) { + h = h_U; + } + } + // Calculate grad-h term + if (force[i].flag == 1) { + PS::F64 drho_dh = 0.0; + for (PS::S32 j = 0; j < n_jp; j++) { + drho_dh += mj[j] * dWdh(rij[j], h); + } + force[i].gradh = 1.0 / (1.0 + (h * drho_dh) / (3.0 * dens)); + } + else { + force[i].gradh = 1.0; // dummy value + } +#if defined(USE_BALSARA_SWITCH) + // Compute \div v & \rot v for Balsara switch + for (PS::S32 j = 0; j < n_jp; j++) { + const PS::F64vec dr = ep_i[i].pos - ep_j[j].pos; + const PS::F64vec dv = ep_i[i].vel - ep_j[j].vel; + force[i].divv -= mj[j] * dv * gradW(dr, force[i].smth); + force[i].rotv -= mj[j] * dv ^ gradW(dr, force[i].smth); + } + force[i].divv /= force[i].dens; + force[i].rotv /= force[i].dens; +#endif + // Release memory + free(mj); + free(rij); + } +#else + for (PS::S32 i = 0; i < n_ip ; i++){ + // Compute density + for (PS::S32 j = 0; j < n_jp; j++){ + const PS::F64vec dr = ep_j[j].pos - ep_i[i].pos; + const PS::F64 rij = std::sqrt(dr * dr); + force[i].dens += ep_j[j].mass * W(rij, ep_i[i].smth); + } + force[i].smth = ep_i[i].smth; + force[i].gradh = 1.0; +#if defined(USE_BALSARA_SWITCH) + // Compute \div v & \rot v for Balsara switch + for (PS::S32 j = 0; j < n_jp; j++) { + const PS::F64vec dr = ep_i[i].pos - ep_j[j].pos; + const PS::F64vec dv = ep_i[i].vel - ep_j[j].vel; + force[i].divv -= ep_j[j].mass * dv * gradW(dr, force[i].smth); + force[i].rotv -= ep_j[j].mass * dv ^ gradW(dr, force[i].smth); + } + force[i].divv /= force[i].dens; + force[i].rotv /= force[i].dens; +#endif + } +#endif + } +}; + +class CalcHydroForce{ +public: + void operator () (const EP_hydro * ep_i, + const PS::S32 n_ip, + const EP_hydro * ep_j, + const PS::S32 n_jp, + Force_hydro * force){ + for (PS::S32 i = 0; i < n_ip; i++){ + const PS::F64vec pos_i = ep_i[i].pos; + const PS::F64vec vel_i = ep_i[i].vel; + const PS::F64 smth_i = ep_i[i].smth; + const PS::F64 dens_i = ep_i[i].dens; + const PS::F64 pres_i = ep_i[i].pres; + const PS::F64 f_i = ep_i[i].gradh; + const PS::F64 snds_i = ep_i[i].snds; + const PS::F64 povrho2_i = pres_i / (dens_i * dens_i); + PS::F64 v_sig_max = 0.0; + for (PS::S32 j = 0; j < n_jp; j++){ + const PS::F64vec dr = pos_i - ep_j[j].pos; + const PS::F64vec dv = vel_i - ep_j[j].vel; + const PS::F64 w_ij = (dv * dr < 0) ? dv * dr / std::sqrt(dr * dr) : 0; + const PS::F64 v_sig = snds_i + ep_j[j].snds - 3.0 * w_ij; + v_sig_max = std::max(v_sig_max, v_sig); + const PS::F64 AV = - 0.5 * alpha_AV * v_sig * w_ij / (0.5 * (dens_i + ep_j[j].dens)) + * 0.5 * (ep_i[i].BalSW + ep_j[j].BalSW); + const PS::F64vec gradW_i = gradW(dr, smth_i); + const PS::F64vec gradW_j = gradW(dr, ep_j[j].smth); + const PS::F64vec gradW_ij = 0.5 * (gradW_i + gradW_j); + const PS::F64 povrho2_j = ep_j[j].pres / (ep_j[j].dens * ep_j[j].dens); + const PS::F64 f_j = ep_j[j].gradh; + force[i].acc -= ep_j[j].mass * (f_i * povrho2_i * gradW_i + +f_j * povrho2_j * gradW_j + +AV * gradW_ij); + force[i].eng_dot += ep_j[j].mass * (f_i * povrho2_i * gradW_i + +0.5 * AV * gradW_ij) * dv; + force[i].ent_dot += 0.5 * ep_j[j].mass * AV * gradW_ij * dv; + } + const PS::F64 p = specific_heat_ratio - 1.0; + force[i].ent_dot *= p/std::pow(dens_i, p); + force[i].dt = CFL_hydro * 2.0 * ep_i[i].smth / v_sig_max; + } + } +}; + +/* Interaction functions */ +#if defined(ENABLE_PHANTOM_GRAPE_X86) +template +void CalcGravity(const EP_grav * iptcl, + const PS::S32 ni, + const TParticleJ * jptcl, + const PS::S32 nj, + Force_grav * force) { + const PS::S32 nipipe = ni; + const PS::S32 njpipe = nj; + PS::F64 (*xi)[3] = (PS::F64 (*)[3])malloc(sizeof(PS::F64) * nipipe * PS::DIMENSION); + PS::F64 (*ai)[3] = (PS::F64 (*)[3])malloc(sizeof(PS::F64) * nipipe * PS::DIMENSION); + PS::F64 *pi = (PS::F64 * )malloc(sizeof(PS::F64) * nipipe); + PS::F64 (*xj)[3] = (PS::F64 (*)[3])malloc(sizeof(PS::F64) * njpipe * PS::DIMENSION); + PS::F64 *mj = (PS::F64 * )malloc(sizeof(PS::F64) * njpipe); + for(PS::S32 i = 0; i < ni; i++) { + xi[i][0] = iptcl[i].getPos()[0]; + xi[i][1] = iptcl[i].getPos()[1]; + xi[i][2] = iptcl[i].getPos()[2]; + ai[i][0] = 0.0; + ai[i][1] = 0.0; + ai[i][2] = 0.0; + pi[i] = 0.0; + } + for(PS::S32 j = 0; j < nj; j++) { + xj[j][0] = jptcl[j].getPos()[0]; + xj[j][1] = jptcl[j].getPos()[1]; + xj[j][2] = jptcl[j].getPos()[2]; + mj[j] = jptcl[j].getCharge(); + } + PS::S32 devid = PS::Comm::getThreadNum(); + g5_set_xmjMC(devid, 0, nj, xj, mj); + g5_set_nMC(devid, nj); + g5_calculate_force_on_xMC(devid, xi, ai, pi, ni); + for(PS::S32 i = 0; i < ni; i++) { + force[i].acc[0] += ai[i][0]; + force[i].acc[1] += ai[i][1]; + force[i].acc[2] += ai[i][2]; + force[i].pot -= pi[i]; + } + free(xi); + free(ai); + free(pi); + free(xj); + free(mj); +} +#else +template +void CalcGravity (const EP_grav * ep_i, + const PS::S32 n_ip, + const TParticleJ * ep_j, + const PS::S32 n_jp, + Force_grav * force) { + //const PS::F64 eps2 = epsilon_gravity * epsilon_gravity; + PS::F64 eps2 = 0.0; + for(PS::S32 i = 0; i < n_ip; i++){ + eps2 = ep_i[i].smth * ep_i[i].smth; + PS::F64vec xi = ep_i[i].getPos(); + PS::F64vec ai = 0.0; + PS::F64 poti = 0.0; + for(PS::S32 j = 0; j < n_jp; j++){ + PS::F64vec rij = xi - ep_j[j].getPos(); + //PS::F64 eps = ((ep_i[i].smth + ep_j[j].smth)/2.); + //eps2 = eps * eps + PS::F64 r3_inv = rij * rij + eps2; + //PS::F64 r3_inv = rij * rij; + //if (r3_inv < 4*eps2) { + // r3_inv += eps2; + //} + PS::F64 r_inv = 1.0/sqrt(r3_inv); + r3_inv = r_inv * r_inv; + r_inv *= ep_j[j].getCharge(); + r3_inv *= r_inv; + ai -= r3_inv * rij; + poti -= r_inv; + } + force[i].acc += ai; + force[i].pot += poti; + } +} +#endif