gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into reconnection

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2020-05-12 11:35:29 -03:00
commit f2f19201aa
19 changed files with 5928 additions and 359 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

52
README.md
View File

@ -18,12 +18,15 @@ following features are already implemented:
reconstructions, reconstructions,
* Riemann solvers of Roe- and HLL-types (HLL, HLLC, and HLLD), * Riemann solvers of Roe- and HLL-types (HLL, HLLC, and HLLD),
* standard boundary conditions: periodic, open, reflective, hydrostatic, etc. * standard boundary conditions: periodic, open, reflective, hydrostatic, etc.
* turbulence driving using Alvelius or Ornstein–Uhlenbeck methods,
* viscous and resistive source terms, * viscous and resistive source terms,
* suppor for passive scalars (up to 100), * support for passive scalars (up to 100),
* data stored in the HDF5 format, * data stored in internal XML+binary or HDF5 format,
* Python interface to read snapshots in both formats,
* MPI parallelization, * MPI parallelization,
* completely written in Fortran 2003, * completely written in Fortran 2003,
* Python interface to read data. * minimum requirements, only Fortran compiler and Python are required to
prepare, run, and analyze your simulations.
This program is free software: you can redistribute it and/or modify it under This program is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software the terms of the GNU General Public License as published by the Free Software
@ -53,19 +56,21 @@ Requirements
version 18.10 or newer, version 18.10 or newer,
- [Intel Fortran](https://software.intel.com/en-us/fortran-compilers) - [Intel Fortran](https://software.intel.com/en-us/fortran-compilers)
compiler version 9.0 or newer. compiler version 9.0 or newer.
* [HDF5 libraries](https://www.hdfgroup.org/solutions/hdf5/), tested with * Optional, but recommended, [OpenMPI](https://www.open-mpi.org/) for parallel
version 1.8 or newer. runs, tested with version 1.8 or newer.
* [OpenMPI](https://www.open-mpi.org/) for parallel runs, tested with version * Optional [HDF5 libraries](https://www.hdfgroup.org/solutions/hdf5/), tested
1.8 or newer. with version 1.10 or newer. The code now uses the new XML-binary snapshot
format. However, if you still want to use older HDF5 snapshot format, you
will need these libraries.
Environment Variables Environment Variables
===================== =====================
If the HDF5 libraries are not installed in the default location, i.e. in the If you need to use the HDF5 libraries and they are not installed in the default
system directory **/usr**, make sure that the environment variable _HDF5DIR_ is location, i.e. in the system directory **/usr**, make sure that the environment
set in your **~/.bashrc** (or **~/.cshrc**) and pointing to the location where variable _HDF5DIR_ is set in your **~/.bashrc** (or **~/.cshrc**) and pointing
the HDF5 libraries have been installed. to the location where the HDF5 libraries have been installed.
Compilation Compilation
@ -104,3 +109,28 @@ In order to run parallel version (after compiling the code with MPI support),
type in your terminal: type in your terminal:
`mpirun -n N ./amun.x -i ./params.in`, `mpirun -n N ./amun.x -i ./params.in`,
where N is the number of processors to use. where N is the number of processors to use.
Reading data
============
By default, the code uses new XML+binary snapshot data format. It can also be
forced by setting parameter **snapshot_format** to **xml**.
In order to read produced data in this format, you will need to install the
provided Python module. Simply change to **python/** directory and run
`python setup.py install --user`
to install the module in your home directory.
Import the module in your python script using
`from amunpy import *`,
and then initiate the interface using
`snapshot = AmunXML(<path to the snapshot directory>)`
and read desired variable using
`var = snapshot.dataset(<variable>)`.
The function **dataset()** returns rescaled uniform mesh variable as NumPy
array.
If you want to read data from HDF5 snapshot, just use
`var = amun_dataset(<snapshot HDF5 file>, <variable>)`.

7
bitbucket-pipelines.yml
View File

@ -1,9 +1,4 @@
# This is a sample build configuration for Other. image: debian
# Check our guides at https://confluence.atlassian.com/x/5Q4SMw for more examples.
# Only use spaces to indent your .yml configuration.
# -----
# You can specify a custom docker image from Docker Hub as your build environment.
image: atlassian/default-image:2
pipelines: pipelines:
default: default:

7
build/make.default
View File

@ -20,11 +20,10 @@ SIGNALS = N
# #
MPI = N MPI = N
# output data format # OUTPUT - the format of the alternative data output;
# supported alternative formats: HDF5
# #
# OUTPUT - the format of output data; so far only HDF5 OUTPUT =
#
OUTPUT = HDF5
#------------------------------------------------------------------------------- #-------------------------------------------------------------------------------
# #

4
build/makefile
View File

@ -68,7 +68,9 @@ FFLAGS += ${CPPPREFIX}-DNDIMS=${NDIMS}
# output data format # output data format
# #
FFLAGS += ${CPPPREFIX}-D${OUTPUT} ifneq (,$(findstring HDF5,$(OUTPUT)))
FFLAGS += ${CPPPREFIX}-DHDF5
endif
# add module path to compiler options # add module path to compiler options
# #

403
python/amunpy.py
View File

@ -43,6 +43,409 @@ try:
scipy_available = True scipy_available = True
except ImportError: except ImportError:
scipy_available = False scipy_available = False
try:
from xxhash import *
xxhash_available = True
except ImportError:
xxhash_available = False
class AmunXML:
"""AMUN XML snapshot class"""
def __init__(self, path, shrink = 1, interpolation = 'rebin'):
import xml.etree.ElementTree as ET
if op.isdir(path):
mfile = op.join(path, 'metadata.xml')
if op.exists(mfile):
self.path = path
self.metafile = mfile
self.format = "xml"
else:
print("%s does not contain 'metadata.xml'!" % path)
return None
else:
mfile = op.basename(path)
if mfile == 'metadata.xml':
self.path = op.dirname(path)
self.metafile = path
self.format = "xml"
else:
print("%s is not 'metadata.xml'!" % path)
return None
tree = ET.parse(self.metafile)
root = tree.getroot()
if root.tag == 'AMUNFile':
self.attributes = dict()
for item in root.iter('Attribute'):
if item.attrib['type'] == 'double':
self.attributes[item.attrib['name']] = float(item.text)
elif item.attrib['type'] == 'integer':
self.attributes[item.attrib['name']] = int(item.text)
else:
self.attributes[item.attrib['name']] = item.text
else:
print("%s is not an AMUN snapshot!" % self.metafile)
return None
if 'variables' in self.attributes:
variables = []
for var in self.attributes['variables'].split():
variables.append(var)
variables.append('level')
if all(v in variables for v in ['velx','vely','velz']):
variables.append('velo')
variables.append('divv')
variables.append('vort')
if all(v in variables for v in ['magx','magy','magz']):
variables.append('magn')
variables.append('divb')
variables.append('curr')
if 'pres' in variables and 'gamma' in self.attributes:
variables.append('eint')
if all(v in variables for v in ['dens','pres']):
variables.append('temp')
if self.attributes['eqsys'] in [ 'hd', 'mhd' ] \
and all(v in variables for v in ['dens','velo']):
variables.append('ekin')
if self.attributes['eqsys'] in [ 'mhd', 'srmhd' ] \
and 'emag' in variables:
variables.append('emag')
if all(v in variables for v in ['ekin','eint']):
variables.append('etot')
if self.attributes['eqsys'] in [ 'srhd', 'srmhd' ] \
and 'velo' in variables:
variables.append('lore')
else:
print("Cannot determine available variables from %s!" % self.metafile)
return None
self.attributes['xlen'] = self.attributes['xmax'] - self.attributes['xmin']
self.attributes['ylen'] = self.attributes['ymax'] - self.attributes['ymin']
if self.attributes['ndims'] == 3:
self.attributes['zlen'] = self.attributes['zmax'] - self.attributes['zmin']
self.cell_size = dict()
for l in range(self.attributes['maxlev']):
n = self.attributes['xblocks'] * self.attributes['ncells'] * 2**l
h = self.attributes['xlen'] / n
self.cell_size[l+1] = h
self.shrink = shrink
self.interpolation = interpolation
self.variables = variables
def readdblocks(self, no, var):
import numpy as np
import xml.etree.ElementTree as ET
dfile = op.join(self.path, 'datablock_%s_%06d.bin' % (var, no))
if op.exists(dfile):
data = np.fromfile(dfile, dtype='float64')
else:
print("Cannot find file %s!" % dfile)
return None
if xxhash_available:
xfile = op.join(self.path, 'datablocks_%06d.xml' % (no))
tree = ET.parse(xfile)
root = tree.getroot()
for item in root.findall('./Hashes/Attribute'):
if item.attrib['name'] == var:
digest = item.text.split(':')[1].lower()
if digest != xxh64(data).hexdigest():
print("File '%s' seems to be corrupted!" % (dfile))
return data
def dataset(self, var):
import numpy as np
import xml.etree.ElementTree as ET
if not var in self.variables:
print("Dataset '%s' not available in the snapshot!" % (var))
print("Available variables are: ", self.variables)
return None
nd = self.attributes['ndims']
nn = self.attributes['bcells']
nc = self.attributes['ncells']
ng = self.attributes['nghosts']
ml = self.attributes['maxlev']
nx = self.attributes['xblocks']
ny = self.attributes['yblocks']
nm = self.attributes['nleafs']
if nd == 3:
nz = self.attributes['zblocks']
rm = np.array([nx, ny, nz])
bm = np.array([nc, nc, nc])
else:
rm = np.array([nx, ny])
bm = np.array([nc, nc])
dm = rm[0:nd] * int(nc * 2**(ml - 1) / self.shrink)
arr = np.zeros(dm[:])
dfile = op.join(self.path, 'metablock_fields.bin')
if op.exists(dfile):
mset = np.fromfile(dfile, dtype='int32')
if xxhash_available:
if 'fields' in self.attributes:
digest = self.attributes['fields'].split(':')[1].lower()
if digest != xxh64(mset).hexdigest():
print("File '%s' seems to be corrupted!" % (dfile))
else:
return None
nf = int(mset.size / nm)
mset = np.reshape(mset, [nf,nm])
level = dict()
coord = dict()
for n in range(nm):
level[mset[0,n]] = mset[ 1,n]
coord[mset[0,n]] = mset[2:5,n]
for n in range(self.attributes['nprocs']):
dfile = op.join(self.path, 'datablocks_%06d.xml' % (n))
tree = ET.parse(dfile)
root = tree.getroot()
nb = - 1
digest = None
for item in root.findall('./DataBlocks/Attribute'):
if item.attrib['name'] == 'dblocks':
nb = int(item.text)
for item in root.findall('./Hashes/Attribute'):
if item.attrib['name'] == 'ids':
digest = item.text.split(':')[1].lower()
if nb > 0:
cm = [ nn ]*nd
cm.append( nb )
dfile = op.join(self.path, 'datablock_ids_%06d.bin' % (n))
if op.exists(dfile):
ids = np.fromfile(dfile, dtype='int32')
if xxhash_available:
if digest != xxh64(ids).hexdigest():
print("File '%s' seems to be corrupted!" % (dfile))
else:
return None
if var == 'level':
dset = np.zeros(cm)
for p in range(nb):
dset[...,p] = level[ids[p]]
elif var == 'velo':
tmp = self.readdblocks(n, 'velx')
dset = tmp**2
tmp = self.readdblocks(n, 'vely')
dset += tmp**2
tmp = self.readdblocks(n, 'velz')
dset += tmp**2
dset = np.reshape(np.sqrt(dset), cm)
elif var == 'magn':
tmp = self.readdblocks(n, 'magx')
dset = tmp**2
tmp = self.readdblocks(n, 'magy')
dset += tmp**2
tmp = self.readdblocks(n, 'magz')
dset += tmp**2
dset = np.reshape(np.sqrt(dset), cm)
elif var == 'ekin':
tmp = self.readdblocks(n, 'velx')
dset = tmp**2
tmp = self.readdblocks(n, 'vely')
dset += tmp**2
tmp = self.readdblocks(n, 'velz')
dset += tmp**2
tmp = self.readdblocks(n, 'dens')
dset *= tmp
dset *= 0.5
dset = np.reshape(dset, cm)
elif var == 'emag':
tmp = self.readdblocks(n, 'magx')
dset = tmp**2
tmp = self.readdblocks(n, 'magy')
dset += tmp**2
tmp = self.readdblocks(n, 'magz')
dset += tmp**2
dset *= 0.5
dset = np.reshape(dset, cm)
elif var == 'eint':
dset = self.readdblocks(n, 'pres')
dset *= 1.0 / (self.attributes('gamma') - 1.0)
dset = np.reshape(dset, cm)
elif var == 'temp':
dset = self.readdblocks(n, 'pres')
tmp = self.readdblocks(n, 'pres')
dset /= tmp
dset = np.reshape(dset, cm)
elif var == 'etot':
tmp = self.readdblocks(n, 'velx')
dset = tmp**2
tmp = self.readdblocks(n, 'vely')
dset += tmp**2
tmp = self.readdblocks(n, 'velz')
dset += tmp**2
tmp = self.readdblocks(n, 'dens')
dset *= tmp
tmp = self.readdblocks(n, 'pres')
dset += 2.0 / (self.attributes('gamma') - 1.0) * tmp
if 'magn' in self.variables:
tmp = self.readdblocks(n, 'magx')
dset = tmp**2
tmp = self.readdblocks(n, 'magy')
dset += tmp**2
tmp = self.readdblocks(n, 'magz')
dset += tmp**2
dset *= 0.5
dset = np.reshape(dset, cm)
elif var == 'lorentz':
tmp = self.readdblocks(n, 'velx')
dset = tmp**2
tmp = self.readdblocks(n, 'vely')
dset += tmp**2
tmp = self.readdblocks(n, 'velz')
dset += tmp**2
dset = 1.0 / np.sqrt(1.0 - dset)
dset = np.reshape(dset, cm)
elif var == 'divv':
m = nd - 1
tmp = np.reshape(self.readdblocks(n, 'velx'), cm)
dset = (np.roll(tmp, -1, axis = m) \
- np.roll(tmp, 1, axis = m))
m -= 1
tmp = np.reshape(self.readdblocks(n, 'vely'), cm)
dset += (np.roll(tmp, -1, axis = m) \
- np.roll(tmp, 1, axis = m))
m -= 1
if m >= 0:
tmp = np.reshape(self.readdblocks(n, 'velz'), cm)
dset += (np.roll(tmp, -1, axis = m) \
- np.roll(tmp, 0, axis = m))
dset *= 0.5
for p in range(nb):
dset[...,p] /= self.cell_size[level[ids[p]]]
elif var == 'divb':
m = nd - 1
tmp = np.reshape(self.readdblocks(n, 'magx'), cm)
dset = (np.roll(tmp, -1, axis = m) \
- np.roll(tmp, 1, axis = m))
m -= 1
tmp = np.reshape(self.readdblocks(n, 'magy'), cm)
dset += (np.roll(tmp, -1, axis = m) \
- np.roll(tmp, 1, axis = m))
m -= 1
if m >= 0:
tmp = np.reshape(self.readdblocks(n, 'magz'), cm)
dset += (np.roll(tmp, -1, axis = m) \
- np.roll(tmp, 0, axis = m))
dset *= 0.5
for p in range(nb):
dset[...,p] /= self.cell_size[level[ids[p]]]
elif var == 'vort':
if nd == 3:
tmp = np.reshape(self.readdblocks(n, 'velx'), cm)
wy = (np.roll(tmp, -1, axis = 0) \
- np.roll(tmp, 1, axis = 0))
wz = (np.roll(tmp, 1, axis = 1) \
- np.roll(tmp, -1, axis = 1))
tmp = np.reshape(self.readdblocks(n, 'vely'), cm)
wz += (np.roll(tmp, -1, axis = 2) \
- np.roll(tmp, 1, axis = 2))
wx = (np.roll(tmp, 1, axis = 0) \
- np.roll(tmp, -1, axis = 0))
tmp = np.reshape(self.readdblocks(n, 'velz'), cm)
wx += (np.roll(tmp, -1, axis = 1) \
- np.roll(tmp, 1, axis = 1))
wy += (np.roll(tmp, 1, axis = 2) \
- np.roll(tmp, -1, axis = 2))
else:
tmp = np.reshape(self.readdblocks(n, 'velx'), cm)
wz = (np.roll(tmp, 1, axis = 0) \
- np.roll(tmp, -1, axis = 0))
tmp = np.reshape(self.readdblocks(n, 'vely'), cm)
wz += (np.roll(tmp, -1, axis = 1) \
- np.roll(tmp, 1, axis = 1))
tmp = np.reshape(self.readdblocks(n, 'velz'), cm)
wx = (np.roll(tmp, -1, axis = 0) \
- np.roll(tmp, 1, axis = 0))
wy = (np.roll(tmp, 1, axis = 1) \
- np.roll(tmp, -1, axis = 1))
dset = 0.5 * np.sqrt(wx**2 + wy**2 + wz**2)
for p in range(nb):
dset[...,p] /= self.cell_size[level[ids[p]]]
elif var == 'curr':
if nd == 3:
tmp = np.reshape(self.readdblocks(n, 'magx'), cm)
wy = (np.roll(tmp, -1, axis = 0) \
- np.roll(tmp, 1, axis = 0))
wz = (np.roll(tmp, 1, axis = 1) \
- np.roll(tmp, -1, axis = 1))
tmp = np.reshape(self.readdblocks(n, 'magy'), cm)
wz += (np.roll(tmp, -1, axis = 2) \
- np.roll(tmp, 1, axis = 2))
wx = (np.roll(tmp, 1, axis = 0) \
- np.roll(tmp, -1, axis = 0))
tmp = np.reshape(self.readdblocks(n, 'magz'), cm)
wx += (np.roll(tmp, -1, axis = 1) \
- np.roll(tmp, 1, axis = 1))
wy += (np.roll(tmp, 1, axis = 2) \
- np.roll(tmp, -1, axis = 2))
else:
tmp = np.reshape(self.readdblocks(n, 'magx'), cm)
wz = (np.roll(tmp, 1, axis = 0) \
- np.roll(tmp, -1, axis = 0))
tmp = np.reshape(self.readdblocks(n, 'magy'), cm)
wz += (np.roll(tmp, -1, axis = 1) \
- np.roll(tmp, 1, axis = 1))
tmp = np.reshape(self.readdblocks(n, 'magz'), cm)
wx = (np.roll(tmp, -1, axis = 0) \
- np.roll(tmp, 1, axis = 0))
wy = (np.roll(tmp, 1, axis = 1) \
- np.roll(tmp, -1, axis = 1))
dset = 0.5 * np.sqrt(wx**2 + wy**2 + wz**2)
for p in range(nb):
dset[...,p] /= self.cell_size[level[ids[p]]]
else:
dset = self.readdblocks(n, var)
dset = np.reshape(dset, cm)
for p in range(nb):
ii = ids[p]
nl = 2**(ml - level[ii])
if nl <= self.shrink:
method = 'rebin'
else:
method = self.interpolation
cm = bm[0:nd] * int(nl / self.shrink)
il = coord[ii][0:nd] * cm[0:nd]
iu = il + cm[0:nd]
if nd == 3:
ib, jb, kb = il[0], il[1], il[2]
ie, je, ke = iu[0], iu[1], iu[2]
arr[kb:ke,jb:je,ib:ie] = interpolate(dset[:,:,:,p], cm, ng, method = method)
else:
ib, jb = il[0], il[1]
ie, je = iu[0], iu[1]
arr[jb:je,ib:ie] = interpolate(dset[:,:,p], cm, ng, method = method)
return arr
def amun_compatible(fname): def amun_compatible(fname):

7
python/setup.py
View File

@ -3,12 +3,15 @@ from setuptools import setup
setup( setup(
name='amunpy', name='amunpy',
description='Python Interface fo AMUN snapshots', description='Python Interface fo AMUN snapshots',
version='0.1', version='0.2',
author='Grzegorz Kowal', author='Grzegorz Kowal',
author_email='grzegorz@amuncode.org', author_email='grzegorz@amuncode.org',
url='https://www.amuncode.org/', url='https://www.amuncode.org/',
license='GPLv3', license='GPLv3',
py_modules=['amunpy'], py_modules=['amunpy'],
install_requires=['h5py', 'numpy'], install_requires=['h5py', 'numpy'],
extra_require=['scipy'] extras_require={
"digest": ['xxhash'],
"interpolation": ['scipy'],
}
) )

7
sources/blocks.F90
View File

@ -1181,6 +1181,13 @@ module blocks
! !
if (status == 0) then if (status == 0) then
! reset arrays
!
pdata%u0 = 0.0d+00
pdata%u1 = 0.0d+00
pdata%q = 0.0d+00
pdata%f = 0.0d+00
! initiate the conserved variable pointer ! initiate the conserved variable pointer
! !
pdata%u => pdata%u0 pdata%u => pdata%u0

39
sources/driver.F90
View File

@ -50,6 +50,8 @@ program amun
use evolution , only : print_evolution use evolution , only : print_evolution
use evolution , only : advance, new_time_step use evolution , only : advance, new_time_step
use evolution , only : step, time, dt use evolution , only : step, time, dt
use forcing , only : initialize_forcing, finalize_forcing
use forcing , only : print_forcing
use gravity , only : initialize_gravity, finalize_gravity use gravity , only : initialize_gravity, finalize_gravity
use helpers , only : print_welcome, print_section, print_parameter use helpers , only : print_welcome, print_section, print_parameter
use integrals , only : initialize_integrals, finalize_integrals use integrals , only : initialize_integrals, finalize_integrals
@ -278,7 +280,7 @@ program amun
! the restart snapshot, otherwise, read them from the parameter file ! the restart snapshot, otherwise, read them from the parameter file
! !
if (restart_from_snapshot()) then if (restart_from_snapshot()) then
call read_snapshot_parameter("problem", problem, status) call read_snapshot_parameter("problem", problem , status)
if (check_status(status /= 0)) go to 3800 if (check_status(status /= 0)) go to 3800
call read_snapshot_parameter("eqsys" , eqsys , status) call read_snapshot_parameter("eqsys" , eqsys , status)
if (check_status(status /= 0)) go to 3800 if (check_status(status /= 0)) go to 3800
@ -288,7 +290,11 @@ program amun
if (check_status(status /= 0)) go to 3800 if (check_status(status /= 0)) go to 3800
call read_snapshot_parameter("maxlev" , toplev , status) call read_snapshot_parameter("maxlev" , toplev , status)
if (check_status(status /= 0)) go to 3800 if (check_status(status /= 0)) go to 3800
call read_snapshot_parameter("bdims" , bdims , status) call read_snapshot_parameter("xblocks", bdims(1), status)
if (check_status(status /= 0)) go to 3800
call read_snapshot_parameter("yblocks", bdims(2), status)
if (check_status(status /= 0)) go to 3800
call read_snapshot_parameter("zblocks", bdims(3), status)
if (check_status(status /= 0)) go to 3800 if (check_status(status /= 0)) go to 3800
call read_snapshot_parameter("xmin" , xmin , status) call read_snapshot_parameter("xmin" , xmin , status)
if (check_status(status /= 0)) go to 3800 if (check_status(status /= 0)) go to 3800
@ -356,8 +362,14 @@ program amun
! initialize the remaining modules ! initialize the remaining modules
! !
call initialize_random(1, 0) call initialize_random(1, 0, nproc, status)
if (check_status(status /= 0)) go to 3700 if (check_status(status /= 0)) then
if (master) then
write(error_unit,"('[AMUN::program]: ', a)") &
"Problem initializing RANDOM module!"
end if
go to 3700
end if
call initialize_equations(eqsys, eos, master, status) call initialize_equations(eqsys, eos, master, status)
if (check_status(status /= 0)) then if (check_status(status /= 0)) then
if (master) then if (master) then
@ -487,19 +499,28 @@ program amun
end if end if
go to 2100 go to 2100
end if end if
call initialize_forcing(master, status)
if (check_status(status /= 0)) then
if (master) then
write(error_unit,"('[AMUN::program]: ', a)") &
"Problem initializing FORCING module!"
end if
go to 2000
end if
call initialize_integrals(master, nrun, status) call initialize_integrals(master, nrun, status)
if (check_status(status /= 0)) then if (check_status(status /= 0)) then
if (master) then if (master) then
write(error_unit,"('[AMUN::program]: ', a)") & write(error_unit,"('[AMUN::program]: ', a)") &
"Problem initializing INTEGRALS module!" "Problem initializing INTEGRALS module!"
end if end if
go to 2000 go to 1900
end if end if
! print module information ! print module information
! !
call print_equations(master) call print_equations(master)
call print_sources(master) call print_sources(master)
call print_forcing(master)
call print_coordinates(master) call print_coordinates(master)
call print_boundaries(master) call print_boundaries(master)
call print_shapes(master) call print_shapes(master)
@ -756,12 +777,18 @@ program amun
! finalize modules ! finalize modules
! !
2000 continue 1900 continue
call finalize_integrals(status) call finalize_integrals(status)
if (check_status(status /= 0) .and. master) then if (check_status(status /= 0) .and. master) then
write(error_unit,"('[AMUN::program]: ', a)") & write(error_unit,"('[AMUN::program]: ', a)") &
"Problem finalizing INTEGRALS module!" "Problem finalizing INTEGRALS module!"
end if end if
2000 continue
call finalize_forcing(status)
if (check_status(status /= 0) .and. master) then
write(error_unit,"('[AMUN::program]: ', a)") &
"Problem finalizing FORCING module!"
end if
2100 continue 2100 continue
call finalize_evolution(status) call finalize_evolution(status)
if (check_status(status /= 0) .and. master) then if (check_status(status /= 0) .and. master) then

2
sources/equations.F90
View File

@ -179,7 +179,7 @@ module equations
! isothermal speed of sound and its second power ! isothermal speed of sound and its second power
! !
real(kind=8) , save :: csnd = 1.0d+00, csnd2 = 1.0d+00 real(kind=8), save :: csnd = 1.0d+00, csnd2 = 1.0d+00
! maximum speed in the system ! maximum speed in the system
! !

7
sources/evolution.F90
View File

@ -414,6 +414,7 @@ module evolution
! !
use blocks , only : set_blocks_update use blocks , only : set_blocks_update
use coordinates, only : toplev use coordinates, only : toplev
use forcing , only : update_forcing, forcing_enabled
use mesh , only : update_mesh use mesh , only : update_mesh
! local variables are not implicit by default ! local variables are not implicit by default
@ -441,7 +442,11 @@ module evolution
! !
call evolve() call evolve()
! chec if we need to perform the refinement step ! add forcing contribution
!
if (forcing_enabled) call update_forcing(time, dt)
! check if we need to perform the refinement step
! !
if (toplev > 1) then if (toplev > 1) then

1817
sources/forcing.F90 Normal file
View File

File diff suppressed because it is too large Load Diff

561
sources/hash.F90 Normal file
View File

@ -0,0 +1,561 @@
!!******************************************************************************
!!
!! This file is part of the AMUN source code, a program to perform
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!! adaptive mesh.
!!
!! Copyright (C) 2012-2020 Yann Collet
!! Copyright (C) 2020 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! module: HASH
!!
!! This module provides 64-bit version of the xxHash64 by Yann Collet.
!! This is a Fortran implementation based on the XXH64 specification
!! published at
!! https://github.com/Cyan4973/xxHash/blob/dev/doc/xxhash_spec.md
!!
!! For additional info, see
!! http://www.xxhash.com or https://github.com/Cyan4973/xxHash
!!
!!******************************************************************************
!
module hash
! module variables are not implicit by default
!
implicit none
! hash parameters
!
integer(kind=8), parameter :: seed = 0_8
integer(kind=8), parameter :: prime1 = -7046029288634856825_8, &
prime2 = -4417276706812531889_8, &
prime3 = 1609587929392839161_8, &
prime4 = -8796714831421723037_8, &
prime5 = 2870177450012600261_8
! by default everything is private
!
private
! declare public subroutines
!
public :: xxh64_integer, xxh64_long, xxh64_double, xxh64_complex
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
contains
!
!===============================================================================
!!
!!*** PRIVATE SUBROUTINES ****************************************************
!!
!===============================================================================
!
!
!===============================================================================
!
! function XXH64_INTEGER:
! ----------------------
!
! Function calculates XXH64 hash for a given integer vector.
!
! Arguments:
!
! n - the size of the input vector;
! data - the input vactor of integer values;
!
!===============================================================================
!
integer(kind=8) function xxh64_integer(n, data) result(hash)
implicit none
! subroutine arguments
!
integer(kind=4) , intent(in) :: n
integer(kind=4), dimension(n), intent(in) :: data
! local variables
!
integer(kind=8) :: remain, offset
! local arrays
!
integer(kind=8), dimension(4) :: lane, chk
!-------------------------------------------------------------------------------
!
hash = 0_8
offset = 1
remain = n
if (remain >= 8) then
lane(1) = seed + prime1 + prime2
lane(2) = seed + prime2
lane(3) = seed + 0_8
lane(4) = seed - prime1
do while (remain >= 8)
chk(1:4) = transfer(data(offset:offset+7), chk)
lane(1) = xxh64_round(lane(1), chk(1))
lane(2) = xxh64_round(lane(2), chk(2))
lane(3) = xxh64_round(lane(3), chk(3))
lane(4) = xxh64_round(lane(4), chk(4))
offset = offset + 8
remain = remain - 8
end do
hash = xxh64_rotl(lane(1), 1) + xxh64_rotl(lane(2), 7) + &
xxh64_rotl(lane(3), 12) + xxh64_rotl(lane(4), 18)
hash = xxh64_merge(hash, lane(1))
hash = xxh64_merge(hash, lane(2))
hash = xxh64_merge(hash, lane(3))
hash = xxh64_merge(hash, lane(4))
else
hash = seed + prime5
end if
hash = hash + int(4 * n, kind = 8)
do while (remain >= 2)
chk(1) = transfer(data(offset:offset+1), chk(1))
hash = ieor(hash, xxh64_round(0_8, chk(1)))
hash = xxh64_rotl(hash, 27)
hash = hash * prime1 + prime4
offset = offset + 2
remain = remain - 2
end do
if (remain == 1) then
hash = ieor(hash, data(offset) * prime1)
hash = xxh64_rotl(hash, 23)
hash = hash * prime2 + prime3
offset = offset + 1
remain = remain - 1
end if
hash = xxh64_aval(hash)
return
!-------------------------------------------------------------------------------
!
end function xxh64_integer
!
!===============================================================================
!
! function XXH64_LONG:
! -------------------
!
! Function calculates XXH64 hash for a given long integer vector.
!
! Arguments:
!
! n - the size of the input vector;
! data - the input vactor of real values;
!
!===============================================================================
!
integer(kind=8) function xxh64_long(n, data) result(hash)
implicit none
! subroutine arguments
!
integer(kind=4) , intent(in) :: n
integer(kind=8), dimension(n), intent(in) :: data
! local variables
!
integer(kind=8) :: remain, offset
! local arrays
!
integer(kind=8), dimension(4) :: lane, chk
!-------------------------------------------------------------------------------
!
hash = 0_8
offset = 1
remain = n
if (remain >= 4) then
lane(1) = seed + prime1 + prime2
lane(2) = seed + prime2
lane(3) = seed + 0_8
lane(4) = seed - prime1
do while (remain >= 4)
chk(1:4) = transfer(data(offset:offset+3), chk)
lane(1) = xxh64_round(lane(1), chk(1))
lane(2) = xxh64_round(lane(2), chk(2))
lane(3) = xxh64_round(lane(3), chk(3))
lane(4) = xxh64_round(lane(4), chk(4))
offset = offset + 4
remain = remain - 4
end do
hash = xxh64_rotl(lane(1), 1) + xxh64_rotl(lane(2), 7) + &
xxh64_rotl(lane(3), 12) + xxh64_rotl(lane(4), 18)
hash = xxh64_merge(hash, lane(1))
hash = xxh64_merge(hash, lane(2))
hash = xxh64_merge(hash, lane(3))
hash = xxh64_merge(hash, lane(4))
else
hash = seed + prime5
end if
hash = hash + int(8 * n, kind = 8)
do while (remain >= 1)
hash = ieor(hash, xxh64_round(0_8, transfer(data(offset), 0_8)))
hash = xxh64_rotl(hash, 27)
hash = hash * prime1 + prime4
offset = offset + 1
remain = remain - 1
end do
hash = xxh64_aval(hash)
return
!-------------------------------------------------------------------------------
!
end function xxh64_long
!
!===============================================================================
!
! function XXH64_DOUBLE:
! ---------------------
!
! Function calculates XXH64 hash for a given double precision vector.
!
! Arguments:
!
! n - the size of the input vector;
! data - the input vactor of real values;
!
!===============================================================================
!
integer(kind=8) function xxh64_double(n, data) result(hash)
implicit none
! subroutine arguments
!
integer(kind=4) , intent(in) :: n
real(kind=8), dimension(n), intent(in) :: data
! local variables
!
integer(kind=8) :: remain, offset
! local arrays
!
integer(kind=8), dimension(4) :: lane, chk
!-------------------------------------------------------------------------------
!
hash = 0_8
offset = 1
remain = n
if (remain >= 4) then
lane(1) = seed + prime1 + prime2
lane(2) = seed + prime2
lane(3) = seed + 0_8
lane(4) = seed - prime1
do while (remain >= 4)
chk(1:4) = transfer(data(offset:offset+3), chk)
lane(1) = xxh64_round(lane(1), chk(1))
lane(2) = xxh64_round(lane(2), chk(2))
lane(3) = xxh64_round(lane(3), chk(3))
lane(4) = xxh64_round(lane(4), chk(4))
offset = offset + 4
remain = remain - 4
end do
hash = xxh64_rotl(lane(1), 1) + xxh64_rotl(lane(2), 7) + &
xxh64_rotl(lane(3), 12) + xxh64_rotl(lane(4), 18)
hash = xxh64_merge(hash, lane(1))
hash = xxh64_merge(hash, lane(2))
hash = xxh64_merge(hash, lane(3))
hash = xxh64_merge(hash, lane(4))
else
hash = seed + prime5
end if
hash = hash + int(8 * n, kind = 8)
do while (remain >= 1)
hash = ieor(hash, xxh64_round(0_8, transfer(data(offset), 0_8)))
hash = xxh64_rotl(hash, 27)
hash = hash * prime1 + prime4
offset = offset + 1
remain = remain - 1
end do
hash = xxh64_aval(hash)
return
!-------------------------------------------------------------------------------
!
end function xxh64_double
!
!===============================================================================
!
! function XXH64_COMPLEX:
! ----------------------
!
! Function calculates XXH64 hash for a given double precision complex vector.
!
! Arguments:
!
! n - the size of the input vector;
! data - the input vactor of real values;
!
!===============================================================================
!
integer(kind=8) function xxh64_complex(n, data) result(hash)
implicit none
! subroutine arguments
!
integer(kind=4) , intent(in) :: n
complex(kind=8), dimension(n), intent(in) :: data
! local variables
!
integer(kind=8) :: remain, offset
! local arrays
!
integer(kind=8), dimension(4) :: lane, chk
!-------------------------------------------------------------------------------
!
hash = 0_8
offset = 1
remain = n
if (remain >= 2) then
lane(1) = seed + prime1 + prime2
lane(2) = seed + prime2
lane(3) = seed + 0_8
lane(4) = seed - prime1
do while (remain >= 2)
chk(1:4) = transfer(data(offset:offset+1), chk)
lane(1) = xxh64_round(lane(1), chk(1))
lane(2) = xxh64_round(lane(2), chk(2))
lane(3) = xxh64_round(lane(3), chk(3))
lane(4) = xxh64_round(lane(4), chk(4))
offset = offset + 2
remain = remain - 2
end do
hash = xxh64_rotl(lane(1), 1) + xxh64_rotl(lane(2), 7) + &
xxh64_rotl(lane(3), 12) + xxh64_rotl(lane(4), 18)
hash = xxh64_merge(hash, lane(1))
hash = xxh64_merge(hash, lane(2))
hash = xxh64_merge(hash, lane(3))
hash = xxh64_merge(hash, lane(4))
else
hash = seed + prime5
end if
hash = hash + int(16 * n, kind = 8)
if (remain == 1) then
hash = ieor(hash, xxh64_round(0_8, transfer(dreal(data(offset)), 0_8)))
hash = xxh64_rotl(hash, 27)
hash = hash * prime1 + prime4
hash = ieor(hash, xxh64_round(0_8, transfer(dimag(data(offset)), 0_8)))
hash = xxh64_rotl(hash, 27)
hash = hash * prime1 + prime4
offset = offset + 1
remain = remain - 1
end if
hash = xxh64_aval(hash)
return
!-------------------------------------------------------------------------------
!
end function xxh64_complex
!
!===============================================================================
!
! function XXH64_ROUND:
! --------------------
!
! Function processes one stripe of the input data updating
! the correponding lane.
!
! Arguments:
!
! lane - the lane;
! input - the 8-byte data to process;
!
!===============================================================================
!
integer(kind=8) function xxh64_round(lane, input)
implicit none
! subroutine arguments
!
integer(kind=8), intent(in) :: lane, input
!-------------------------------------------------------------------------------
!
xxh64_round = lane + (input * prime2)
xxh64_round = xxh64_rotl(xxh64_round, 31)
xxh64_round = xxh64_round * prime1
return
!-------------------------------------------------------------------------------
!
end function xxh64_round
!
!===============================================================================
!
! function XXH64_MERGE:
! --------------------
!
! Function performs merging of the given lane in to the hash.
!
! Arguments:
!
! hash - the hash to merge to;
! lane - the lane being merged;
!
!===============================================================================
!
integer(kind=8) function xxh64_merge(hash, lane)
implicit none
! subroutine arguments
!
integer(kind=8), intent(in) :: hash, lane
!-------------------------------------------------------------------------------
!
xxh64_merge = ieor(hash, xxh64_round(0_8, lane))
xxh64_merge = xxh64_merge * prime1 + prime4
return
!-------------------------------------------------------------------------------
!
end function xxh64_merge
!
!===============================================================================
!
! function XXH64_AVAL:
! -------------------
!
! Function calculates the final mix of the hash.
!
! Arguments:
!
! hash - the hash to mix;
!
!===============================================================================
!
integer(kind=8) function xxh64_aval(hash)
implicit none
! subroutine arguments
!
integer(kind=8), intent(in) :: hash
!-------------------------------------------------------------------------------
!
xxh64_aval = hash
xxh64_aval = ieor(xxh64_aval, ishft(xxh64_aval, -33)) * prime2
xxh64_aval = ieor(xxh64_aval, ishft(xxh64_aval, -29)) * prime3
xxh64_aval = ieor(xxh64_aval, ishft(xxh64_aval, -32))
return
!-------------------------------------------------------------------------------
!
end function xxh64_aval
!
!===============================================================================
!
! function XXH64_ROTL:
! -------------------
!
! Function calculates the rotation of the input 8-byte word by a given amount.
!
! Arguments:
!
! byte - the byte to be rotates;
! amount - the amount by which rotate the input byte;
!
!===============================================================================
!
integer(kind=8) function xxh64_rotl(byte, amount)
implicit none
! subroutine arguments
!
integer(kind=8), intent(in) :: byte
integer(kind=4), intent(in) :: amount
!-------------------------------------------------------------------------------
!
xxh64_rotl = ior(ishft(byte, amount), ishft(byte, amount - 64))
return
!-------------------------------------------------------------------------------
!
end function xxh64_rotl
!===============================================================================
!
end module hash

4
sources/helpers.F90
View File

@ -215,9 +215,9 @@ module helpers
msg = trim(adjustl(description)) msg = trim(adjustl(description))
if (value >= 0.0d+00) then if (value >= 0.0d+00) then
write(*,"(4x,a26,1x,'=',es9.2)") msg, value write(*,"(4x,a26,1x,'=',es12.5)") msg, value
else else
write(*,"(4x,a26,1x,'=',1x,es9.2)") msg, value write(*,"(4x,a26,1x,'=',1x,es12.5)") msg, value
end if end if
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------

13
sources/integrals.F90
View File

@ -180,9 +180,10 @@ module integrals
! write the integral file header ! write the integral file header
! !
write(funit,"('#',a8,10(1x,a18))") 'step', 'time', 'dt' & write(funit,"('#',a8,13(1x,a18))") 'step', 'time', 'dt' &
, 'mass', 'momx', 'momy', 'momz' & , 'mass', 'momx', 'momy', 'momz' &
, 'ener', 'ekin', 'emag', 'eint' , 'ener', 'ekin', 'emag', 'eint' &
, 'einj', 'erat', 'arms'
write(funit,"('#')") write(funit,"('#')")
! depending on the append parameter, choose the right file initialization for ! depending on the append parameter, choose the right file initialization for
@ -369,6 +370,7 @@ module integrals
use equations , only : ien, imx, imy, imz use equations , only : ien, imx, imy, imz
use equations , only : magnetized, gamma, csnd use equations , only : magnetized, gamma, csnd
use evolution , only : step, time, dtn use evolution , only : step, time, dtn
use forcing , only : einj, rinj, arms
use mpitools , only : master use mpitools , only : master
#ifdef MPI #ifdef MPI
use mpitools , only : reduce_sum_real_array use mpitools , only : reduce_sum_real_array
@ -518,6 +520,11 @@ module integrals
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
end if end if
! sum up the injected energy and injection rate
!
inarr( 9) = einj
inarr(10) = rinj
! get average, minimum and maximum values of density ! get average, minimum and maximum values of density
! !
#if NDIMS == 3 #if NDIMS == 3
@ -881,7 +888,7 @@ module integrals
! write down the integrals and statistics to appropriate files ! write down the integrals and statistics to appropriate files
! !
if (master) then if (master) then
write(funit,"(i9,10(1x,1es18.8e3))") step, time, dtn, inarr(1:8) write(funit,"(i9,13(1x,1es18.8e3))") step, time, dtn, inarr(1:10), arms
write(sunit,"(i9,23(1x,1es18.8e3))") step, time & write(sunit,"(i9,23(1x,1es18.8e3))") step, time &
, avarr(1), mnarr(1), mxarr(1) & , avarr(1), mnarr(1), mxarr(1) &
, avarr(2), mnarr(2), mxarr(2) & , avarr(2), mnarr(2), mxarr(2) &

2729
sources/io.F90
View File

File diff suppressed because it is too large Load Diff

13
sources/mpitools.F90
View File

@ -1345,7 +1345,7 @@ module mpitools
! include external procedures and variables ! include external procedures and variables
! !
use iso_fortran_env, only : error_unit use iso_fortran_env, only : error_unit
use mpi , only : mpi_real8, mpi_sum, mpi_success use mpi , only : mpi_double_complex, mpi_sum, mpi_success
! local variables are not implicit by default ! local variables are not implicit by default
! !
@ -1354,12 +1354,12 @@ module mpitools
! subroutine arguments ! subroutine arguments
! !
integer , intent(in) :: n integer , intent(in) :: n
complex, dimension(n), intent(inout) :: cbuf complex(kind=8), dimension(n), intent(inout) :: cbuf
integer , intent(out) :: iret integer , intent(out) :: iret
! local variables ! local variables
! !
real(kind=8), dimension(n) :: rbuf, ibuf, tbuf complex(kind=8), dimension(n) :: tbuf
! local parameters ! local parameters
! !
@ -1377,10 +1377,7 @@ module mpitools
call start_timer(imm) call start_timer(imm)
#endif /* PROFILE */ #endif /* PROFILE */
tbuf(:) = real(cbuf(:)) call mpi_allreduce(cbuf, tbuf, n, mpi_double_complex, mpi_sum, comm, iret)
call mpi_allreduce(tbuf, rbuf, n, mpi_real8, mpi_sum, comm, iret)
tbuf(:) = aimag(cbuf(:))
call mpi_allreduce(tbuf, ibuf, n, mpi_real8, mpi_sum, comm, iret)
#ifdef PROFILE #ifdef PROFILE
! stop time accounting for the MPI reduce ! stop time accounting for the MPI reduce
@ -1390,7 +1387,7 @@ module mpitools
! substitute the result ! substitute the result
! !
cbuf(1:n) = cmplx(rbuf(1:n), ibuf(1:n)) cbuf(:) = tbuf(:)
! check if the operation was successful ! check if the operation was successful
! !

56
sources/parameters.F90
View File

@ -75,7 +75,7 @@ module parameters
! declare public subroutines ! declare public subroutines
! !
public :: read_parameters, finalize_parameters public :: read_parameters, finalize_parameters
public :: get_parameter public :: get_parameter_file, get_parameter
#ifdef MPI #ifdef MPI
public :: redistribute_parameters public :: redistribute_parameters
#endif /* MPI */ #endif /* MPI */
@ -231,6 +231,60 @@ module parameters
! !
!=============================================================================== !===============================================================================
! !
! subroutine GET_PARAMETER_FILE:
! -----------------------------
!
! Subroutine returns the full path to the parameter file.
!
! Arguments:
!
! pfile - the parameter full file path;
! status - the status value, 0 for success;
!
!===============================================================================
!
subroutine get_parameter_file(pfile, status)
! import external procedures
!
use iso_fortran_env, only : error_unit
! local variables are not implicit by default
!
implicit none
! input and output variables
!
character(len=*), intent(out) :: pfile
integer , intent(out) :: status
! local variables
!
character(len=mlen) :: tfile
! local parameters
!
character(len=*), parameter :: loc = 'PARAMETERS::get_parameter_file()'
!
!-------------------------------------------------------------------------------
!
status = 0
if (len(pfile) <= mlen) then
write(pfile,"(a)") trim(adjustl(fname))
else
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Parameter file path too long for subroutine argument!"
write(tfile,"(a)") trim(adjustl(fname))
write(pfile,"(a)") tfile(1:len(pfile))
status = 1
end if
!-------------------------------------------------------------------------------
!
end subroutine get_parameter_file
!
!===============================================================================
!
! subroutine GET_PARAMETERS_NUMBER: ! subroutine GET_PARAMETERS_NUMBER:
! -------------------------------- ! --------------------------------
! !

12
sources/problems.F90
View File

@ -1515,7 +1515,7 @@ module problems
use equations , only : nv use equations , only : nv
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl
use parameters , only : get_parameter use parameters , only : get_parameter
use random , only : randomn use random , only : randsym
! local variables are not implicit by default ! local variables are not implicit by default
! !
@ -1658,10 +1658,10 @@ module problems
! add a random seed velocity component ! add a random seed velocity component
! !
do i = 1, nn do i = 1, nn
q(ivx,i) = q(ivx,i) + vper * randomn() q(ivx,i) = q(ivx,i) + vper * randsym()
q(ivy,i) = q(ivy,i) + vper * randomn() q(ivy,i) = q(ivy,i) + vper * randsym()
#if NDIMS == 3 #if NDIMS == 3
q(ivz,i) = q(ivz,i) + vper * randomn() q(ivz,i) = q(ivz,i) + vper * randsym()
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
end do end do
@ -1729,7 +1729,7 @@ module problems
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl
use equations , only : csnd2 use equations , only : csnd2
use parameters , only : get_parameter use parameters , only : get_parameter
use random , only : randomn use random , only : randsym
! local variables are not implicit by default ! local variables are not implicit by default
! !
@ -1882,7 +1882,7 @@ module problems
! !
if (vper /= 0.0d+00) then if (vper /= 0.0d+00) then
do i = 1, nn do i = 1, nn
q(ivy,i) = q(ivy,i) + vper * randomn() q(ivy,i) = q(ivy,i) + vper * randsym()
end do end do
end if end if

467
sources/random.F90
View File

@ -27,8 +27,10 @@
!! !!
!! References: !! References:
!! !!
!! [1] Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating !! [1] David Blackman and Sebastiano Vigna,
!! random variables', J. Statist. Software, v5(8) !! "Scrambled linear pseudorandom number generators", 2019
!! http://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf
!! [2] http://prng.di.unimi.it/
!! !!
!!****************************************************************************** !!******************************************************************************
! !
@ -47,7 +49,7 @@ module random
#ifdef PROFILE #ifdef PROFILE
! timer indices ! timer indices
! !
integer , save :: iri, irc integer, save :: iri, irg
#endif /* PROFILE */ #endif /* PROFILE */
! random generator type ! random generator type
@ -56,9 +58,15 @@ module random
! variables to store seeds for random generator ! variables to store seeds for random generator
! !
integer , save :: kp = 0 integer , save :: kp = 0, np = 0
integer , save :: nseeds, lseed integer , save :: nseeds, lseed
integer(kind=4), dimension(:), allocatable, save :: seeds integer , save :: state
integer(kind=8), dimension(:,:), allocatable, save :: seeds
integer(kind=4), dimension(:), allocatable, save :: ss
! integer to real conversion parameters
!
real(kind=8), parameter :: i64tor8 = 2.0d+00**(-53)
! by default everything is private ! by default everything is private
! !
@ -67,7 +75,8 @@ module random
! declare public subroutines ! declare public subroutines
! !
public :: initialize_random, finalize_random public :: initialize_random, finalize_random
public :: nseeds, set_seeds, get_seeds, randomu, randomz, randomn public :: nseeds, set_seeds, get_seeds, gentype
public :: randuni, randsym, randnorz
contains contains
! !
@ -76,77 +85,89 @@ module random
! subroutine INITIALIZE_RANDOM: ! subroutine INITIALIZE_RANDOM:
! ---------------------------- ! ----------------------------
! !
! subroutine initializes random number generator; ! Subroutine initializes random number generator.
!
! Arguments:
!
! nthreads - the number of threads;
! nthread - the thread number;
! nproc - the process number;
! !
!=============================================================================== !===============================================================================
! !
subroutine initialize_random(nthreads, nthread) subroutine initialize_random(nthreads, nthread, nproc, status)
! obtain required variables from other modules
!
use parameters, only : get_parameter use parameters, only : get_parameter
! declare all variables as implicit
!
implicit none implicit none
! subroutine arguments ! subroutine arguments
! !
integer, intent(in) :: nthreads, nthread integer, intent(in) :: nthreads, nthread, nproc
integer, intent(out) :: status
! local variables ! local variables
! !
integer :: i integer :: i
real(kind=4) :: r
! local arrays
!
integer(kind=8), dimension(4) :: s
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE status = 0
! set timer descriptions
! #ifdef PROFILE
call set_timer('random:: initialization' , iri) call set_timer('random:: initialization' , iri)
call set_timer('random:: number generation', irc) call set_timer('random:: number generation', irg)
! start accounting time for the random number generator
!
call start_timer(iri) call start_timer(iri)
#endif /* PROFILE */ #endif /* PROFILE */
! set the thread number ! set local copies of the thread and process numbers,
! determine the number of seeds
! !
kp = nthread kp = nthread
np = nproc
! obtain the generator type
!
call get_parameter("gentype", gentype)
! calculate the number of seeds
!
nseeds = nthreads nseeds = nthreads
lseed = nseeds - 1 lseed = nseeds - 1
! allocate seeds for random number generator ! allocate seeds
! !
allocate(seeds(0:lseed)) allocate(seeds(4, 0:lseed), stat = status)
! prepare the seeds depending on the type of generator if (status == 0) then
! prepare seeds depending on the type of generator
! !
call get_parameter("gentype", gentype)
select case(gentype) select case(gentype)
case('random') case('random')
do i = 0, lseed state = 1234567890
call random_number(r) s(1) = splitmix64()
seeds(i) = 123456789 * r s(2) = splitmix64()
s(3) = splitmix64()
s(4) = splitmix64()
do i = 1, nseeds * np
call jump(s(:))
end do
seeds(:,0) = s(:)
do i = 1, lseed
call jump(s(:))
seeds(:,i) = s(:)
end do end do
case default case default
r = 0.1234567890123456789 state = 1234567890
do i = 0, lseed seeds(1,:) = splitmix64()
seeds(i) = 123456789 * r seeds(2,:) = splitmix64()
end do seeds(3,:) = splitmix64()
seeds(4,:) = splitmix64()
end select end select
end if
#ifdef PROFILE #ifdef PROFILE
! stop accounting time for the random number generator
!
call stop_timer(iri) call stop_timer(iri)
#endif /* PROFILE */ #endif /* PROFILE */
@ -159,31 +180,23 @@ module random
! subroutine FINALIZE_RANDOM: ! subroutine FINALIZE_RANDOM:
! -------------------------- ! --------------------------
! !
! subroutine releases memory allocated by random number generator variables; ! Subroutine releases memory allocated by random number generator variables.
! !
!=============================================================================== !===============================================================================
! !
subroutine finalize_random() subroutine finalize_random()
! declare all variables as implicit
!
implicit none implicit none
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE #ifdef PROFILE
! start accounting time for the random number generator
!
call start_timer(iri) call start_timer(iri)
#endif /* PROFILE */ #endif /* PROFILE */
! deallocate seeds if they are allocated
!
if (allocated(seeds)) deallocate(seeds) if (allocated(seeds)) deallocate(seeds)
#ifdef PROFILE #ifdef PROFILE
! stop accounting time for the random number generator
!
call stop_timer(iri) call stop_timer(iri)
#endif /* PROFILE */ #endif /* PROFILE */
@ -196,59 +209,70 @@ module random
! subroutine SET_SEEDS: ! subroutine SET_SEEDS:
! -------------------- ! --------------------
! !
! subroutine sets the seeds from the given array; ! Subroutine sets the seeds.
!
! Arguments:
!
! n - the number of seeds;
! sds - the new seeds;
! gen - the flag indication if the seeds have to be regenerated;
! !
!=============================================================================== !===============================================================================
! !
subroutine set_seeds(np, nseeds, generate) subroutine set_seeds(n, sds, gen)
! declare all variables as implicit
!
implicit none implicit none
! input arguments ! input arguments
! !
integer , intent(in) :: np integer , intent(in) :: n
integer(kind=4), dimension(np), intent(in) :: nseeds integer(kind=8), dimension(4,n), intent(in) :: sds
logical , intent(in) :: generate logical , intent(in) :: gen
! local variables ! local variables
! !
integer :: i, l integer :: i
real(kind=4) :: r
! local arrays
!
integer(kind=8), dimension(4) :: s
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE #ifdef PROFILE
! start accounting time for the random number generator
!
call start_timer(iri) call start_timer(iri)
#endif /* PROFILE */ #endif /* PROFILE */
! set the seeds only if the input array and seeds have the same sizes ! set the seeds only if the input array and seeds have the same sizes
! !
l = min(lseed, np - 1) if (gen .or. n /= nseeds) then
seeds(0:l) = nseeds(1:l+1)
select case(gentype) select case(gentype)
case('random') case('random')
if (generate) then state = 1234567890
do i = 0, lseed s(1) = splitmix64()
call random_number(r) s(2) = splitmix64()
seeds(i) = 123456789 * r s(3) = splitmix64()
s(4) = splitmix64()
do i = 1, nseeds * np
call jump(s(:))
end do end do
else seeds(:,0) = s(:)
do i = l + 1, lseed do i = 1, lseed
call random_number(r) call jump(s(:))
seeds(i) = 123456789 * r seeds(:,i) = s(:)
end do end do
end if
case default case default
seeds(l+1:lseed) = seeds(0) state = 1234567890
seeds(1,:) = splitmix64()
seeds(2,:) = splitmix64()
seeds(3,:) = splitmix64()
seeds(4,:) = splitmix64()
end select end select
else
seeds(:,:) = sds(:,:)
end if
#ifdef PROFILE #ifdef PROFILE
! stop accounting time for the random number generator
!
call stop_timer(iri) call stop_timer(iri)
#endif /* PROFILE */ #endif /* PROFILE */
@ -261,33 +285,33 @@ module random
! subroutine GET_SEEDS: ! subroutine GET_SEEDS:
! -------------------- ! --------------------
! !
! subroutine returns the seeds through an array; ! Subroutine returns the seeds.
!
! Arguments:
!
! sds - the return seeds array;
! !
!=============================================================================== !===============================================================================
! !
subroutine get_seeds(rseeds) subroutine get_seeds(sds)
! declare all variables as implicit
!
implicit none implicit none
! output arguments ! output arguments
! !
integer(kind=4), dimension(nseeds), intent(out) :: rseeds integer(kind=8), dimension(4,nseeds), intent(out) :: sds
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE #ifdef PROFILE
! start accounting time for the random number generator
!
call start_timer(iri) call start_timer(iri)
#endif /* PROFILE */ #endif /* PROFILE */
rseeds(1:nseeds) = seeds(0:lseed) if (size(sds,1) == 4 .and. size(sds,2) == nseeds) then
sds(:,:) = seeds(:,:)
end if
#ifdef PROFILE #ifdef PROFILE
! stop accounting time for the random number generator
!
call stop_timer(iri) call stop_timer(iri)
#endif /* PROFILE */ #endif /* PROFILE */
@ -297,165 +321,268 @@ module random
! !
!=============================================================================== !===============================================================================
! !
! function RANDOMU: ! function SPLITMIX64:
! ---------------- ! -------------------
! !
! function generates uniformly distributed random numbers in range 0.0..1.0; ! Pseudorandom Number Generator for floating point numbers based on
! SplitMix64 algorithm [1], used for initialization of xoshiro256p()
! generator.
!
! References:
!
! [1] David Blackman and Sebastiano Vigna,
! "Scrambled linear pseudorandom number generators", 2019
! http://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf
! !
!=============================================================================== !===============================================================================
! !
function randomu() result(val) integer(kind=8) function splitmix64()
! declare all variables as implicit
!
implicit none implicit none
! output variables ! local parameters
! !
real(kind=8) :: val integer(kind=8), parameter :: a = -7046029254386353131_8
integer(kind=8), parameter :: b = -4658895280553007687_8
integer(kind=8), parameter :: c = -7723592293110705685_8
!
!-------------------------------------------------------------------------------
!
state = state + a
splitmix64 = state
splitmix64 = ieor(splitmix64, ishft(splitmix64, -30)) * b
splitmix64 = ieor(splitmix64, ishft(splitmix64, -27)) * c
splitmix64 = ieor(splitmix64, ishft(splitmix64, -31))
return
!-------------------------------------------------------------------------------
!
end function splitmix64
!
!===============================================================================
!
! function XOSHIRO256P:
! --------------------
!
! Pseudorandom Number Generator for floating point numbers based on
! xoshiro256+ algorithm [1]. It return a 256-bit random integer.
!
! References:
!
! [1] David Blackman and Sebastiano Vigna,
! "Scrambled linear pseudorandom number generators", 2019
! http://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf
!
!===============================================================================
!
integer(kind=8) function xoshiro256p()
implicit none
! local variables ! local variables
! !
integer(kind=4) :: jz, jsr integer(kind=8) :: temp
!
!-------------------------------------------------------------------------------
!
xoshiro256p = ishft(seeds(1,kp) + seeds(4,kp), -11)
temp = ishft(seeds(2,kp), 17)
seeds(3,kp) = ieor(seeds(3,kp), seeds(1,kp))
seeds(4,kp) = ieor(seeds(4,kp), seeds(2,kp))
seeds(2,kp) = ieor(seeds(2,kp), seeds(3,kp))
seeds(1,kp) = ieor(seeds(1,kp), seeds(4,kp))
seeds(3,kp) = ieor(seeds(3,kp), temp)
seeds(4,kp) = ior(ishft(seeds(4,kp), 45), ishft(seeds(4,kp), -19))
return
!-------------------------------------------------------------------------------
!
end function xoshiro256p
!
!===============================================================================
!
! subroutine JUMP:
! ---------------
!
! The jump function for the xoshiro256p() generator used to generate
! non-overlapping seed subsequences for parallel computations.
!
! Arguments:
!
! sds - temporary seeds array;
!
! References:
!
! [1] David Blackman and Sebastiano Vigna,
! "Scrambled linear pseudorandom number generators", 2019
! http://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf
!
!===============================================================================
!
subroutine jump(sds)
implicit none
! arguments
!
integer(kind=8), dimension(4), intent(inout) :: sds
! local variables
!
integer(kind=8), dimension(4) :: s
integer(kind=8) :: temp
integer :: i, b
! local parameters
!
integer(kind=8), dimension(4), parameter :: jmp = (/ 1733541517147835066_8 &
, -3051731464161248980_8 &
, -6244198995065845334_8 &
, 4155657270789760540_8 /)
!
!-------------------------------------------------------------------------------
!
s(:) = 0_8
do i = 1, 4
do b = 1, 64
if (iand(jmp(i), ishft(1_8, b)) /= 0_8) then
s(1) = ieor(s(1), sds(1))
s(2) = ieor(s(2), sds(2))
s(3) = ieor(s(3), sds(3))
s(4) = ieor(s(4), sds(4))
end if
temp = ishft(sds(2), 17)
sds(3) = ieor(sds(3), sds(1))
sds(4) = ieor(sds(4), sds(2))
sds(2) = ieor(sds(2), sds(3))
sds(1) = ieor(sds(1), sds(4))
sds(3) = ieor(sds(3), temp)
sds(4) = ior(ishft(sds(4), 45), ishft(sds(4), -19))
end do
end do
sds(:) = s(:)
return
!-------------------------------------------------------------------------------
!
end subroutine jump
!
!===============================================================================
!
! function RANDUNI:
! ----------------
!
! Function generates uniformly distributed floating point real random number
! in range 0.0..1.0 using xoshiro256p() generator.
!
!===============================================================================
!
real(kind=8) function randuni()
implicit none
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE #ifdef PROFILE
! start accounting time for the random number generation call start_timer(irg)
!
call start_timer(irc)
#endif /* PROFILE */ #endif /* PROFILE */
jsr = seeds(kp) randuni = i64tor8 * xoshiro256p()
jz = jsr
jsr = ieor( jsr, ishft( jsr, 13 ) )
jsr = ieor( jsr, ishft( jsr, -17 ) )
jsr = ieor( jsr, ishft( jsr, 5 ) )
seeds(kp) = jsr
val = 0.5d+00 + 0.23283064365d-09 * (jz + jsr)
#ifdef PROFILE #ifdef PROFILE
! stop accounting time for the random number generation call stop_timer(irg)
!
call stop_timer(irc)
#endif /* PROFILE */ #endif /* PROFILE */
return return
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end function randomu end function randuni
! !
!=============================================================================== !===============================================================================
! !
! function RANDOMZ: ! function RANDSYM:
! ---------------- ! ----------------
! !
! function generates uniformly distributed random numbers in range -0.5..0.5; ! Function generates uniformly distributed floating point random number
! in range -1.0..1.0 using xoshiro256p() generator.
! !
!=============================================================================== !===============================================================================
! !
function randomz() result(val) real(kind=8) function randsym()
! declare all variables as implicit
!
implicit none implicit none
! output variables
!
real(kind=8) :: val
! local variables
!
integer(kind=4) :: jz, jsr
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE #ifdef PROFILE
! start accounting time for the random number generation call start_timer(irg)
!
call start_timer(irc)
#endif /* PROFILE */ #endif /* PROFILE */
jsr = seeds(kp) randsym = 2.0d+00 * (i64tor8 * xoshiro256p() - 5.0d-01)
jz = jsr
jsr = ieor( jsr, ishft( jsr, 13 ) )
jsr = ieor( jsr, ishft( jsr, -17 ) )
jsr = ieor( jsr, ishft( jsr, 5 ) )
seeds(kp) = jsr
val = 0.23283064365d-09 * (jz + jsr)
#ifdef PROFILE #ifdef PROFILE
! stop accounting time for the random number generation call stop_timer(irg)
!
call stop_timer(irc)
#endif /* PROFILE */ #endif /* PROFILE */
return return
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end function randomz end function randsym
! !
!=============================================================================== !===============================================================================
! !
! function RANDOMN: ! function RANDNORZ:
! ---------------- ! -----------------
! !
! function generates uniformly distributed random numbers in range -1.0..1.0; ! Function generates complex floating point number with a normal
! distribution (μ = 0, σ = 1) using Marsaglia polar method [1].
!
! References:
!
! [1] Marsaglia, G. & Bray, T. A.,
! "A Convenient Method for Generating Normal Variables",
! SIAM Review, vol. 6, no. 3, pp. 260-264,
! https://doi.org/10.1137/1006063
! !
!=============================================================================== !===============================================================================
! !
function randomn() result(val) complex(kind=8) function randnorz()
! declare all variables as implicit
!
implicit none implicit none
! output variables logical :: c
! real(kind=8) :: u, v, s
real(kind=8) :: val
! local variables
!
integer(kind=4) :: jz, jsr
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
#ifdef PROFILE #ifdef PROFILE
! start accounting time for the random number generation call start_timer(irg)
!
call start_timer(irc)
#endif /* PROFILE */ #endif /* PROFILE */
jsr = seeds(kp) c = .true.
jz = jsr do while(c)
u = randsym()
v = randsym()
s = u**2 + v**2
c = s <= 0.0d+00 .or. s >= 1.0d+00
end do
jsr = ieor( jsr, ishft( jsr, 13 ) ) randnorz = dsqrt(- 2.0d+00 * dlog(s) / s) * dcmplx(u, v)
jsr = ieor( jsr, ishft( jsr, -17 ) )
jsr = ieor( jsr, ishft( jsr, 5 ) )
seeds(kp) = jsr
val = 0.46566128730d-09 * (jz + jsr)
#ifdef PROFILE #ifdef PROFILE
! stop accounting time for the random number generation call stop_timer(irg)
!
call stop_timer(irc)
#endif /* PROFILE */ #endif /* PROFILE */
return return
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end function randomn end function randnorz
!=============================================================================== !===============================================================================
! !