UCLchem interface

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 <name>.code  - clean & build the community code <name> (or matching name*)"
+	@echo "make clean        - clean codes and libraries"
+	@echo "make distclean    - clean codes and libraries and all configuration files"

doc/interactive_tutorial/amuserc

@@ -1,2 +1,2 @@
-[output]
-printing_strategy=simple
+#[output]
+#printing_strategy=simple

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,16 +53,17 @@ 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))
     print("cloud n_H:", parts.number_density.max().in_(units.cm**-3))
     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()
 

src/amuse/amuserc

@@ -0,0 +1,2 @@
+[channel]
+# debugger=gdb

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
+
+
+

src/amuse/community/uclchem/__init__.py

@@ -0,0 +1 @@
+from .interface import UCLchem

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