diff --git a/Makefile b/Makefile new file mode 100644 index 0000000000..7f18080df2 --- /dev/null +++ b/Makefile @@ -0,0 +1,96 @@ +-include config.mk + +PYTHON ?= python +CLEAN ?= yes + +export PATH := ${PATH}:$(PWD)/bin +export PYTHONPATH := $(PYTHONPATH):$(PWD)/src:$(PWD)/test + +python_version_full := $(wordlist 2,4,$(subst ., ,$(shell $(PYTHON) --version 2>&1))) +python_version_major := $(word 1,${python_version_full}) +python_version_minor := $(word 2,${python_version_full}) +python_version_patch := $(word 3,${python_version_full}) + +all: config.mk + @-mkdir -p test_results + $(PYTHON) setup.py generate_main + $(PYTHON) setup.py build_codes --inplace + +framework: config.mk + @-mkdir -p test_results + $(PYTHON) setup.py generate_main + $(PYTHON) setup.py build_libraries --inplace + +allinbuild: + $(PYTHON) setup.py build + +build: + $(PYTHON) setup.py build + +# should pick up prefix from configure? +install: + $(PYTHON) setup.py install + +docclean: + make -C doc clean + +clean: + $(PYTHON) setup.py clean + $(PYTHON) setup.py clean_codes --inplace + +oclean: + $(PYTHON) setup.py clean + +distclean: + -rm -f src/amuse/config.mk + -rm -f amuse.sh + -rm -f iamuse.sh + -rm -f ibis-deploy.sh + -rm -f bin/amusifier + -rm -rf test_results src/amuse.egg-info + + -rm -f test/*.000 test/fort.* test/perr test/pout test/test.h5 test/*.log + -rm -f test/codes_tests/perr test/codes_tests/pout + -rm -f test/core_tests/plummer_back_100.ini + -rm -f test/test_python_implementation test/twobody + + $(PYTHON) setup.py clean + $(PYTHON) setup.py dist_clean + $(PYTHON) setup.py clean_codes --inplace + $(PYTHON) setup.py dist_clean --inplace + + make -C doc clean + -find ./ -name "*.pyc" -exec rm \{} \; + -find ./ -type d -name "__pycache__" -exec rm -Rf \{} \; + -find ./ -type d -name "ccache" -exec rm -Rf \{} \; + -rm -Rf build + -rm -f config.mk + -rm -f config.log build.log config.status + -rm -f amuse.cfg + -rm -f test*.pickle test.csv + +tests: + $(PYTHON) setup.py tests + +doc: + $(PYTHON) setup.py -q build_latex + +html: + make -C doc html + +latexpdf: + make -C doc latexpdf + +%.code: +ifneq (,$(findstring s,$(MAKEFLAGS))) + $(PYTHON) setup.py build_code --inplace --clean=$(CLEAN) --code-name=$* +else + $(PYTHON) setup.py -v build_code --inplace --clean=$(CLEAN) --code-name=$* +endif + +help: + @echo "brief overview of most important make options:" + @echo "make - build all AMUSE libraries and community codes " + @echo "make .code - clean & build the community code (or matching name*)" + @echo "make clean - clean codes and libraries" + @echo "make distclean - clean codes and libraries and all configuration files" diff --git a/doc/interactive_tutorial/amuserc b/doc/interactive_tutorial/amuserc index 0b081ed9d5..3bfaed74c5 100644 --- a/doc/interactive_tutorial/amuserc +++ b/doc/interactive_tutorial/amuserc @@ -1,2 +1,2 @@ -[output] -printing_strategy=simple +#[output] +#printing_strategy=simple diff --git a/examples/simple/molecular_cloud_chemistry.py b/examples/simple/molecular_cloud_chemistry.py index f0f65595b8..21a5fd26ed 100644 --- a/examples/simple/molecular_cloud_chemistry.py +++ b/examples/simple/molecular_cloud_chemistry.py @@ -16,8 +16,10 @@ from amuse.community.fi.interface import Fi from amuse.community.krome.interface import Krome +from amuse.community.uclchem.interface import UCLchem from amuse.ext.molecular_cloud import molecular_cloud +from amuse.io import write_set_to_file gamma = 1.4 meanmwt = 1.35 | units.amu @@ -27,11 +29,12 @@ def get_n_T_xi(density, u): number_density=density/meanmwt temperature=((gamma - 1) * meanmwt * u / constants.kB) ionrate=ionization_rate*numpy.ones_like(density) - return (number_density, temperature, ionrate) + radfield = 1.0|units.habing + return (number_density, temperature, ionrate, radfield) def update_chem(sph_parts, chem_parts): channel = sph_parts.new_channel_to(chem_parts) - channel.transform(["number_density", "temperature", "ionrate"], get_n_T_xi, ["density", "u"]) + channel.transform(["number_density", "temperature", "ionrate","radfield"], get_n_T_xi, ["density", "u"]) def evolve_sph_with_chemistry(sph, chem, tend): sph.evolve_model(tend) @@ -50,7 +53,7 @@ def run_mc(N=5000, Mcloud=10000. | units.MSun, Rcloud=1. | units.parsec): tff = 1 / (4 * numpy.pi * constants.G * rho_cloud)**0.5 parts.density = rho_cloud - + parts.radfield = 1.0|units.habing update_chem(parts, parts) print("Tcloud:", parts.temperature.max().in_(units.K)) @@ -58,8 +61,9 @@ def run_mc(N=5000, Mcloud=10000. | units.MSun, Rcloud=1. | units.parsec): print("freefall time:", tff.in_(units.Myr)) sph = Fi(conv) - chem = Krome(redirection="none") - + #chem = Krome(redirection="none") + chem = UCLchem() + print(parts) sph.parameters.self_gravity_flag = True sph.parameters.use_hydro_flag = True sph.parameters.isothermal_flag = True @@ -71,9 +75,15 @@ def run_mc(N=5000, Mcloud=10000. | units.MSun, Rcloud=1. | units.parsec): sph.gas_particles.add_particles(parts) chem.particles.add_particles(parts) - + chem.out_species = ["OH", "OCS", "CO", "CS", "CH3OH"] + H2_index = chem.get_index_of_species('H2') + CO_index = chem.get_index_of_species('CO') tnow = sph.model_time + sph_channel = sph.particles.new_channel_to(parts) + chem_channel = chem.particles.new_channel_to(parts) + + f = pyplot.figure() pyplot.ion() pyplot.show() @@ -82,13 +92,15 @@ def run_mc(N=5000, Mcloud=10000. | units.MSun, Rcloud=1. | units.parsec): while i < (end_time / timestep + 0.5): evolve_sph_with_chemistry(sph, chem, i * timestep) tnow = sph.model_time + sph_channel.copy() + chem_channel.copy() print("done with step:", i, tnow.in_(units.Myr)) i += 1 n = (sph.particles.density / meanmwt).value_in(units.cm**-3) - fh2 = chem.particles.abundances[:, chem.species["H2"]] - co = chem.particles.abundances[:, chem.species["CO"]] - + fh2 = chem.particles.abundances[:, H2_index] + co = chem.particles.abundances[:, CO_index] + pyplot.clf() pyplot.loglog(n, fh2, 'r.') pyplot.loglog(n, co, 'g.') @@ -97,7 +109,8 @@ def run_mc(N=5000, Mcloud=10000. | units.MSun, Rcloud=1. | units.parsec): pyplot.xlabel("density (cm**-3)") pyplot.ylabel("H_2,CO abundance") f.canvas.flush_events() - + particles_to_save = parts.copy() + write_set_to_file(particles_to_save, "mol_cloud_uclchem_{}.txt".format(i), format='amuse',overwrite_file=True) print("done. press key to exit") input() diff --git a/src/amuse/amuserc b/src/amuse/amuserc index e69de29bb2..e44bc69eaa 100644 --- a/src/amuse/amuserc +++ b/src/amuse/amuserc @@ -0,0 +1,2 @@ +[channel] +# debugger=gdb diff --git a/src/amuse/community/tupan/src b/src/amuse/community/tupan/src new file mode 160000 index 0000000000..67d3aa103d --- /dev/null +++ b/src/amuse/community/tupan/src @@ -0,0 +1 @@ +Subproject commit 67d3aa103d77248a04e8f112930ba7bdb55024b2 diff --git a/src/amuse/community/uclchem/Makefile b/src/amuse/community/uclchem/Makefile new file mode 100644 index 0000000000..58cd24cd5d --- /dev/null +++ b/src/amuse/community/uclchem/Makefile @@ -0,0 +1,39 @@ +ifeq ($(origin AMUSE_DIR), undefined) + AMUSE_DIR := $(shell amusifier --get-amuse-dir) +endif +-include ${AMUSE_DIR}/config.mk + + +CLASSNAME=UCLchemInterface + +FFLAGS += $(FCFLAGS) -std=legacy -O2 -ffree-line-length-none +INCL = -I./src/ +export FFLAGS + +all: uclchem_worker + +.PHONY: src/libchem.a +src/libchem.a: + $(MAKE) -C src/ libchem.a + +worker_code.f90: interface.py + $(CODE_GENERATOR) --type=f90 $< $(CLASSNAME) -o $@ + +uclchem_worker: src/libchem.a worker_code.f90 interface.o + $(MPIFC) $(FFLAGS) $(INCL) $(FS_FLAGS) $(LDFLAGS) worker_code.f90 interface.o -o $@ src/libchem.a $(LIBS) $(FS_LIBS) $(LIBS) + +interface.o: src/libchem.a + +%.o: %.f90 + $(MPIFC) $(FFLAGS) $(INCL) -c -o $@ $< + +clean: + make -C src/ clean + rm -f *.pyc + rm -f interface.o uclchem_worker.f90 worker_code.f90 + rm -f uclchem_worker + +distclean: clean + + + diff --git a/src/amuse/community/uclchem/__init__.py b/src/amuse/community/uclchem/__init__.py new file mode 100644 index 0000000000..fbdb1dba26 --- /dev/null +++ b/src/amuse/community/uclchem/__init__.py @@ -0,0 +1 @@ +from .interface import UCLchem diff --git a/src/amuse/community/uclchem/interface.f90 b/src/amuse/community/uclchem/interface.f90 new file mode 100644 index 0000000000..0596b56fe6 --- /dev/null +++ b/src/amuse/community/uclchem/interface.f90 @@ -0,0 +1,126 @@ + function initialize_code() result(ret) + use uclchemhelper + integer :: ret + ret=chem_initialize() + end function + + function cleanup_code() result(ret) + use uclchemhelper + integer :: ret + ret=chem_end() + end function + + function commit_particles() result(ret) + use uclchemhelper + integer :: ret + ret=chem_commit_particles() + end function + + function recommit_particles() result(ret) + use uclchemhelper + integer :: ret + ret=chem_commit_particles() + end function + + function commit_parameters() result(ret) + use uclchemhelper + integer :: ret + ret=chem_commit_parameters() + end function + + function new_particle(id,dens,temperature,ionrate,uvrad) result(ret) + use uclchemhelper + integer :: ret,id + double precision :: dens,temperature,ionrate,uvrad + ret=add_particle(id,dens,temperature,ionrate,uvrad) + end function + + function set_state(id,dens,temperature,ionrate,uvrad) result(ret) + use uclchemhelper + integer :: ret,id + double precision :: dens,temperature,ionrate,uvrad + ret=set_particle_state(id,dens,temperature,ionrate,uvrad) + end function + + function get_state(id,dens,temperature,ionrate,uvrad) result(ret) + use uclchemhelper + integer :: ret,id + double precision :: dens,temperature,ionrate,uvrad + ret=get_particle_state(id,dens,temperature,ionrate,uvrad) + end function + + function get_abundance(id,aid,x) result(ret) + use uclchemhelper + integer :: ret,id,aid + double precision :: x + ret=get_particle_abundance(id,aid,x) + end function + + function set_abundance(id,aid,x) result(ret) + use uclchemhelper + integer ret,id,aid + double precision x + ret=set_particle_abundance(id,aid,x) + end function + + function get_firstlast_abundance(first,last) result(ret) + use network + integer :: ret,first,last + first=1 + last=nSpec + ret=0 + end function + + function get_name_of_species(index,s) result(ret) + use network + integer ret,index + character*16 :: ss(nSpec),s + index = index + 1 + if(index.LT.1.OR.index.GT.nSpec) then + ret=-1 + return + endif + ss=specname + s=ss(index) + ret=0 + end function + + function get_index_of_species(s,index) result(ret) + use network + integer ret,index + character*16 s + if (any(specname == s)) then + index = FINDLOC(specname, s, dim=1) - 1 + ret = 0 + else + ret = 1 + end if + end function + + function run_model(dictionary) result(ret) + use uclchemhelper + integer :: ret + character(len=*) :: dictionary(nparticle) + ret=simple_evolution(dictionary, out_species) + end function + + function get_time(time) result(ret) + use uclchemhelper + integer :: ret + double precision :: time + ret=get_current_time(time) + end function + + function delete_particle(id) result(ret) + use uclchemhelper + integer :: ret,id + ret=remove_particle(id) + end function + + function set_species(species) result(ret) + use uclchemhelper + integer :: ret + character(len=*) :: species + out_species = species + ret = 0 + end function diff --git a/src/amuse/community/uclchem/interface.py b/src/amuse/community/uclchem/interface.py new file mode 100644 index 0000000000..5a6b9d1ed9 --- /dev/null +++ b/src/amuse/community/uclchem/interface.py @@ -0,0 +1,356 @@ +""" +Interface for the UCLChem code +""" + +from amuse.units import units + +# from amuse.community.interface.chem import ChemicalModelingInterface +# from amuse.community.interface.chem import ChemicalModeling +from amuse.community.interface.common import CommonCode, CommonCodeInterface +from amuse.community import ( + CodeInterface, + LiteratureReferencesMixIn, + legacy_function, + LegacyFunctionSpecification, + InCodeComponentImplementation, +) +from pathlib import Path + + +class UCLchemInterface(CodeInterface, CommonCodeInterface, LiteratureReferencesMixIn): + """ + UCLCHEM: A Gas-Grain Chemical Code for astrochemical modelling + + .. [#] ADS:2017AJ....154...38H (Holdship, J. ; Viti, S, ; Jiménez-Serra, I.; Makrymallis, A. ; Priestley, F. , 2017, AJ) + """ + + def __init__(self, mode="cpu", **options): + CodeInterface.__init__(self, name_of_the_worker="uclchem_worker", **options) + LiteratureReferencesMixIn.__init__(self) + # New units used in chemical simulation + # cr_ion = named("cosmic ray ionisation rate", "cr_ion", 1.3e-17 * units.s**-1) + # habing = named("habing", "hab", 1.6e-3 * units.erg * units.cm**-2 * units.s**-1) + + @legacy_function + def commit_particles(): + # Standard function + function = LegacyFunctionSpecification() + function.result_type = "i" + return function + + @legacy_function + def recommit_particles(): + # Standard function + function = LegacyFunctionSpecification() + function.result_type = "i" + return function + + @legacy_function + def commit_parameters(): + # Standard function; doesn't actually do anything, but is required + function = LegacyFunctionSpecification() + function.result_type = "i" + return function + + def recommit_parameters(): + # Standard function; for UCLchem it's just the same as commit_parameters() + return self.commit_parameters() + + @legacy_function + def set_state(): + # Standard function; needs the same parameters as new_particles() + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter("id", dtype="i", direction=function.IN) + for x in ["number_density", "temperature", "ionrate", "radfield"]: + function.addParameter(x, dtype="d", direction=function.IN) + function.result_type = "i" + return function + + @legacy_function + def get_state(): + # Standard function; needs the same parameters as new_particles() + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter("id", dtype="i", direction=function.IN) + for x in ["number_density", "temperature", "ionrate", "radfield"]: + function.addParameter(x, dtype="d", direction=function.OUT) + function.result_type = "i" + return function + + @legacy_function + def get_abundance(): + # Standard function + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter("id", dtype="i", direction=function.IN) + function.addParameter("aid", dtype="i", direction=function.IN) + function.addParameter("abundance", dtype="d", direction=function.OUT) + function.result_type = "i" + return function + + @legacy_function + def set_abundance(): + # Standard function + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter("id", dtype="i", direction=function.IN) + function.addParameter("aid", dtype="i", direction=function.IN) + function.addParameter("abundance", dtype="d", direction=function.IN) + function.result_type = "i" + return function + + @legacy_function + def get_firstlast_abundance(): + # Standard function; mostly used internally, not really useful to most users + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter("first", dtype="i", direction=function.OUT) + function.addParameter("last", dtype="i", direction=function.OUT) + function.result_type = "i" + return function + + @legacy_function + def get_name_of_species(): + # Standard function + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter("index", dtype="i", direction=function.IN) + function.addParameter("name", dtype="s", direction=function.OUT) + function.result_type = "i" + return function + + @legacy_function + def get_index_of_species(): + # Standard function + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter("name", dtype="s", direction=function.IN) + function.addParameter("index", dtype="i", direction=function.OUT) + function.result_type = "i" + return function + + @legacy_function + def new_particle(): + # Standard function; sets required parameters for the particle set + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter("id", dtype="int32", direction=function.OUT) + for x in ["number_density", "temperature", "ionrate", "radfield"]: + function.addParameter(x, dtype="d", direction=function.IN) + function.result_type = "i" + return function + + @legacy_function + def delete_particle(): + # Standard function + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter("id", dtype="int32", direction=function.IN) + function.result_type = "i" + return function + + @legacy_function + def get_time(): + # Standard function + function = LegacyFunctionSpecification() + function.addParameter("time", dtype="float64", direction=function.OUT) + function.result_type = "i" + return function + + @legacy_function + def run_model(): + # Implementation here is UCLchem-specific; evolve_model is the function to call + function = LegacyFunctionSpecification() + function.addParameter("dictionary", dtype="s", direction=function.IN) + function.result_type = "i" + return function + + @legacy_function + def set_species(): + # Standard function + function = LegacyFunctionSpecification() + function.addParameter("species", dtype="s", direction=function.IN) + function.result_type = "i" + return function + + def _reform_inputs(self, param_dict, out_species): + # UCLchem-specific function; reforms the dictionary input to a string that is parsed in the FORTRAN layer + # Copied verbatim from the UCLchem source code + if param_dict is None: + param_dict = {} + else: + # lower case (and conveniently copy so we don't edit) the user's dictionary + # this is key to UCLCHEM's "case insensitivity" + new_param_dict = {} + for k, v in param_dict.items(): + assert ( + k.lower() not in new_param_dict + ), f"Lower case key {k} is already in the dict, stopping" + if isinstance(v, Path): + v = str(v) + new_param_dict[k.lower()] = v + param_dict = new_param_dict.copy() + del new_param_dict + if out_species is not None: + n_out = len(out_species) + param_dict["outspecies"] = n_out + out_species = " ".join(out_species) + else: + out_species = "" + n_out = 0 + return n_out, param_dict, out_species + + +class UCLchem(CommonCode): + def __init__(self, convert_nbody=None, **options): + legacy_interface = UCLchemInterface(**options) + self.uclchem_time = 0.0 | units.yr + InCodeComponentImplementation.__init__(self, legacy_interface) + # New units used in chemical simulation + # cr_ion = named("cosmic ray ionisation rate", "cr_ion", 1.3e-17 * units.s**-1) + # habing = named("habing", "hab", 1.6e-3 * units.erg * units.cm**-2 * units.s**-1) + + def evolve_model(self, tend): + # Standard function; implementation is UCLchem-specific + assert tend >= self.uclchem_time, "end time must be larger than uclchem_time" + dictionary, out_species = self._build_dict(tend=tend) + self.set_species(out_species) + output = self.run_model(dictionary) + self.uclchem_time = tend + return output + + def _build_dict(self, tend): + # UCLchem specific function; gets all relevant variables and creates a dictionary from them + dictionary_list = [] + outSpecies = self.out_species + attributes = self.particles.get_attribute_names_defined_in_store() + for i in range(len(self.particles.key)): + dictionary = dict() + if "temperature" in attributes: + dictionary["initialTemp"] = self.particles.temperature.value_in( + units.K + )[i] + if "number_density" in attributes: + dictionary["initialDens"] = self.particles.number_density.value_in( + units.cm**-3 + )[i] + if "ionrate" in attributes: + dictionary["zeta"] = self.particles.ionrate.value_in(units.cr_ion)[i] + if "radfield" in attributes: + dictionary["radfield"] = self.particles.radfield.value_in(units.habing)[ + i + ] + dictionary["finalTime"] = tend.value_in( + units.yr + ) - self.uclchem_time.value_in(units.yr) + _, dictionary, outSpecies_out = self._reform_inputs(dictionary, outSpecies) + dictionary_list.append(str(dictionary)) + return dictionary_list, str(outSpecies_out) + + def define_parameters(self, handler): + # UCLchem-specific parameter + handler.add_interface_parameter( + "out_species", "Array of molecules to use", default_value=["H", "H2"] + ) + + def define_methods(self, handler): + CommonCode.define_methods(self, handler) + handler.add_method( + "set_state", + ( + handler.INDEX, + units.cm**-3, + units.K, + units.s**-1, + units.habing, + ), + (handler.ERROR_CODE,), + ) + handler.add_method( + "get_state", + (handler.INDEX,), + ( + units.cm**-3, + units.K, + units.s**-1, + units.habing, + handler.ERROR_CODE, + ), + ) + handler.add_method( + "get_abundance", + ( + handler.INDEX, + handler.INDEX, + ), + ( + handler.NO_UNIT, + handler.ERROR_CODE, + ), + ) + handler.add_method( + "set_abundance", + ( + handler.INDEX, + handler.INDEX, + handler.NO_UNIT, + ), + (handler.ERROR_CODE,), + ) + handler.add_method( + "get_firstlast_abundance", + (), + ( + handler.NO_UNIT, + handler.NO_UNIT, + handler.ERROR_CODE, + ), + ) + handler.add_method( + "new_particle", + (units.cm**-3, units.K, units.s**-1, units.habing), + ( + handler.INDEX, + handler.ERROR_CODE, + ), + ) + handler.add_method("delete_particle", (handler.INDEX,), (handler.ERROR_CODE,)) + handler.add_method("evolve_model", (units.yr), (handler.ERROR_CODE)) + handler.add_method("get_time", (), (units.yr, handler.ERROR_CODE)) + + def define_particle_sets(self, handler): + handler.define_set("particles", "id") + 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_gridded_getter( + "particles", + "get_abundance", + "get_firstlast_abundance", + names=("abundances",), + ) + handler.add_gridded_setter( + "particles", + "set_abundance", + "get_firstlast_abundance", + names=("abundances",), + ) + + def define_state(self, handler): + CommonCode.define_state(self, handler) + handler.add_transition("INITIALIZED", "EDIT", "commit_parameters") + handler.add_transition("RUN", "PARAMETER_CHANGE_A", "invoke_state_change2") + handler.add_transition("EDIT", "PARAMETER_CHANGE_B", "invoke_state_change2") + handler.add_transition("PARAMETER_CHANGE_A", "RUN", "recommit_parameters") + handler.add_transition("PARAMETER_CHANGE_B", "EDIT", "recommit_parameters") + handler.add_method("EDIT", "new_particle") + handler.add_method("EDIT", "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_method("RUN", "evolve_model") + handler.add_method("RUN", "get_state") + handler.add_method("RUN", "get_abundance") diff --git a/src/amuse/community/uclchem/src/Makefile b/src/amuse/community/uclchem/src/Makefile new file mode 100755 index 0000000000..bf2990ed00 --- /dev/null +++ b/src/amuse/community/uclchem/src/Makefile @@ -0,0 +1,53 @@ +SHELL = /bin/sh +################################################################################ +# +# User Options +################################################################################ +FC =gfortran +#FC=ifort + +#f2pyFC=gnu95 #this must match FC (ifort=intelem, gfortran=gnu95) + +#Unforgiving debugging flags +#FFLAGS =-g -fbacktrace -Wall -fcheck=all -fPIC +#Fast optimizing flags +FFLAGS = -O3 -fPIC -ffree-line-length-0 -g + + +############################################################################### +EXT_SUFFIX = $(shell python3-config --extension-suffix) + + +PHYSICS_OBS= physics-core.o cloud.o hotcore.o sputtering.o cshock.o jshock.o collapse.o +CHEMISTRY_OBS = photoreactions.o surfacereactions.o chemistry.o +##simple makefile for uclchem +##user must point the FC variable to their preferred fortran compiler +##builds ode solver, physics module and chemistry module before linking together for main + +##User can also compile a python module which contains the subroutines in wrap.f90 as functions +##to do this run "python -m numpy.f2py -c --help-fcompiler" to find your fortran compiler +##and edit the python makef + +##physics module selected by changing physics variable to chosen fortran file. +main: constants.o network.o ${PHYSICS_OBS} dvode.o ${CHEMISTRY_OBS} io.o wrap.f90 defaultparameters.f90 + ${FC} ${FFLAGS} -o ../../uclchem io.o ${PHYSICS_OBS} dvode.o ${CHEMISTRY_OBS} network.o wrap.f90 main.f90 + +##module for compiling in amuse +libchem.a: constants.o network.o ${PHYSICS_OBS} dvode.o ${CHEMISTRY_OBS} io.o hash.o amuse_helpers.o defaultparameters.f90 + ${AR} crs libchem.a io.o ${PHYSICS_OBS} dvode.o ${CHEMISTRY_OBS} network.o hash.o amuse_helpers.o + +# The order here is the compilation order and it's vital that modules are compiled after their dependencies +python: constants.o network.o ${PHYSICS_OBS} dvode.o ${CHEMISTRY_OBS} io.o wrap.f90 defaultparameters.f90 + python3 -m numpy.f2py -c --fcompiler=${f2pyFC} io.o ${PHYSICS_OBS} dvode.o ${CHEMISTRY_OBS} network.o -m uclchemwrap wrap.f90 + +chemistry.o: odes.f90 rates.f90 chemistry.f90 + ${FC} ${FFLAGS} -c chemistry.f90 + +sputtering.o: surfacereactions.o sputtering.f90 + ${FC} ${FFLAGS} -c sputtering.f90 + +%.o: %.f90 + ${FC} ${FFLAGS} -c $< + +clean: + rm -rf *.o *.mod *.so ../uclchem/*.so *.a \ No newline at end of file diff --git a/src/amuse/community/uclchem/src/amuse_helpers.f90 b/src/amuse/community/uclchem/src/amuse_helpers.f90 new file mode 100644 index 0000000000..4527e0a336 --- /dev/null +++ b/src/amuse/community/uclchem/src/amuse_helpers.f90 @@ -0,0 +1,685 @@ +MODULE uclchemhelper + USE physicscore + USE chemistry + USE io + USE constants + use network + IMPLICIT NONE + type particle_type + !Particles are defined here + integer :: id + double precision :: density !unit = cm^-3 + double precision :: temperature !unit = K + double precision :: ionrate !unit = galactic cr ionization rate + double precision :: uvrad !unit = Habing + double precision :: abundances(nSpec+1) + end type + integer :: nmols + type(particle_type), allocatable :: particles(:) + + double precision :: tcurrent ! time unit = yr + + integer :: nparticle + integer :: tot_id + + character(len=500) :: out_species + + integer, parameter :: NMAX=1000000 + + logical :: particles_searcheable=.FALSE. + integer, allocatable :: pid(:) + +CONTAINS + function chem_initialize() result(ret) + integer :: ret + tcurrent=0. + nparticle=0 + tot_id=0 + nmols=nspec + if(.not.allocated(particles)) allocate(particles(NMAX)) + particles(:)%density=0. + + ret=0 + end function + + function chem_commit_parameters() result(ret) + integer :: ret + ret=0 + end function + + function chem_commit_particles() result (ret) + integer :: ret + integer :: n + integer :: i + + particles_searcheable=.FALSE. + nparticle=clean_particles(particles) + ret=0 + end function + + function chem_end() result(ret) + integer :: ret + if(allocated(particles)) deallocate(particles) + ret=0 + end function + + function set_particle_state(id,density,temperature,ionrate,uvrad) result(ret) + integer :: ret + integer :: id,index + double precision :: density, temperature, ionrate, uvrad + index=find_particle(id) + if(index.LT.0) then + ret=index + return + endif + + particles(index)%density=density + particles(index)%temperature=temperature + particles(index)%ionrate=ionrate + particles(index)%uvrad=uvrad + ret=0 + + end function + + function get_particle_state(id,density,temperature,ionrate,uvrad) result(ret) + integer :: ret + integer :: id,index + double precision :: density, temperature, ionrate, uvrad + index=find_particle(id) + if(index.LT.0) then + ret=index + return + endif + + density=particles(index)%density + temperature=particles(index)%temperature + ionrate=particles(index)%ionrate + uvrad=particles(index)%uvrad + ret=0 + + end function + + function get_particle_abundance(id, aid, abundance) result(ret) + integer :: ret + integer :: id,index,aid + double precision :: abundance + index=find_particle(id) + if(index.LT.0) then + ret=index + return + endif + if(aid.LT.1.OR.aid.GT.500) then + ret=-1 + return + endif + abundance=particles(index)%abundances(aid) + ret=0 + end function + + function set_particle_abundance(id, aid, abundance) result(ret) + integer :: ret + integer :: id,index,aid + double precision :: abundance + index=find_particle(id) + if(index.LT.0) then + ret=index + return + endif + if(aid.LT.1.OR.aid.GT.500) then + ret=-1 + return + endif + particles(index)%abundances(aid)=abundance + ret=0 + end function + + function add_particle(id,density,temperature,ionrate,uvrad) result(ret) + use network + integer :: ret + integer :: i,id + double precision :: density, temperature, ionrate, uvrad + double precision :: x(500) + particles_searcheable=.FALSE. + id=new_id() + i=nparticle+1 + + if(i.GT.NMAX) then + ret=-1 + return + endif + particles(i)%id=id + particles(i)%density=density + particles(i)%temperature=temperature + particles(i)%ionrate=ionrate + particles(i)%uvrad=uvrad + particles(i)%abundances=1.d-40 + + if(density.GT.0) then + !values taken from the UCLchem default parameters + particles(i)%abundances(nh) = 0.5 !H + particles(i)%abundances(nh2) = 0.25 !H2 + particles(i)%abundances(nhe) = 0.1 !He + particles(i)%abundances(ncx) = 1.77d-04 !C+ (C is fully ionised by default) + particles(i)%abundances(nc) = 1.d-10 !C + particles(i)%abundances(no) = 3.34d-04 !O + particles(i)%abundances(nn) = 6.18d-05 !N + particles(i)%abundances(nmg) = 2.256d-06 !Mg + particles(i)%abundances(np) = 7.78d-08 !P + particles(i)%abundances(nf) = 3.6d-08 ! + particles(i)%abundances(nsx) = 3.51d-6 !S + particles(i)%abundances(nsix) = 1.78d-06 !Si + particles(i)%abundances(nclx) = 3.39d-08 !Cl + particles(i)%abundances(nelec) = 1.77d-04 + 3.51d-6 + 1.78d-06 + 3.39d-08 !electrons; sum of ions + endif + + nparticle=nparticle+1 + ret=0 + end function + + function remove_particle(id) result(ret) + integer :: ret + integer :: i,id + i=find_particle(id) + if(i.LE.0) then + ret=i + return + endif + if(particles(i)%density.LT.0) then + ret=-4 + return + endif + particles(i)%density=-1. + ret=0 + end function + + function new_id() + integer new_id + tot_id=tot_id+1 + new_id=tot_id + end function + + function find_particle(id_) result(index) + !Function to find a particle with a given id. Taken from the Krome interface. + use hashMod + type(hash_type),save :: hash + integer id_,index + integer, save :: nbod=0 + + if(.NOT.particles_searcheable) then + nbod=nparticle + if(allocated(pid)) deallocate(pid) + allocate(pid(nbod)) + pid(1:nbod)=particles(1:nbod)%id + call initHash(nbod/2+1,nbod, pid,hash) + particles_searcheable=.TRUE. + endif + + + index=find(id_,pid,hash) + + if(index.LE.0) then + index=-1 + return + endif + if(index.GT.nbod) then + index=-2 + return + endif + if(pid(index).NE.id_) then + index=-3 + return + endif + + end function + + function get_current_time(time) result(ret) + integer :: ret + double precision :: time + time = tcurrent + ret = 0 + end function + + function clean_particles(par) result(np) + integer :: left,right,np + type(particle_type), allocatable :: par(:) + type(particle_type) :: tmp + left=1 + if(.NOT.allocated(par)) then + np = 0 + return + endif + right=size(par) + if(right.EQ.0) then + np=0 + return + endif + do while(.TRUE.) + do while(par(left)%density.GT.0.AND.left.LT.right) + left=left+1 + enddo + do while(par(right)%density.LE.0.AND.left.LT.right) + right=right-1 + enddo + if(left.LT.right) then + tmp=par(left) + par(left)=par(right) + par(right)=tmp + else + exit + endif + enddo + if(par(left)%density.GT.0) left=left+1 + np=left-1 + end function + + function simple_evolution(dictionary, outSpeciesIn) result(ret) + !Evolves each particle seperately + integer :: ret + integer :: i, iret + CHARACTER(LEN=*) :: dictionary(:) + CHARACTER(LEN=*) :: outSpeciesIn + + ret = 0 + INCLUDE 'defaultparameters.f90' + print *, outSpeciesIn + do i=1,nparticle + CALL dictionaryParser(dictionary(i), outSpeciesIn,ret) + iret = evolve_1_particle(particles(i)) + ret = min(iret,ret) + enddo + tcurrent = timeInYears + return + + end function + + function evolve_1_particle(part) result(ret) + use cloud_mod + type(particle_type) :: part + integer :: ret + + dstep=1 + + + call coreInitializePhysics(ret) + CALL coreInitializePhysics(ret) + + CALL initializeChemistry(readabunds=.FALSE.) + + dstep=1 + !Make sure the saved abundances from previous steps are used + abund(:,1) = part%abundances + DO WHILE (((endAtFinalDensity) .and. (density(1) < finalDens)) .or. & + &((.not. endAtFinalDensity) .and. (timeInYears < finalTime))) + currentTimeold=currentTime + !In standard UCLchem, each physics module has its own timestep scheme. + CALL updateTargetTime + !loop over parcels, counting from centre out to edge of cloud + DO dstep=1,points + !reset time if this isn't first depth point + + currentTime=currentTimeold + !update chemistry from currentTime to targetTime + CALL updateChemistry(ret) + IF (ret .lt. 0) THEN + write(*,*) 'Error updating chemistry' + RETURN + END IF + + !get time in years for output, currentTime is now equal to targetTime + timeInYears= currentTime/SECONDS_PER_YEAR + + !Update physics so it's correct for new currentTime and start of next time step + Call coreUpdatePhysics + !Sublimation checks if Sublimation should happen this time step and does it + CALL sublimation(abund) + + !write this depth step now time, chemistry and physics are consistent + CALL output + END DO + END DO + !Save the computed abundances to the particle set + part%abundances=abund(:,1) + currentTime = 0.0 + end function + + SUBROUTINE get_rates(dictionary,abundancesIn,speciesIndx,rateIndxs,& + &speciesRates,successFlag,transfer,swap,bulk_layers) + !Given a species of interest, some parameters and abundances, this subroutine + !return the rate of all reactions that include that species plus some extra variables + !to allow for the calculation of the rate of bulk/surface ice transfer. + USE cloud_mod + USE network, only : nspec + CHARACTER(LEN=*):: dictionary + DOUBLE PRECISION :: abundancesIn(500),speciesRates(500) + DOUBLE PRECISION :: transfer,swap,bulk_layers + INTEGER:: rateIndxs(500),speciesIndx, successFlag + DOUBLE PRECISION :: ydot(nspec+1) + INTEGER :: speci,bulk_version,surface_version + !f2py intent(in) dictionary,abundancesIn,speciesIndx,rateIndxs + !f2py intent(out) speciesRates,successFlag,transfer,swap,bulk_layers + INCLUDE 'defaultparameters.f90' + + CALL dictionaryParser(dictionary, "",successFlag) + IF (successFlag .lt. 0) THEN + WRITE(*,*) 'Error reading parameter dictionary' + RETURN + END IF + CALL coreInitializePhysics(successFlag) + CALL initializePhysics(successFlag) + IF (successFlag .lt. 0) then + WRITE(*,*) 'Error initializing physics' + RETURN + END IF + + CALL initializeChemistry(readAbunds) + dstep=1 + successFlag=1 + abund(:nspec,dstep)=abundancesIn(:nspec) + abund(neq,dstep)=initialDens + currentTime=0.0 + timeInYears=0.0 + + targetTime=1.0d-7 + CALL updateChemistry(successFlag) + + CALL F(NEQ,currentTime,abund(:,dstep),ydot) + + speciesRates=rate(rateIndxs) + + IF ((specname(speciesIndx)(1:1) .eq. "#") .or.& + & (specname(speciesIndx)(1:1) .eq. "@")) THEN + DO speci=1,nSpec + IF (specname(speci) .eq. "@"//specname(speciesIndx)(2:)) bulk_version=speci + IF (specname(speci) .eq. "#"//specname(speciesIndx)(2:)) surface_version=speci + END DO + IF (YDOT(nsurface) .lt. 0) THEN + transfer=YDOT(nsurface)*surfaceCoverage*abund(bulk_version,1)/safeBulk + ELSE + transfer=YDOT(nsurface)*surfaceCoverage*abund(surface_version,1) + END If + swap=totalSwap + bulk_layers=bulkLayersReciprocal + ELSE + swap=0.0 + transfer=0.0 + bulk_layers=0.0 + END IF + + END SUBROUTINE get_rates + + SUBROUTINE dictionaryParser(dictionary, outSpeciesIn,successFlag) + !Reads the input parameters from a string containing a python dictionary/JSON format + !set of parameter names and values. + !dictionary - lowercase keys matching the names of the parameters in the select case below + !OutSpeciesIn - string containing the species to output + !successFlag - integer flag to indicate success + CHARACTER(LEN=*) :: dictionary, outSpeciesIn + INTEGER, INTENT(OUT) :: successFlag + INTEGER, ALLOCATABLE, DIMENSION(:) :: locations + LOGICAL :: ChemicalDuplicateCheck + INTEGER :: posStart, posEnd, whileInteger,inputindx + CHARACTER(LEN=100) :: inputParameter, inputValue + close(10) + close(11) + close(7) + + !always deallocate these so that if user didn't specify them, + ! they don't remain from previous run + IF (ALLOCATED(outIndx)) DEALLOCATE(outIndx) + IF (ALLOCATED(outSpecies)) DEALLOCATE(outSpecies) + + !All reads use IOSTAT which will change successFlag from 0 if an error occurs + !so set zero and check for non-zero each loop. + successFlag=0 + + whileInteger = 0 + + posStart = scan(dictionary, '{') + DO WHILE (whileInteger .NE. 1) + posEnd = scan(dictionary, ':') + inputParameter = dictionary(posStart+2:posEnd-2) + dictionary = dictionary(posEnd:) + posStart = scan(dictionary, ' ') + IF (scan(dictionary, ',') .EQ. 0) THEN + posEnd = scan(dictionary, '}') + whileInteger = 1 + ELSE + posEnd = scan(dictionary, ',') + END IF + inputValue = dictionary(posStart+1:posEnd-1) + + SELECT CASE (inputParameter) + CASE('alpha') + !To provide alphas, set keyword alpha in inputdictionary with a dictionary value + !that dictionary should be index:value pairs for the alpha array + posStart=scan(dictionary,'{') + posEnd=scan(dictionary,'}') + CALL coefficientParser(dictionary(posStart+1:posEnd),alpha) + CASE('beta') + !To provide alphas, set keyword alpha in inputdictionary with a dictionary value + !that dictionary should be index:value pairs for the alpha array + posStart=scan(dictionary,'{') + posEnd=scan(dictionary,'}') + CALL coefficientParser(dictionary(posStart+1:posEnd),beta) + CASE('gamma') + !To provide alphas, set keyword alpha in inputdictionary with a dictionary value + !that dictionary should be index:value pairs for the alpha array + posStart=scan(dictionary,'{') + posEnd=scan(dictionary,'}') + CALL coefficientParser(dictionary(posStart+1:posEnd),gama) + CASE('initialtemp') + READ(inputValue,*,iostat=successFlag) initialTemp + CASE('initialdens') + READ(inputValue,*,iostat=successFlag) initialDens + CASE('finaldens') + READ(inputValue,*,iostat=successFlag) finalDens + CASE('currenttime') + READ(inputValue,*,iostat=successFlag) currentTime + CASE('finaltime') + READ(inputValue,*,iostat=successFlag) finalTime + CASE('radfield') + READ(inputValue,*,iostat=successFlag) radfield + CASE('zeta') + READ(inputValue,*,iostat=successFlag) zeta + CASE('freezefactor') + READ(inputValue,*,iostat=successFlag) freezeFactor + CASE('rout') + READ(inputValue,*,iostat=successFlag) rout + CASE('rin') + READ(inputValue,*,iostat=successFlag) rin + CASE('baseav') + READ(inputValue,*,iostat=successFlag) baseAv + CASE('points') + READ(inputValue,*,iostat=successFlag) points + CASE('endatfinaldensity') + Read(inputValue,*,iostat=successFlag) endAtFinalDensity + CASE('freefall') + READ(inputValue,*,iostat=successFlag) freefall + CASE('freefallfactor') + READ(inputValue,*,iostat=successFlag) freefallFactor + CASE('desorb') + READ(inputValue,*,iostat=successFlag) desorb + CASE('h2desorb') + READ(inputValue,*,iostat=successFlag) h2desorb + CASE('crdesorb') + READ(inputValue,*,iostat=successFlag) crdesorb + CASE('uvdesorb') + READ(inputValue,*,iostat=successFlag) uvdesorb + CASE('thermdesorb') + READ(inputValue,*,iostat=successFlag) uvdesorb + CASE('instantsublimation') + READ(inputValue,*,iostat=successFlag) instantSublimation + CASE('cosmicrayattenuation') + READ(inputValue,*,iostat=successFlag) cosmicRayAttenuation + CASE('ionmodel') + READ(inputValue,*,iostat=successFlag) ionModel + CASE('improvedh2crpdissociation') + READ(inputValue,*,iostat=successFlag) improvedH2CRPDissociation + CASE('ion') + READ(inputValue,*,iostat=successFlag) ion + CASE('fhe') + READ(inputValue,*,iostat=successFlag) fhe + CASE('fc') + READ(inputValue,*,iostat=successFlag) fc + CASE('fo') + READ(inputValue,*,iostat=successFlag) fo + CASE('fn') + READ(inputValue,*,iostat=successFlag) fn + CASE('fs') + READ(inputValue,*,iostat=successFlag) fs + CASE('fmg') + READ(inputValue,*,iostat=successFlag) fmg + CASE('fsi') + READ(inputValue,*,iostat=successFlag) fsi + CASE('fcl') + READ(inputValue,*,iostat=successFlag) fcl + CASE('fp') + READ(inputValue,*,iostat=successFlag) fp + CASE('ff') + READ(inputValue,*,iostat=successFlag) ff + CASE('outspecies') + READ(inputValue,*,iostat=successFlag) nout + ALLOCATE(outIndx(nout)) + ALLOCATE(outSpecies(nout)) + IF (outSpeciesIn .eq. "") THEN + write(*,*) "Outspecies parameter set but no outspecies string given" + write(*,*) "general(parameter_dict,outSpeciesIn) requires a delimited string of species names" + write(*,*) "if outSpecies or columnFlag is set in the parameter dictionary" + successFlag=-1 + RETURN + ELSE + READ(outSpeciesIn,*, END=22) outSpecies + IF (outSpeciesIn .eq. "") THEN +22 write(*,*) "mismatch between outSpeciesIn and number given in dictionary" + write(*,*) "Number:",nout + write(*,*) "Species list:",outSpeciesIn + successFlag=-1 + RETURN + END IF + END IF + !assign array indices for important species to the integers used to store them. + DO i=1,nspec + DO j=1,nout + IF (specname(i).eq.outSpecies(j)) outIndx(j)=i + END DO + END DO + CASE('writestep') + READ(inputValue,*,iostat=successFlag) writeStep + CASE('ebmaxh2') + READ(inputValue,*,iostat=successFlag) ebmaxh2 + CASE('epsilon') + READ(inputValue,*,iostat=successFlag) epsilon + CASE('uvcreff') + READ(inputValue,*,iostat=successFlag) uvcreff + CASE('ebmaxcr') + READ(inputValue,*,iostat=successFlag) ebmaxcr + CASE('phi') + READ(inputValue,*,iostat=successFlag) phi + CASE('ebmaxuvcr') + READ(inputValue,*,iostat=successFlag) ebmaxuvcr + CASE('uv_yield') + READ(inputValue,*,iostat=successFlag) uv_yield + CASE('metallicity') + READ(inputValue,*,iostat=successFlag) metallicity + CASE('omega') + READ(inputValue,*,iostat=successFlag) omega + CASE('reltol') + READ(inputValue,*,iostat=successFlag) reltol + CASE('abstol_factor') + READ(inputValue,*,iostat=successFlag) abstol_factor + CASE('abstol_min') + READ(inputValue,*,iostat=successFlag) abstol_min + ! CASE('jacobian') + ! READ(inputValue,*) jacobian + CASE('abundsavefile') + READ(inputValue,*,iostat=successFlag) abundSaveFile + abundSaveFile = TRIM(abundSaveFile) + open(abundSaveID,file=abundSaveFile,status="unknown") + CASE('abundloadfile') + READ(inputValue,*,iostat=successFlag) abundLoadFile + abundLoadFile = TRIM(abundLoadFile) + open(abundLoadID,file=abundLoadFile,status='old') + CASE('outputfile') + READ(inputValue,*,iostat=successFlag) outFile + outputFile = trim(outFile) + fullOutput=.True. + open(outputId,file=outputFile,status='unknown',iostat=successFlag) + IF (successFlag .ne. 0) THEN + write(*,*) "An error occured when opening the output file!"//& + & NEW_LINE('A')//& + &" The failed file was ",outputFile& + &, NEW_LINE('A')//"A common error is that the directory doesn't exist"& + &//NEW_LINE('A')//"************************" + successFlag=-1 + RETURN + END IF + CASE('columnfile') + IF (trim(outSpeciesIn) .NE. '' ) THEN + columnOutput=.True. + READ(inputValue,*,iostat=successFlag) columnFile + columnFile = trim(columnFile) + open(columnId,file=columnFile,status='unknown') + ELSE + WRITE(*,*) "Error in output species. No species were given but a column file was given." + WRITE(*,*) "columnated output requires output species to be chosen." + successFlag=-1 + RETURN + END IF + + CASE DEFAULT + WRITE(*,*) "Problem with given parameter: '", trim(inputParameter), "'." + WRITE(*,*) "This is either not supported yet, or invalid." + successFlag=-1 + RETURN + END SELECT + dictionary = dictionary(posEnd:) + IF (SCAN(dictionary,',') .eq. 0) whileInteger=1 + + !check for failure + IF (successFlag .ne. 0) THEN + WRITE(*,*) "Error reading ",inputParameter + write(*,*) "This is usually due to wrong type." + successFlag=PARAMETER_READ_ERROR + RETURN + END IF + END DO + + END SUBROUTINE dictionaryParser + + SUBROUTINE coefficientParser(coeffDictString,coeffArray) + !Similar to dictionaryParser, it reads a python dictionary + !however, it's intended to read pairs of reaction indices and coefficient values + !for the alpha, beta, and gama arrays. + ! No return value, just modifies the coeffArray + CHARACTER(LEN=*) :: coeffDictString + REAL(dp), INTENT(INOUT) :: coeffArray(*) + INTEGER :: inputIndx,posStart,posEnd + CHARACTER(LEN=100) :: inputValue + LOGICAL :: continue_flag + + continue_flag=.True. + DO WHILE (continue_flag) + !substring containing integer key + posStart=1 + posEnd=SCAN(coeffDictString,':') + !read it into index integer + READ(coeffDictString(posStart:posEnd-1),*) inputindx + + !substring including alpha value for the index. + posStart=posEnd+1 + posEnd=SCAN(coeffDictString,',') + !last value will have a } instead of , so grab index and tell loop to finish + IF (posEnd .eq. 0) THEN + posEnd=SCAN(coeffDictString,"}") + continue_flag=.False. + END IF + + !read that substring + inputValue=coeffDictString(posStart:posEnd-1) + READ(inputValue,*) coeffArray(inputIndx) + !update string to remove this entry + coeffDictString=coeffDictString(posEnd+1:) + END DO + END SUBROUTINE coefficientParser + +END MODULE uclchemhelper \ No newline at end of file diff --git a/src/amuse/community/uclchem/src/chemistry.f90 b/src/amuse/community/uclchem/src/chemistry.f90 new file mode 100755 index 0000000000..844a687eed --- /dev/null +++ b/src/amuse/community/uclchem/src/chemistry.f90 @@ -0,0 +1,318 @@ +! Chemistry module of UCL_CHEM. ! +! Contains all the core machinery of the code, not really intended to be altered in standard ! +! use. Use a (custom) physics module to alter temp/density behaviour etc. ! +! ! +! chemistry module contains rates.f90, a series of subroutines to calculate all reaction rates! +! when updateChemistry is called from main, these rates are calculated, the ODEs are solved ! +! from currentTime to targetTime to get abundances at targetTime and then all abundances are ! +! written to the fullOutput file. ! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +MODULE chemistry +USE physicscore +USE dvode_f90_m +USE network +USE photoreactions +USE surfacereactions +USE constants +IMPLICIT NONE + !These integers store the array index of important species and reactions, x is for ions + !loop counters + INTEGER :: i,j,l,writeStep,writeCounter=0,loopCounter,failedIntegrationCounter + INTEGER, PARAMETER :: maxLoops=10,maxConsecutiveFailures=10 + + !Flags to control desorption processes + LOGICAL :: desorb,h2desorb,crdesorb,uvdesorb,thermdesorb + + + !Array to store reaction rates + REAL(dp) :: rate(nreac) + + + !DLSODE variables + INTEGER :: ITASK,ISTATE,NEQ,MXSTEP + REAL(dp) :: reltol,abstol_factor,abstol_min + REAL(dp), ALLOCATABLE :: abstol(:) + TYPE(VODE_OPTS) :: OPTIONS + !initial fractional elemental abudances and arrays to store abundances + REAL(dp) :: fh,fd,fhe,fc,fo,fn,fs,fmg,fsi,fcl,fp,ff,ffe,fli,fna,fpah,f15n,f13c,f18O,metallicity + REAL(dp) :: h2col,cocol,ccol,h2colToCell,cocolToCell,ccolToCell + REAL(dp),ALLOCATABLE :: abund(:,:) + + !Variables controlling chemistry + LOGICAL :: PARAMETERIZE_H2FORM=.True. + REAL(dp) :: radfield,freezeFactor,omega,grainArea,cion,h2dis,lastTemp=0.0 + REAL(dp) :: ebmaxh2,epsilon,ebmaxcr,phi,ebmaxuvcr,uv_yield,uvcreff + REAL(dp), PARAMETER :: h2StickingZero=0.87d0,hStickingZero=1.0d0, h2StickingTemp=87.0d0,hStickingTemp=52.0d0 + + + REAL(dp) :: turbVel=1.0 + REAL(dp) :: MIN_ABUND = 1.0d-30 !Minimum abundance allowed +CONTAINS + SUBROUTINE initializeChemistry(readAbunds) + LOGICAL, INTENT(IN) :: readAbunds + ! Sets variables at the start of every run. + ! Since python module persists, it's not enough to set initial + ! values in module definitions above. Reset here. + NEQ=nspec+1 + IF (ALLOCATED(abund)) DEALLOCATE(abund,vdiff) + ALLOCATE(abund(NEQ,points),vdiff(SIZE(iceList))) + !Set abundances to initial elemental if not reading them in. + IF (.NOT. readAbunds) THEN + !ensure abund is initially zero + abund= MIN_ABUND + + !Start by filling all metallicity scaling elements + !neutral atoms + abund(no,:) = fo + abund(nn,:) = fn + abund(nmg,:) = fmg + abund(np,:) = fp + abund(nf,:) = ff + !abund(nfe,:) = ffe + abund(nna,:) = fna + abund(nli,:) = fli + abund(npah,:) = fpah + !default to ions + abund(nsx,:) = fs + abund(nsix,:) = fsi + abund(nclx,:) = fcl + !Decide how much carbon is initiall ionized using parameters.f90 + SELECT CASE (ion) + CASE(0) + abund(nc,:)=fc + abund(ncx,:)=1.d-10 + CASE(1) + abund(nc,:)=fc*0.5 + abund(ncx,:)=fc*0.5 + CASE(2) + abund(nc,:)=1.d-10 + abund(ncx,:)=fc + END SELECT + + !isotopes + abund(n18o,:) = f18o + abund(n15n,:) = f15n + abund(n13c,:) = f13c + + abund(nelec,:)=abund(ncx,:)+abund(nsix,:)+abund(nsx,:)+abund(nclx,:) + + abund=abund*metallicity + + !Total H nuclei is always 1 so put fh into H and whatever is left over in H2 + abund(nh,:) = fh + abund(nh2,:) = 0.5*(1.0e0-fh) + abund(nd,:)=fd + + abund(nhe,:) = fhe + ENDIF + abund(neq,:)=density + !Initial calculations of diffusion frequency for each species bound to grain + !and other parameters required for diffusion reactions + DO i=lbound(iceList,1),ubound(iceList,1) + j=iceList(i) + vdiff(i)=VDIFF_PREFACTOR*bindingEnergy(i)/mass(j) + vdiff(i)=dsqrt(vdiff(i)) + END DO + + !DVODE SETTINGS + ISTATE=1 + ITASK=1 + + !set integration counts + loopCounter=0 + failedIntegrationCounter=0 + + IF (.NOT. ALLOCATED(abstol)) THEN + ALLOCATE(abstol(NEQ)) + END IF + !OPTIONS = SET_OPTS(METHOD_FLAG=22, ABSERR_VECTOR=abstol, RELERR=reltol,USER_SUPPLIED_JACOBIAN=.FALSE.) + + !Set rates to zero to ensure they don't hold previous values or random ones if we don't set them in calculateReactionRates + rate=0.0 + !We typically don't recalculate rates that only depend on temperature if the temp hasn't changed + !use arbitrarily high value to make sure they are calculated at least once. + lastTemp=99.0d99 + END SUBROUTINE initializeChemistry + + + + SUBROUTINE updateChemistry(successFlag) + !Updates the abundances for the next time step, first updating chemical variables and reaction rates, + !then by solving the ODE system to obtain new abundances. + !Solving ODEs is complex so we have two checks to try to automatically overcome difficulties and end stalled models + !Firstly, the integration subroutine is called up to maxLoops times whilst adjusting variables to help integration converge. + !If it succeeds before maxLoops, we continue as normal, otherwise we'll call it a fail. + !Secondly, we check for stalls caused by the the solver loop reducing the targetTime to overcome difficulties. + !That reduction the possibility of the code "succeeding" by integrating tiny target times. We have a counter that resets each time + !the code integrates to the planned targetTime rather than a reduced one. If the counter reaches maxConsecutiveFailures, we end the code. + + INTEGER, INTENT(OUT) :: successFlag + real(dp) :: originalTargetTime !targetTime can be altered by integrator but we'd like to know if it was changed + + !Integration can fail in a way that we can manage. Allow maxLoops tries before giving up. + loopCounter=0 + successFlag=1 + originalTargetTime=targetTime + DO WHILE((currentTime .lt. targetTime) .and. (loopCounter .lt. maxLoops)) + !allow option for dens to have been changed elsewhere. + IF (.not. freefall) abund(nspec+1,dstep)=density(dstep) + + !First sum the total column density over all points further towards edge of cloud + IF (dstep.gt.1) THEN + h2ColToCell=(sum(abund(nh2,:dstep-1)*density(:dstep-1)))*(cloudSize/real(points)) + coColToCell=(sum(abund(nco,:dstep-1)*density(:dstep-1)))*(cloudSize/real(points)) + cColToCell=(sum(abund(nc,:dstep-1)*density(:dstep-1)))*(cloudSize/real(points)) + ELSE + h2ColToCell=0.0 + coColToCell=0.0 + cColToCell=0.0 + ENDIF + !then add half the column density of the current point to get average in this "cell" + h2Col=h2ColToCell+0.5*abund(nh2,dstep)*density(dstep)*(cloudSize/real(points)) + coCol=coColToCell+0.5*abund(nco,dstep)*density(dstep)*(cloudSize/real(points)) + cCol=cColToCell+0.5*abund(nc,dstep)*density(dstep)*(cloudSize/real(points)) + + !Reset surface and bulk values in case of integration error or sputtering + abund(nBulk,dstep)=sum(abund(bulkList,dstep)) + abund(nSurface,dstep)=sum(abund(surfaceList,dstep)) + !recalculate coefficients for ice processes + safeMantle=MAX(1d-30,abund(nSurface,dstep)) + safeBulk=MAX(1d-30,abund(nBulk,dstep)) + + if (refractoryList(1) .gt. 0) safeBulk=safeBulk-SUM(abund(refractoryList,dstep)) + bulkLayersReciprocal=MIN(1.0,NUM_SITES_PER_GRAIN/(GAS_DUST_DENSITY_RATIO*safeBulk)) + surfaceCoverage=bulkGainFromMantleBuildUp() + + + CALL calculateReactionRates + + !Integrate chemistry, and return fail if unrecoverable error was reached + CALL integrateODESystem(successFlag) + IF (successFlag .lt. 0) THEN + write(*,*) "Integration failed, exiting" + RETURN + END IF + + + + !1.d-30 stops numbers getting too small for fortran. + WHERE(abund