diff --git a/README.md b/README.md index b406f4a..0b49a21 100644 --- a/README.md +++ b/README.md @@ -18,12 +18,15 @@ following features are already implemented: reconstructions, * Riemann solvers of Roe- and HLL-types (HLL, HLLC, and HLLD), * standard boundary conditions: periodic, open, reflective, hydrostatic, etc. +* turbulence driving using Alvelius or Ornstein–Uhlenbeck methods, * viscous and resistive source terms, -* suppor for passive scalars (up to 100), -* data stored in the HDF5 format, +* support for passive scalars (up to 100), +* data stored in internal XML+binary or HDF5 format, +* Python interface to read snapshots in both formats, * MPI parallelization, * 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 the terms of the GNU General Public License as published by the Free Software @@ -53,19 +56,21 @@ Requirements version 18.10 or newer, - [Intel Fortran](https://software.intel.com/en-us/fortran-compilers) compiler version 9.0 or newer. -* [HDF5 libraries](https://www.hdfgroup.org/solutions/hdf5/), tested with - version 1.8 or newer. -* [OpenMPI](https://www.open-mpi.org/) for parallel runs, tested with version - 1.8 or newer. +* Optional, but recommended, [OpenMPI](https://www.open-mpi.org/) for parallel + runs, tested with version 1.8 or newer. +* Optional [HDF5 libraries](https://www.hdfgroup.org/solutions/hdf5/), tested + 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 ===================== -If the HDF5 libraries are not installed in the default location, i.e. in the -system directory **/usr**, make sure that the environment variable _HDF5DIR_ is -set in your **~/.bashrc** (or **~/.cshrc**) and pointing to the location where -the HDF5 libraries have been installed. +If you need to use the HDF5 libraries and they are not installed in the default +location, i.e. in the system directory **/usr**, make sure that the environment +variable _HDF5DIR_ is set in your **~/.bashrc** (or **~/.cshrc**) and pointing +to the location where the HDF5 libraries have been installed. Compilation @@ -103,4 +108,29 @@ In order to run serial version, just type in your terminal: In order to run parallel version (after compiling the code with MPI support), type in your terminal: `mpirun -n N ./amun.x -i ./params.in`, -where N is the number of processors to use. \ No newline at end of file +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()` +and read desired variable using + `var = snapshot.dataset()`. + +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(, )`. diff --git a/bitbucket-pipelines.yml b/bitbucket-pipelines.yml index 003e634..a52627f 100644 --- a/bitbucket-pipelines.yml +++ b/bitbucket-pipelines.yml @@ -1,9 +1,4 @@ -# This is a sample build configuration for Other. -# 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 +image: debian pipelines: default: @@ -22,4 +17,4 @@ pipelines: - make clean - make MPI=Y NDIMS=2 - make clean - - make MPI=Y NDIMS=3 \ No newline at end of file + - make MPI=Y NDIMS=3 diff --git a/build/make.default b/build/make.default index 144ac77..46af636 100644 --- a/build/make.default +++ b/build/make.default @@ -20,11 +20,10 @@ SIGNALS = 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 = HDF5 +OUTPUT = #------------------------------------------------------------------------------- # @@ -39,4 +38,4 @@ NDIMS = 2 # # path for temporary files during compilation (uncomment if you wish to use it) # -#OBJSDIR := /tmp/${USER}/amun-code/objects \ No newline at end of file +#OBJSDIR := /tmp/${USER}/amun-code/objects diff --git a/build/makefile b/build/makefile index df7c65d..d5fa593 100644 --- a/build/makefile +++ b/build/makefile @@ -68,7 +68,9 @@ FFLAGS += ${CPPPREFIX}-DNDIMS=${NDIMS} # output data format # -FFLAGS += ${CPPPREFIX}-D${OUTPUT} +ifneq (,$(findstring HDF5,$(OUTPUT))) +FFLAGS += ${CPPPREFIX}-DHDF5 +endif # add module path to compiler options # diff --git a/python/amunpy.py b/python/amunpy.py index 853d304..7c21985 100644 --- a/python/amunpy.py +++ b/python/amunpy.py @@ -39,10 +39,413 @@ import numpy as np import os.path as op import sys try: - from scipy.interpolate import splrep, splev, interp1d, PchipInterpolator, pchip_interpolate - scipy_available = True + from scipy.interpolate import splrep, splev, interp1d, PchipInterpolator, pchip_interpolate + scipy_available = True 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): diff --git a/python/setup.py b/python/setup.py index 5a87b83..a0392a0 100644 --- a/python/setup.py +++ b/python/setup.py @@ -3,12 +3,15 @@ from setuptools import setup setup( name='amunpy', description='Python Interface fo AMUN snapshots', - version='0.1', + version='0.2', author='Grzegorz Kowal', author_email='grzegorz@amuncode.org', url='https://www.amuncode.org/', license='GPLv3', py_modules=['amunpy'], install_requires=['h5py', 'numpy'], - extra_require=['scipy'] + extras_require={ + "digest": ['xxhash'], + "interpolation": ['scipy'], + } ) diff --git a/sources/blocks.F90 b/sources/blocks.F90 index d24aecc..42e2091 100644 --- a/sources/blocks.F90 +++ b/sources/blocks.F90 @@ -1181,6 +1181,13 @@ module blocks ! 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 ! pdata%u => pdata%u0 diff --git a/sources/driver.F90 b/sources/driver.F90 index 0c3f7f5..bcdafaa 100644 --- a/sources/driver.F90 +++ b/sources/driver.F90 @@ -50,6 +50,8 @@ program amun use evolution , only : print_evolution use evolution , only : advance, new_time_step 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 helpers , only : print_welcome, print_section, print_parameter use integrals , only : initialize_integrals, finalize_integrals @@ -278,30 +280,34 @@ program amun ! the restart snapshot, otherwise, read them from the parameter file ! 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 - call read_snapshot_parameter("eqsys" , eqsys , status) + call read_snapshot_parameter("eqsys" , eqsys , status) if (check_status(status /= 0)) go to 3800 - call read_snapshot_parameter("eos" , eos , status) + call read_snapshot_parameter("eos" , eos , status) if (check_status(status /= 0)) go to 3800 - call read_snapshot_parameter("ncells" , ncells , status) + call read_snapshot_parameter("ncells" , ncells , status) 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 - 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("xmin" , xmin , status) + call read_snapshot_parameter("yblocks", bdims(2), status) if (check_status(status /= 0)) go to 3800 - call read_snapshot_parameter("xmax" , xmax , status) + call read_snapshot_parameter("zblocks", bdims(3), status) if (check_status(status /= 0)) go to 3800 - call read_snapshot_parameter("ymin" , ymin , status) + call read_snapshot_parameter("xmin" , xmin , status) if (check_status(status /= 0)) go to 3800 - call read_snapshot_parameter("ymax" , ymax , status) + call read_snapshot_parameter("xmax" , xmax , status) + if (check_status(status /= 0)) go to 3800 + call read_snapshot_parameter("ymin" , ymin , status) + if (check_status(status /= 0)) go to 3800 + call read_snapshot_parameter("ymax" , ymax , status) if (check_status(status /= 0)) go to 3800 #if NDIMS == 3 - call read_snapshot_parameter("zmin" , zmin , status) + call read_snapshot_parameter("zmin" , zmin , status) if (check_status(status /= 0)) go to 3800 - call read_snapshot_parameter("zmax" , zmax , status) + call read_snapshot_parameter("zmax" , zmax , status) if (check_status(status /= 0)) go to 3800 #endif /* NDIMS == 3 */ else @@ -356,8 +362,14 @@ program amun ! initialize the remaining modules ! - call initialize_random(1, 0) - if (check_status(status /= 0)) go to 3700 + call initialize_random(1, 0, nproc, status) + 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) if (check_status(status /= 0)) then if (master) then @@ -487,19 +499,28 @@ program amun end if go to 2100 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) if (check_status(status /= 0)) then if (master) then write(error_unit,"('[AMUN::program]: ', a)") & "Problem initializing INTEGRALS module!" end if - go to 2000 + go to 1900 end if ! print module information ! call print_equations(master) call print_sources(master) + call print_forcing(master) call print_coordinates(master) call print_boundaries(master) call print_shapes(master) @@ -756,12 +777,18 @@ program amun ! finalize modules ! - 2000 continue + 1900 continue call finalize_integrals(status) if (check_status(status /= 0) .and. master) then write(error_unit,"('[AMUN::program]: ', a)") & "Problem finalizing INTEGRALS module!" 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 call finalize_evolution(status) if (check_status(status /= 0) .and. master) then diff --git a/sources/equations.F90 b/sources/equations.F90 index a2e02ff..4551481 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -179,7 +179,7 @@ module equations ! 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 ! diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 1c62072..12e8b73 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -414,6 +414,7 @@ module evolution ! use blocks , only : set_blocks_update use coordinates, only : toplev + use forcing , only : update_forcing, forcing_enabled use mesh , only : update_mesh ! local variables are not implicit by default @@ -441,7 +442,11 @@ module evolution ! 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 diff --git a/sources/forcing.F90 b/sources/forcing.F90 new file mode 100644 index 0000000..fe80363 --- /dev/null +++ b/sources/forcing.F90 @@ -0,0 +1,1817 @@ +!!****************************************************************************** +!! +!! 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) 2017-2020 Grzegorz Kowal +!! +!! 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 . +!! +!!***************************************************************************** +!! +!! module: FORCING +!! +!! This module provides energy injection, e.g. turbulence driving, supernova +!! explosions, etc. +!! +!!***************************************************************************** +! +module forcing + +#ifdef PROFILE +! import external subroutines +! + use timers, only : set_timer, start_timer, stop_timer +#endif /* PROFILE */ + +! module variables are not implicit by default +! + implicit none + +#ifdef PROFILE +! timer indices +! + integer, save :: ifi, ifu +#endif /* PROFILE */ + +! interfaces for procedure pointers +! + abstract interface + subroutine update_forcing_iface(t, dt) + real(kind=8), intent(in) :: t, dt + end subroutine + end interface + +! pointers to the forcing methods +! + procedure(update_forcing_iface), pointer, save :: update_forcing => null() + +! flag indicating if the energy injection is enabled and if it is done in +! the Fourier space +! + logical, save :: forcing_enabled = .false. + logical, save :: forcing_fourier = .false. + logical, save :: energy_profile = .false. + +! implemented injection methods +! + integer, parameter :: injection_none = 0 + integer, parameter :: injection_eddy = 1 + integer, parameter :: injection_oh = 2 + integer, parameter :: injection_alvelius = 3 + +! implemented driving force spectral profiles +! + integer, parameter :: profile_normal = 1 + integer, parameter :: profile_parabolic = 2 + integer, parameter :: profile_powerlaw = 3 + +! driving force method and spectral profile +! + integer, save :: injection_method = injection_none + integer, save :: injection_profile = profile_normal + +! forcing parameters +! + real(kind=8), save :: power = 1.0d+00 + real(kind=8), save :: rate = 1.0d+03 + real(kind=8), save :: amp = 1.0d-03 + real(kind=8), save :: del = 1.0d-01 + +! Fourier profile parameters +! + integer , save :: nmodes = 0 + real(kind=8), save :: fpower = 1.0d+00 + real(kind=8), save :: lscale = 5.0d-01 + real(kind=8), save :: tscale = 1.0d+00 + real(kind=8), save :: vscale = 5.0d-01 + real(kind=8), save :: ascale = 5.0d-01 + real(kind=8), save :: kf = 2.0d+00 + real(kind=8), save :: kl = 0.0d+00 + real(kind=8), save :: ku = 4.0d+00 + real(kind=8), save :: kc = 5.0d-01 + real(kind=8), save :: pidx = - 5.0d+00 / 3.0d+00 + real(kind=8), save :: fmin = 1.0d-08 + real(kind=8), save :: norm = 3.0d+00 / sqrt(2.0d+00) + +! module parameters +! + real(kind=8), save :: dinj = 0.0d+00 + real(kind=8), save :: einj = 0.0d+00 + real(kind=8), save :: rinj = 0.0d+00 + real(kind=8), save :: arms = 0.0d+00 + +! bound for the driving mode wave number +! + integer, save :: kmax = 4 + +! arrays for driving mode positions and amplitudes +! + integer , dimension(:,: ) , allocatable :: kmodes + real(kind=8), dimension(:) , allocatable :: fmodes + +! vector of driving directions in Alfvelius method +! + real(kind=8), dimension(:,:) , allocatable :: e1vecs, e2vecs + +! array for driving mode coefficients +! + complex(kind=8), dimension(:,:), allocatable :: fcoefs + +! velocity Fourier coefficients in Alfvelius method +! + complex(kind=8), dimension(:,:), allocatable :: vcoefs + +! by default everything is private +! + private + +! declare public subroutines +! + public :: initialize_forcing, finalize_forcing, print_forcing + public :: update_forcing + public :: forcing_enabled, einj, rinj, arms, nmodes, fcoefs + +!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! + contains +! +!=============================================================================== +!! +!!*** PUBLIC SUBROUTINES ***************************************************** +!! +!=============================================================================== +! +!=============================================================================== +! +! subroutine INITIALIZE_FORCING: +! ----------------------------- +! +! Subroutine initializes module FORCING. +! +! Arguments: +! +! verbose - flag determining if the subroutine should be verbose; +! status - return flag of the procedure execution status; +! +!=============================================================================== +! + subroutine initialize_forcing(verbose, status) + +! import external procedures and variables +! + use constants , only : pi2 + use iso_fortran_env, only : error_unit + use parameters , only : get_parameter + use random , only : randuni, randnorz + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + logical, intent(in) :: verbose + integer, intent(out) :: status + +! local variables +! + character(len=32) :: wfmt = "(1x,'WARNING:',1x,a)" + character(len=64) :: injection = "none" + character(len=64) :: profile_type = "gauss" + character(len=64) :: profile_energy = "off" + logical :: test + integer :: i, j, k = 0, l, k2 + real(kind=8) :: kl2, ku2, kv2, kv + real(kind=8) :: fa, fi, uu, phi + +! local vectors +! + real(kind=8), dimension(NDIMS) :: u, v + +! allocatable arrays +! + integer, dimension(:), allocatable :: kcount + +! local parameters +! + character(len=*), parameter :: loc = 'FORCING::initialize_forcing()' +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! set timer descriptions +! + call set_timer('forcing:: initialization', ifi) + call set_timer('forcing:: update' , ifu) + +! start accounting time for initialization/finalization +! + call start_timer(ifi) +#endif /* PROFILE */ + +! reset the status flag +! + status = 0 + +! obtain the chosen injection method +! + call get_parameter("driving_method", injection) + +! select the energy injection method +! + select case(trim(injection)) + case("eddies", "eddy") + forcing_enabled = .true. + injection_method = injection_eddy + update_forcing => update_forcing_eddies + + call get_parameter("injection_power", power) + call get_parameter("injection_rate" , rate ) + call get_parameter("eddy_amplitude" , amp ) + call get_parameter("eddy_width" , del ) + + case("Ornstein–Uhlenbeck", "ornstein–uhlenbeck", "oh") + forcing_enabled = .true. + forcing_fourier = .true. + injection_method = injection_oh + update_forcing => update_forcing_oh + +! get the driving power +! + call get_parameter("driving_timescale", tscale) + if (tscale <= 0.0d+00) then + if (verbose) then + write(*,*) + write(*,wfmt) "'driving_timescale' must be larger than zero!" + write(*,wfmt) "Resetting it to the default value: 1.0!" + end if + tscale = 1.0d+00 + end if + + case("Alvelius", "alvelius", "al") + forcing_enabled = .true. + forcing_fourier = .true. + injection_method = injection_alvelius + update_forcing => update_forcing_alvelius + +! get the driving power +! + call get_parameter("driving_power", fpower) + if (fpower <= 0.0d+00) then + if (verbose) then + write(*,*) + write(*,wfmt) "'driving_power' must be larger than zero!" + write(*,wfmt) "Resetting it to the default value: 1.0!" + end if + fpower = 1.0d+00 + end if + + case default + injection_method = injection_none + if (verbose .and. injection /= "none") then + write(*,*) + write(*,wfmt) "Unknown injection method! Forcing disabled!" + end if + end select + +! initialization of the Fourier space driving +! + if (forcing_fourier) then + + call get_parameter("driving_wavenumber" , kf) + call get_parameter("driving_profile_width", kc) + call get_parameter("driving_profile_lower", kl) + call get_parameter("driving_profile_upper", ku) + call get_parameter("amplitude_threshold" , fmin) + + if (kf < 1.0d+00) then + if (verbose) then + write(*,*) + write(*,wfmt) "'driving_wavenumber' must be larger than 1.0!" + write(*,wfmt) "Resetting it to the default value: 2.0!" + end if + kf = 2.0d+00 + end if + if (kc <= 0.0d+00) then + if (verbose) then + write(*,*) + write(*,wfmt) "'driving_profile_width' must be larger than 0.0!" + write(*,wfmt) "Resetting it to the default value: 1.0!" + end if + kc = 1.0d+00 + end if + if (kc > kf) then + if (verbose) then + write(*,*) + write(*,wfmt) "'driving_profile_width' must be smaller than " // & + "'driving_wavenumber'!" + write(*,wfmt) "Resetting it to the half of 'driving_wavenumber'!" + end if + kc = 5.0d-01 * kf + end if + if (kl >= kf) then + if (verbose) then + write(*,*) + write(*,wfmt) "'driving_profile_lower' must be smaller than " // & + "'driving_wavenumber'!" + write(*,wfmt) "Setting it to 'driving_wavenumber' " // & + "- 'driving_profile_width'!" + end if + kl = kf - kc + end if + if (ku <= kf) then + if (verbose) then + write(*,*) + write(*,wfmt) "'driving_profile_upper' must be larger than " // & + "'driving_wavenumber'!" + write(*,wfmt) "Setting it to 'driving_wavenumber' " // & + "+ 'driving_profile_width'!" + end if + ku = kf + kc + end if + fmin = max(1.0d-08, fmin) + kl2 = kl**2 + ku2 = ku**2 + +! get the upper bound for the wave number +! + kmax = ceiling(ku) + +! allocate arrays to count the modes with the same wave number +! + allocate(kcount(NDIMS * kmax**2), stat = status) + + if (status == 0) then + +! initialize mode counts +! + kcount(:) = 0 + +! count the modes with the same k, get the minimum sufficient number of +! driving modes +! +#if NDIMS == 3 + do i = 1, kmax + do j = -kmax, kmax + do k = -kmax, kmax + k2 = i**2 + j**2 + k**2 + kv2 = 1.0d+00 * k2 + if (kl2 <= kv2 .and. kv2 <= ku2) then + nmodes = nmodes + 1 + kcount(k2) = kcount(k2) + 2 + end if + end do + end do + end do + do j = 1, kmax + do k = -kmax, kmax + k2 = j**2 + k**2 + kv2 = 1.0d+00 * k2 + if (kl2 <= kv2 .and. kv2 <= ku2) then + nmodes = nmodes + 1 + kcount(k2) = kcount(k2) + 2 + end if + end do + end do + do k = 1, kmax + k2 = k**2 + kv2 = 1.0d+00 * k2 + if (kl2 <= kv2 .and. kv2 <= ku2) then + nmodes = nmodes + 1 + kcount(k2) = kcount(k2) + 2 + end if + end do +#else /* NDIMS == 3 */ + do i = 1, kmax + do j = -kmax, kmax + k2 = i**2 + j**2 + kv2 = 1.0d+00 * k2 + if (kl2 <= kv2 .and. kv2 <= ku2) then + nmodes = nmodes + 1 + kcount(k2) = kcount(k2) + 2 + end if + end do + end do + do j = 1, kmax + k2 = j**2 + kv2 = 1.0d+00 * k2 + if (kl2 <= kv2 .and. kv2 <= ku2) then + nmodes = nmodes + 1 + kcount(k2) = kcount(k2) + 2 + end if + end do +#endif /* NDIMS == 3 */ + +! allocate arrays for driving modes (positions, amplitudes and coefficients) +! + allocate(e1vecs(nmodes,NDIMS), e2vecs(nmodes,NDIMS), & + kmodes(nmodes,NDIMS), fmodes(nmodes), & + fcoefs(nmodes,NDIMS), stat = status) + + if (status == 0) then + +! prepare wave vectors +! + l = 0 +#if NDIMS == 3 + do i = 1, kmax + do j = -kmax, kmax + do k = -kmax, kmax + k2 = i**2 + j**2 + k**2 + kv2 = 1.0d+00 * k2 + if (kl2 <= kv2 .and. kv2 <= ku2) then + l = l + 1 + kmodes(l,1) = i + kmodes(l,2) = j + kmodes(l,3) = k + end if + end do + end do + end do + do j = 1, kmax + do k = -kmax, kmax + k2 = j**2 + k**2 + kv2 = 1.0d+00 * k2 + if (kl2 <= kv2 .and. kv2 <= ku2) then + l = l + 1 + kmodes(l,1) = 0 + kmodes(l,2) = j + kmodes(l,3) = k + end if + end do + end do + do k = 1, kmax + k2 = k**2 + kv2 = 1.0d+00 * k2 + if (kl2 <= kv2 .and. kv2 <= ku2) then + l = l + 1 + kmodes(l,1) = 0 + kmodes(l,2) = 0 + kmodes(l,3) = k + end if + end do +#else /* NDIMS == 3 */ + do i = 1, kmax + do j = -kmax, kmax + k2 = i**2 + j**2 + kv2 = 1.0d+00 * k2 + if (kl2 <= kv2 .and. kv2 <= ku2) then + l = l + 1 + kmodes(l,1) = i + kmodes(l,2) = j + end if + end do + end do + do j = 1, kmax + k2 = j**2 + kv2 = 1.0d+00 * k2 + if (kl2 <= kv2 .and. kv2 <= ku2) then + l = l + 1 + kmodes(l,1) = 0 + kmodes(l,2) = j + end if + end do +#endif /* NDIMS == 3 */ + +! generate random vectors, perpendicular to the wave vector +! + do l = 1, nmodes + + u(:) = kmodes(l,:) +#if NDIMS == 3 + if (abs(kmodes(l,3)) < abs(kmodes(l,1))) then + v(:) = (/ kmodes(l,2), - kmodes(l,1), 0 /) + else + v(:) = (/ 0, - kmodes(l,3), kmodes(l,2) /) + end if +#else /* NDIMS == 3 */ + v(:) = (/ kmodes(l,2), - kmodes(l,1) /) +#endif /* NDIMS == 3 */ + + uu = dot_product(v(:), v(:)) + if (uu > 0.0d+00) then + + e1vecs(l,:) = v(:) +#if NDIMS == 3 + e2vecs(l,1) = u(2) * v(3) - u(3) * v(2) + e2vecs(l,2) = u(3) * v(1) - u(1) * v(3) + e2vecs(l,3) = u(1) * v(2) - u(2) * v(1) +#endif /* NDIMS == 3 */ + + e1vecs(l,:) = e1vecs(l,:) / sqrt(sum(e1vecs(l,:)**2)) +#if NDIMS == 3 + e2vecs(l,:) = e2vecs(l,:) / sqrt(sum(e2vecs(l,:)**2)) +#else /* NDIMS == 3 */ + e2vecs(l,:) = 0.0d+00 +#endif /* NDIMS == 3 */ + + else + write(error_unit,"('[', a, ']: ', a)") trim(loc), "v(:) is ZERO!" + end if + + end do + +! determine if the profile is set by energy or amplitude +! + call get_parameter("profile_by_energy", profile_energy) + select case(trim(profile_energy)) + case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") + energy_profile = .true. + case default + energy_profile = .false. + kcount(:) = 1 + end select + +! initialize driving mode amplitudes depending on the chosen profile +! + call get_parameter("driving_profile", profile_type) + select case(trim(profile_type)) + +! --- power-law profile --- +! + case("power-law", "powerlaw", "band") + + injection_profile = profile_powerlaw + +! get the power index of the power-law distribution +! + call get_parameter("powerlaw_index", pidx) + +! calculate amplitudes of driving modes +! + fi = 0.0d+00 + do l = 1, nmodes + k2 = sum(kmodes(l,:)**2) + kv = sqrt(1.0d+00 * k2) + fa = kv**pidx / kcount(k2) + fi = fi + fa + fmodes(l) = sqrt(fa) + end do + fmodes(:) = fmodes(:) / sqrt(fi) + +! --- normal distribution profile --- +! + case("gauss", "gaussian", "normal") + + injection_profile = profile_normal + +! calculate amplitudes of driving modes +! + fi = 0.0d+00 + do l = 1, nmodes + k2 = sum(kmodes(l,:)**2) + kv = sqrt(1.0d+00 * k2) + fa = exp(- 0.5d+00 * ((kv - kf) / kc)**2) / kcount(k2) + fi = fi + fa + fmodes(l) = sqrt(fa) + end do + fmodes(:) = fmodes(:) / sqrt(fi) + +! --- parabolic profile --- +! + case("parabolic") + + injection_profile = profile_parabolic + +! correct the injection scale and width +! + kf = 0.5d+00 * (ku + kl) + kc = 0.5d+00 * (ku - kl) + +! calculate amplitudes of driving modes +! + fi = 0.0d+00 + do l = 1, nmodes + k2 = sum(kmodes(l,:)**2) + kv = sqrt(1.0d+00 * k2) + fa = 4.0d+00 * (kv - kl) * (ku - kv) / (ku - kl)**2 + fa = max(0.0d+00, fa)**2 / kcount(k2) + fi = fi + fa + fmodes(l) = sqrt(fa) + end do + fmodes(:) = fmodes(:) / sqrt(fi) + +! --- normal distribution profile --- +! + case default + +! warn if the specified profile is not implemented and the normal one is used +! + if (verbose .and. .not. (profile_type == "normal" .or. & + profile_type == "gauss" .or. & + profile_type == "gaussian")) then + write(*,*) + write(*,wfmt) "Unknown spectral profile of driving!" + write(*,wfmt) "Available profiles: " // & + "'power-law', 'normal', 'parabolic'." + write(*,wfmt) "Using default normal distribution profile " // & + "'gaussian'." + + end if + + injection_profile = profile_normal + +! calculate amplitudes of driving modes +! + fi = 0.0d+00 + do l = 1, nmodes + k2 = sum(kmodes(l,:)**2) + kv = sqrt(1.0d+00 * k2) + fa = exp(- 0.5d+00 * ((kv - kf) / kc)**2) / kcount(k2) + fi = fi + fa + fmodes(l) = sqrt(fa) + end do + do l = 1, nmodes + k2 = sum(kmodes(l,:)**2) + fmodes(l) = fmodes(l) / sqrt(fi) + end do + + end select + +! === Ornstein–Uhlenbeck driving === +! + if (injection_method == injection_oh) then + +! get the characteristic driving length, velocity, acceleration scales, and +! determine the estimated driving power (since the main injection is through +! the forcing-velocity correlation, it may not represent the actual +! instantenous injection rate) +! + lscale = 1.0d+00 / kf + vscale = lscale / tscale + ascale = sqrt(2.0d+00) * lscale / tscale**2 / sqrt(tscale) + fpower = vscale**2 / tscale + +! normalize the driving amplitude profile to the driving power +! + fmodes(:) = ascale * fmodes(:) + +! generate the initial driving mode coefficients +! + do l = 1, nmodes + +#if NDIMS == 3 + phi = pi2 * randuni() + fcoefs(l,:) = fmodes(l) * (e1vecs(l,:) * cos(phi) & + + e2vecs(l,:) * sin(phi)) * randnorz() +#else /* NDIMS == 3 */ + fcoefs(l,:) = fmodes(l) * e1vecs(l,:) * randnorz() +#endif /* NDIMS == 3 */ + + end do ! l = 1, nmodes + +! calculate the rms of acceleration +! + arms = sqrt(5.0d-01 * sum(dreal(fcoefs(:,:))**2 & + + dimag(fcoefs(:,:))**2)) + + end if ! injection = 'OH' + +! === Alvelius driving === +! + if (injection_method == injection_alvelius) then + +! get the characteristic driving length, velocity, time, and acceleration scales +! + lscale = 1.0d+00 / kf + tscale = sqrt(lscale / sqrt(2.0d+00 * fpower)) + vscale = lscale / tscale + ascale = vscale / tscale + +! normalize the driving amplitude profile to the driving power +! + fmodes(:) = sqrt(2.0d+00) * ascale * fmodes(:) + +! reset the driving coefficients and calculate the rms of acceleration +! + fcoefs(:,:) = cmplx(0.0d+00, 0.0d+00) + arms = 0.0d+00 + +! allocate vectors of velocity Fourier coefficients +! + allocate(vcoefs(nmodes,NDIMS), stat = status) + if (status /= 0) then + write(error_unit,"('[', a, ']: ', a)") trim(loc), & + "Cannot allocate vcoefs(:,:)!" + end if + + end if ! injection = 'Alvelius' + + end if ! allocation of driving mode arrays + end if ! allocation of kcount + +! deallocate mode counter +! + if (allocated(kcount)) deallocate(kcount) + + end if ! Fourier forcing + +#ifdef PROFILE +! stop accounting time +! + call stop_timer(ifi) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine initialize_forcing +! +!=============================================================================== +! +! subroutine FINALIZE_FORCING: +! --------------------------- +! +! Subroutine releases memory used by the module variables. +! +! Arguments: +! +! status - return flag of the procedure execution status; +! +!=============================================================================== +! + subroutine finalize_forcing(status) + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + integer, intent(out) :: status +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for initialization/finalization +! + call start_timer(ifi) +#endif /* PROFILE */ + + if (allocated(kmodes)) deallocate(kmodes) + if (allocated(fmodes)) deallocate(fmodes) + if (allocated(fcoefs)) deallocate(fcoefs) + if (allocated(e1vecs)) deallocate(e1vecs) + if (allocated(e2vecs)) deallocate(e2vecs) + if (allocated(vcoefs)) deallocate(vcoefs) + +#ifdef PROFILE +! stop accounting time +! + call stop_timer(ifi) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine finalize_forcing +! +!=============================================================================== +! +! subroutine PRINT_FORCING: +! ------------------------ +! +! Subroutine prints module parameters. +! +! Arguments: +! +! verbose - a logical flag turning the information printing; +! +!=============================================================================== +! + subroutine print_forcing(verbose) + +! import external procedures and variables +! + use helpers , only : print_section, print_parameter + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + logical, intent(in) :: verbose + +! local variables +! + character(len=64) :: injection = "none" + character(len=64) :: profile_type = "gaussian" +! +!------------------------------------------------------------------------------- +! + if (verbose) then + + select case(injection_method) + case(injection_eddy) + injection = "random eddies" + case(injection_oh) + injection = "Ornstein–Uhlenbeck" + case(injection_alvelius) + injection = "Alvelius" + case default + injection = "none" + end select + + call print_section(verbose, "Forcing") + call print_parameter(verbose, "method", injection) + + select case(injection_method) + case(injection_eddy) + call print_parameter(verbose, "injection power", power) + call print_parameter(verbose, "injection rate" , rate) + call print_parameter(verbose, "eddy amplitude" , amp) + call print_parameter(verbose, "eddy width" , del) + end select + + if (forcing_fourier) then + if (injection_method == injection_oh) then + call print_parameter(verbose, "power (estimated)", fpower) + end if + if (injection_method == injection_alvelius) then + call print_parameter(verbose, "power" , fpower) + end if + call print_parameter(verbose, "length scale" , lscale) + call print_parameter(verbose, "time scale" , tscale) + call print_parameter(verbose, "velocity scale" , vscale) + call print_parameter(verbose, "acceleration scale", ascale) + select case(injection_profile) + case(profile_powerlaw) + profile_type = "power-law" + case(profile_normal) + profile_type = "gaussian" + case(profile_parabolic) + profile_type = "parabolic" + case default + profile_type = "unknown" + end select + + call print_parameter(verbose, "spectrum profile", profile_type) + if (energy_profile) then + call print_parameter(verbose, "profile of", "energy") + else + call print_parameter(verbose, "profile of", "amplitude") + end if + call print_parameter(verbose, "number of driving modes", nmodes) + call print_parameter(verbose, "driving wavenumber" , kf) + call print_parameter(verbose, "lower wavenumber limit" , kl) + call print_parameter(verbose, "upper wavenumber limit" , ku) + if (injection_profile == profile_powerlaw) then + call print_parameter(verbose, "power-law index", pidx) + end if + if (injection_profile == profile_normal) then + call print_parameter(verbose, "spectrum profile width", kc) + end if + call print_parameter(verbose, "amplitude threshold", fmin) + end if + + end if + +!------------------------------------------------------------------------------- +! + end subroutine print_forcing +! +!=============================================================================== +! +! subroutine UPDATE_FORCING_EDDIES: +! -------------------------------- +! +! Subroutine adds the energy injection terms. +! +! Arguments: +! +! t, dt - time and its increment; +! +!=============================================================================== +! + subroutine update_forcing_eddies(t, dt) + +! import external procedures and variables +! + use coordinates, only : xmin, ymin, zmin, xlen, ylen, zlen + use random , only : randuni, randsym + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8), intent(in) :: t, dt + +! local variables +! + integer :: ni, n + real(kind=8) :: tmp + real(kind=8), dimension(3) :: xp, ap +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for forcing term update +! + call start_timer(ifu) +#endif /* PROFILE */ + +! calculate the number of eddies to be injected during this interval +! + dinj = dinj + rate * dt + ni = int(floor(dinj)) + dinj = dinj - 1.0d+00 * ni + +! inject the required number of eddies +! + if (ni > 0) then + + do n = 1, ni + +! get random position +! + xp(1) = xmin + xlen * randuni() + xp(2) = ymin + ylen * randuni() + xp(3) = zmin + zlen * randuni() + +! get random orientation +! +#if NDIMS == 3 + tmp = 0.0d+00 + do while(tmp == 0.0d+00) + ap(1) = randsym() + ap(2) = randsym() + ap(3) = randsym() + tmp = sqrt(ap(1)**2 + ap(2)**2 + ap(3)**2) + end do + ap(:) = amp * (ap(:) / tmp) / del +#else /* NDIMS == 3 */ + ap(:) = sign(1.0d+00, randsym()) * (/ 0.0d+00, 0.0d+00, amp / del /) +#endif /* NDIMS == 3 */ + ap(:) = ap(:) * dt + +! iterate over data blocks and add forcing components +! + call inject_eddy(xp(:), ap(:)) + + end do + + end if + +#ifdef PROFILE +! stop accounting time +! + call stop_timer(ifu) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine update_forcing_eddies +! +!=============================================================================== +! +! subroutine UPDATE_FORCING_OH: +! ---------------------------- +! +! Subroutine adds the energy injection terms for Ornstein–Uhlenbeck driving. +! +! Arguments: +! +! t, dt - time and its increment; +! +!=============================================================================== +! + subroutine update_forcing_oh(t, dt) + +! import external procedures and variables +! + use constants, only : pi2 + use random , only : randuni, randnorz + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8), intent(in) :: t, dt + +! local variables +! + integer :: l + real(kind=8) :: acoeff, dcoeff, phi + real(kind=8) :: dinj + +! local vectors +! + complex(kind=8), dimension(NDIMS) :: finj +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for forcing term update +! + call start_timer(ifu) +#endif /* PROFILE */ + +! calculate drift and diffusion coefficients +! + acoeff = exp(- dt / tscale) + dcoeff = sqrt(1.0d+00 - acoeff**2) + +! iterate over driving modes +! + do l = 1, nmodes + +! generate gaussian random vector along the driving modes +! +#if NDIMS == 3 + phi = pi2 * randuni() + finj(:) = fmodes(l) * (e1vecs(l,:) * cos(phi) & + + e2vecs(l,:) * sin(phi)) * randnorz() +#else /* NDIMS == 3 */ + finj(:) = fmodes(l) * e1vecs(l,:) * randnorz() +#endif /* NDIMS == 3 */ + +! advance the driving force modes +! + fcoefs(l,:) = acoeff * fcoefs(l,:) + dcoeff * finj(:) + + end do + +! calculate the rms of acceleration +! + arms = sqrt(5.0d-01 * sum(dreal(fcoefs(:,:))**2 + dimag(fcoefs(:,:))**2)) + +! store previously injected energy +! + dinj = einj + +! inject driving modes +! + call inject_fmodes(dt) + +! calculate the energy injected during this step +! + dinj = einj - dinj + +! calculate the injection rate +! + rinj = dinj / dt + +#ifdef PROFILE +! stop accounting time +! + call stop_timer(ifu) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine update_forcing_oh +! +!=============================================================================== +! +! subroutine UPDATE_FORCING_ALVELIUS: +! ---------------------------------- +! +! Subroutine adds the energy injection terms for Alvelius driving. +! +! Arguments: +! +! t, dt - time and its increment; +! +!=============================================================================== +! + subroutine update_forcing_alvelius(t, dt) + +! import external procedures and variables +! + use constants, only : pi2 + use random , only : randuni + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8), intent(in) :: t, dt + +! local variables +! + integer :: l + real(kind=8) :: th1, th2, phi, psi, ga, gb, dinj, sqdt + complex(kind=8) :: aran, bran, xi1, xi2 +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for forcing term update +! + call start_timer(ifu) +#endif /* PROFILE */ + +! determine velocify coefficients in Fourier space +! + call get_vcoefs() + +! square root of dt +! + sqdt = sqrt(dt) + +! iterate over driving modes +! + do l = 1, nmodes + +#if NDIMS == 3 +! get xi1 and xi2 +! + xi1 = dot_product(vcoefs(l,:), e1vecs(l,:)) + xi2 = dot_product(vcoefs(l,:), e2vecs(l,:)) + +! get random angles phi and psi and determine th1 and th2 to reduce +! velocity-force correlation +! + phi = pi2 * randuni() + ga = sin(phi) + gb = cos(phi) + psi = pi2 * randuni() + th1 = - ga * dimag(xi1) + gb * (sin(psi) * dreal(xi2) & + - cos(psi) * dimag(xi2)) + if (th1 /= 0.0d+00) then + th1 = atan((ga * dreal(xi1) + gb * (sin(psi) * dimag(xi2) & + + cos(psi) * dreal(xi2))) / th1) + end if + th2 = psi + th1 + +! calculate amplitudes +! + aran = cmplx(cos(th1), sin(th1)) * ga + bran = cmplx(cos(th2), sin(th2)) * gb + +! calculate driving mode coefficient +! + fcoefs(l,:) = fmodes(l) * (aran * e1vecs(l,:) & + + bran * e2vecs(l,:)) / sqdt +#else /* NDIMS == 3 */ +! get random phase +! + th1 = pi2 * randuni() + aran = cmplx(cos(th1), sin(th1)) + +! calculate driving mode coefficient +! + fcoefs(l,:) = fmodes(l) * aran * e1vecs(l,:) / sqdt +#endif /* NDIMS == 3 */ + + end do + +! calculate the rms of acceleration +! + arms = sqrt(5.0d-01 * sum(dreal(fcoefs(:,:))**2 + dimag(fcoefs(:,:))**2)) + +! store previously injected energy +! + dinj = einj + +! inject driving modes +! + call inject_fmodes(dt) + +! calculate the energy injected during this step +! + dinj = einj - dinj + +! calculate the injection rate +! + rinj = dinj / dt + +#ifdef PROFILE +! stop accounting time +! + call stop_timer(ifu) +#endif /* PROFILE */ + +!------------------------------------------------------------------------------- +! + end subroutine update_forcing_alvelius +! +!=============================================================================== +!! +!!*** PRIVATE SUBROUTINES **************************************************** +!! +!=============================================================================== +! +!=============================================================================== +! +! subroutine INJECT_EDDY: +! ---------------------- +! +! Subroutine injects one eddy to the domain at given position. +! +! Arguments: +! +! xp, ap - eddy position and amplitude vectors; +! +!=============================================================================== +! + subroutine inject_eddy(xp, ap) + +! include external variables +! + use blocks, only : block_data, list_data + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8), dimension(3), intent(in) :: xp, ap + +! local pointers +! + type(block_data), pointer :: pdata +! +!------------------------------------------------------------------------------- +! +! assign pdata with the first block on the data block list +! + pdata => list_data + +! iterate over all data blocks +! + do while (associated(pdata)) + +! inject eddy into the current block +! + call inject_eddy_block(pdata, xp(:), ap(:)) + +! assign pdata to the next block +! + pdata => pdata%next + + end do ! over data blocks + +!------------------------------------------------------------------------------- +! + end subroutine inject_eddy +! +!=============================================================================== +! +! subroutine INJECT_EDDY_BLOCK: +! ---------------------------- +! +! Subroutine injects one eddy to the local block. +! +! Arguments: +! +! xp, ap - eddy position and amplitude vectors; +! +!=============================================================================== +! + subroutine inject_eddy_block(pdata, xp, ap) + +! include external variables +! + use blocks , only : block_data + use coordinates, only : nm => bcells + use coordinates, only : ax, ay, az + use coordinates, only : xlen, ylen, zlen + use coordinates, only : periodic + use equations , only : nv + use equations , only : idn, imx, imy, imz, ien + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + type(block_data), pointer , intent(inout) :: pdata + real(kind=8), dimension(3), intent(in) :: xp, ap + +! local variables +! + integer :: i, j, k = 1 + real(kind=8) :: x2, y2, z2, r2 + real(kind=8) :: fx, fy, fz, fp, e1, e2 + +! local arrays +! + real(kind=8), dimension(nm) :: x, y, z +! +!------------------------------------------------------------------------------- +! +! prepare block coordinates +! + if (periodic(1)) then + fx = 0.5d+00 * xlen + x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) - xp(1) + x(:) = abs(x(:) / fx) + x(:) = (abs(1.0d+00 - abs(1.0d+00 - abs(x(:))))) * fx / del + else + x(:) = (pdata%meta%xmin + ax(pdata%meta%level,:) - xp(1)) / del + end if + if (periodic(2)) then + fy = 0.5d+00 * ylen + y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) - xp(2) + y(:) = abs(y(:) / fy) + y(:) = (abs(1.0d+00 - abs(1.0d+00 - abs(y(:))))) * fy / del + else + y(:) = (pdata%meta%ymin + ay(pdata%meta%level,:) - xp(2)) / del + end if +#if NDIMS == 3 + if (periodic(3)) then + fz = 0.5d+00 * zlen + z(:) = pdata%meta%zmin + az(pdata%meta%level,:) - xp(3) + z(:) = abs(z(:) / fz) + z(:) = (abs(1.0d+00 - abs(1.0d+00 - abs(z(:))))) * fz / del + else + z(:) = (pdata%meta%zmin + az(pdata%meta%level,:) - xp(3)) / del + end if +#else /* NDIMS == 3 */ + z(:) = 0.0d+00 +#endif /* NDIMS == 3 */ + +! iterate over the block coordinates +! + if (ien > 0) then +#if NDIMS == 3 + do k = 1, nm + z2 = z(k)**2 +#endif /* NDIMS == 3 */ + do j = 1, nm + y2 = y(j)**2 + do i = 1, nm + x2 = x(i)**2 +#if NDIMS == 3 + r2 = x2 + y2 + z2 +#else /* NDIMS == 3 */ + r2 = x2 + y2 +#endif /* NDIMS == 3 */ + + fp = pdata%u(idn,i,j,k) * exp(- 0.5d+00 * r2) + +#if NDIMS == 3 + fx = (- ap(3) * y(j) + ap(2) * z(k)) * fp + fy = (- ap(1) * z(k) + ap(3) * x(i)) * fp + fz = (- ap(2) * x(i) + ap(1) * y(j)) * fp +#else /* NDIMS == 3 */ + fx = (- ap(3) * y(j)) * fp + fy = (+ ap(3) * x(i)) * fp + fz = 0.0d+00 +#endif /* NDIMS == 3 */ + + e1 = 0.5d+00 * sum(pdata%u(imx:imx,i,j,k)**2) / pdata%u(idn,i,j,k) + + pdata%u(imx,i,j,k) = pdata%u(imx,i,j,k) + fx + pdata%u(imy,i,j,k) = pdata%u(imy,i,j,k) + fy + pdata%u(imz,i,j,k) = pdata%u(imz,i,j,k) + fz + + e2 = 0.5d+00 * sum(pdata%u(imx:imz,i,j,k)**2) / pdata%u(idn,i,j,k) + + pdata%u(ien,i,j,k) = pdata%u(ien,i,j,k) + (e2 - e1) + + end do + end do +#if NDIMS == 3 + end do +#endif /* NDIMS == 3 */ + else +#if NDIMS == 3 + do k = 1, nm + z2 = z(k)**2 +#endif /* NDIMS == 3 */ + do j = 1, nm + y2 = y(j)**2 + do i = 1, nm + x2 = x(i)**2 +#if NDIMS == 3 + r2 = x2 + y2 + z2 +#else /* NDIMS == 3 */ + r2 = x2 + y2 +#endif /* NDIMS == 3 */ + + fp = pdata%u(idn,i,j,k) * exp(- 0.5d+00 * r2) + +#if NDIMS == 3 + fx = (- ap(3) * y(j) + ap(2) * z(k)) * fp + fy = (- ap(1) * z(k) + ap(3) * x(i)) * fp + fz = (- ap(2) * x(i) + ap(1) * y(j)) * fp +#else /* NDIMS == 3 */ + fx = (- ap(3) * y(j)) * fp + fy = (+ ap(3) * x(i)) * fp + fz = 0.0d+00 +#endif /* NDIMS == 3 */ + + pdata%u(imx,i,j,k) = pdata%u(imx,i,j,k) + fx + pdata%u(imy,i,j,k) = pdata%u(imy,i,j,k) + fy + pdata%u(imz,i,j,k) = pdata%u(imz,i,j,k) + fz + + end do + end do +#if NDIMS == 3 + end do +#endif /* NDIMS == 3 */ + end if + +!------------------------------------------------------------------------------- +! + end subroutine inject_eddy_block +! +!=============================================================================== +! +! subroutine INJECT_FMODES: +! ------------------------ +! +! Subroutine injects Fourier driving modes. +! +! Arguments: +! +! dt - the time increment; +! +!=============================================================================== +! + subroutine inject_fmodes(dt) + +! include external variables +! + use blocks, only : block_data, list_data + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + real(kind=8), intent(in) :: dt + +! local pointers +! + type(block_data), pointer :: pdata +! +!------------------------------------------------------------------------------- +! +! assign pdata with the first block on the data block list +! + pdata => list_data + +! iterate over all data blocks +! + do while (associated(pdata)) + +! inject eddy into the current block +! + call inject_fmodes_block(pdata, dt) + +! assign pdata to the next block +! + pdata => pdata%next + + end do ! over data blocks + +!------------------------------------------------------------------------------- +! + end subroutine inject_fmodes +! +!=============================================================================== +! +! subroutine INJECT_FMODES_BLOCK: +! ------------------------------ +! +! Subroutine injects Fourier driving modes to the local block. +! +! Arguments: +! +! pdata - a pointer to the data block; +! dt - the time increment; +! +!=============================================================================== +! + subroutine inject_fmodes_block(pdata, dt) + +! include external variables +! + use blocks , only : block_data + use constants , only : pi2 + use coordinates, only : nm => bcells, nb, ne + use coordinates, only : ax, ay, az, advol + use equations , only : idn, imx, imy, imz, ien, ivx, ivy, ivz + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + type(block_data), pointer, intent(inout) :: pdata + real(kind=8) , intent(in) :: dt + +! local variables +! + integer :: i, j, k = 1, l, n + real(kind=8) :: cs, sn, tt + real(kind=8) :: dvol + +! local arrays +! + real(kind=8), dimension(nm):: x, y, z, kx, ky, kz + real(kind=8), dimension(nm):: snkx, snky, snkz + real(kind=8), dimension(nm):: cskx, csky, cskz +#if NDIMS == 3 + real(kind=8), dimension(3,nm,nm,nm) :: acc + real(kind=8), dimension( nm,nm,nm) :: den +#else /* NDIMS == 3 */ + real(kind=8), dimension(2,nm,nm, 1) :: acc + real(kind=8), dimension( nm,nm, 1) :: den +#endif /* NDIMS == 3 */ +! +!------------------------------------------------------------------------------- +! +! prepare block coordinates +! + x(:) = - pi2 * (pdata%meta%xmin + ax(pdata%meta%level,:)) + y(:) = - pi2 * (pdata%meta%ymin + ay(pdata%meta%level,:)) +#if NDIMS == 3 + z(:) = - pi2 * (pdata%meta%zmin + az(pdata%meta%level,:)) +#else /* NDIMS == 3 */ + z(:) = 0.0d+00 +#endif /* NDIMS == 3 */ + dvol = advol(pdata%meta%level) + +! reset the acceleration +! + acc(:,:,:,:) = 0.0d+00 + +! iterate over driving modes and compute the acceleration in the real space +! + do l = 1, nmodes + + if (fmodes(l) > fmin) then + + kx(:) = kmodes(l,1) * x(:) + ky(:) = kmodes(l,2) * y(:) +#if NDIMS == 3 + kz(:) = kmodes(l,3) * z(:) +#endif /* NDIMS == 3 */ + cskx(:) = cos(kx(:)) + snkx(:) = sin(kx(:)) + csky(:) = cos(ky(:)) + snky(:) = sin(ky(:)) +#if NDIMS == 3 + cskz(:) = cos(kz(:)) + snkz(:) = sin(kz(:)) +#endif /* NDIMS == 3 */ + +#if NDIMS == 3 + do k = 1, nm +#endif /* NDIMS == 3 */ + do j = 1, nm + do i = 1, nm + + cs = cskx(i) * csky(j) - snkx(i) * snky(j) + sn = cskx(i) * snky(j) + snkx(i) * csky(j) +#if NDIMS == 3 + tt = cs + cs = tt * cskz(k) - sn * snkz(k) + sn = tt * snkz(k) + sn * cskz(k) +#endif /* NDIMS == 3 */ + + acc(:,i,j,k) = acc(:,i,j,k) + (dreal(fcoefs(l,:)) * cs & + + dimag(fcoefs(l,:)) * sn) + end do + end do +#if NDIMS == 3 + end do +#endif /* NDIMS == 3 */ + + end if + + end do ! l = 1, nmodes + +! multiply the acceleration by density and time step +! + do n = 1, NDIMS + acc(n,:,:,:) = pdata%u(idn,:,:,:) * acc(n,:,:,:) * dt + end do + +! calculate the kinetic energy before the injection +! + den(:,:,:) = sum(pdata%u(imx:imz,:,:,:)**2, 1) + +! add the momentum increment +! + pdata%u(imx,:,:,:) = pdata%u(imx,:,:,:) + acc(1,:,:,:) + pdata%u(imy,:,:,:) = pdata%u(imy,:,:,:) + acc(2,:,:,:) +#if NDIMS == 3 + pdata%u(imz,:,:,:) = pdata%u(imz,:,:,:) + acc(3,:,:,:) +#endif /* NDIMS == 3 */ + +! calculate the injected energy +! + den(:,:,:) = 5.0d-01 * (sum(pdata%u(imx:imz,:,:,:)**2, 1) - den(:,:,:)) & + / pdata%u(idn,:,:,:) + +! calculate the injected energy +! +#if NDIMS == 3 + einj = einj + sum(den(nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + einj = einj + sum(den(nb:ne,nb:ne, 1 )) * dvol +#endif /* NDIMS == 3 */ + +! add the energy increment +! + if (ien > 0) then + pdata%u(ien,:,:,:) = pdata%u(ien,:,:,:) + den(:,:,:) + end if + +!------------------------------------------------------------------------------- +! + end subroutine inject_fmodes_block +! +!=============================================================================== +! +! subroutine GET_VCOEFS: +! --------------------- +! +! Subroutine determines velocity Fourier coefficients. +! +! +!=============================================================================== +! + subroutine get_vcoefs() + +! include external variables +! + use blocks , only : block_data, list_data +#ifdef MPI + use mpitools, only : reduce_sum_complex_array +#endif /* MPI */ + +! local variables are not implicit by default +! + implicit none + +! local variables +! + integer :: status + +! local pointers +! + type(block_data), pointer :: pdata +! +!------------------------------------------------------------------------------- +! +! reset vcoefs +! + vcoefs(:,:) = cmplx(0.0d+00, 0.0d+00) + +! assign pdata with the first block on the data block list +! + pdata => list_data + +! iterate over all data blocks +! + do while (associated(pdata)) + +! get contribution of velocity coefficients from the current block +! + call get_vcoefs_block(pdata) + +! assign pdata to the next block +! + pdata => pdata%next + + end do ! over data blocks + +#ifdef MPI +! reduce velocity coefficients over all processes +! + call reduce_sum_complex_array(size(vcoefs), vcoefs, status) +#endif /* MPI */ + +!------------------------------------------------------------------------------- +! + end subroutine get_vcoefs +! +!=============================================================================== +! +! subroutine GET_VCOEFS_BLOCK: +! --------------------------- +! +! Subroutine determines Fourier coefficients of velocity. +! +! Arguments: +! +! pdata - a pointer to the data block; +! +!=============================================================================== +! + subroutine get_vcoefs_block(pdata) + +! include external variables +! + use blocks , only : block_data + use constants , only : pi2 + use coordinates, only : nm => bcells, nb, ne + use coordinates, only : ax, ay, az, advol + use equations , only : ivx, ivy, ivz + +! local variables are not implicit by default +! + implicit none + +! subroutine arguments +! + type(block_data), pointer, intent(inout) :: pdata + +! local variables +! + integer :: i, j, k = 1, l + real(kind=8) :: cs, sn, tt, dvol + complex(kind=8) :: cf + +! local arrays +! + real(kind=8), dimension(nm):: x, y, z, kx, ky, kz + real(kind=8), dimension(nm):: snkx, snky, snkz + real(kind=8), dimension(nm):: cskx, csky, cskz +! +!------------------------------------------------------------------------------- +! +! prepare block coordinates +! + x(:) = pi2 * (pdata%meta%xmin + ax(pdata%meta%level,:)) + y(:) = pi2 * (pdata%meta%ymin + ay(pdata%meta%level,:)) +#if NDIMS == 3 + z(:) = pi2 * (pdata%meta%zmin + az(pdata%meta%level,:)) +#else /* NDIMS == 3 */ + z(:) = 0.0d+00 +#endif /* NDIMS == 3 */ + dvol = advol(pdata%meta%level) + +! iterate over driving modes and compute the velocity Fourier coefficients +! + do l = 1, nmodes + + kx(:) = kmodes(l,1) * x(:) + ky(:) = kmodes(l,2) * y(:) +#if NDIMS == 3 + kz(:) = kmodes(l,3) * z(:) +#endif /* NDIMS == 3 */ + cskx(:) = cos(kx(:)) + snkx(:) = sin(kx(:)) + csky(:) = cos(ky(:)) + snky(:) = sin(ky(:)) +#if NDIMS == 3 + cskz(:) = cos(kz(:)) + snkz(:) = sin(kz(:)) +#endif /* NDIMS == 3 */ + +#if NDIMS == 3 + do k = nb, ne +#endif /* NDIMS == 3 */ + do j = nb, ne + do i = nb, ne + + cs = cskx(i) * csky(j) - snkx(i) * snky(j) + sn = cskx(i) * snky(j) + snkx(i) * csky(j) +#if NDIMS == 3 + tt = cs + cs = tt * cskz(k) - sn * snkz(k) + sn = tt * snkz(k) + sn * cskz(k) +#endif /* NDIMS == 3 */ + + cf = cmplx(cs, sn) * dvol + + vcoefs(l,1) = vcoefs(l,1) + pdata%q(ivx,i,j,k) * cf + vcoefs(l,2) = vcoefs(l,2) + pdata%q(ivy,i,j,k) * cf +#if NDIMS == 3 + vcoefs(l,3) = vcoefs(l,3) + pdata%q(ivz,i,j,k) * cf +#endif /* NDIMS == 3 */ + + end do + end do +#if NDIMS == 3 + end do +#endif /* NDIMS == 3 */ + + end do ! l = 1, nmodes + +!------------------------------------------------------------------------------- +! + end subroutine get_vcoefs_block + +!=============================================================================== +! +end module forcing diff --git a/sources/hash.F90 b/sources/hash.F90 new file mode 100644 index 0000000..c5e5a24 --- /dev/null +++ b/sources/hash.F90 @@ -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 +!! +!! 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 . +!! +!!****************************************************************************** +!! +!! 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 diff --git a/sources/helpers.F90 b/sources/helpers.F90 index 497b9de..d9e4429 100644 --- a/sources/helpers.F90 +++ b/sources/helpers.F90 @@ -215,9 +215,9 @@ module helpers msg = trim(adjustl(description)) if (value >= 0.0d+00) then - write(*,"(4x,a26,1x,'=',es9.2)") msg, value + write(*,"(4x,a26,1x,'=',es12.5)") msg, value else - write(*,"(4x,a26,1x,'=',1x,es9.2)") msg, value + write(*,"(4x,a26,1x,'=',1x,es12.5)") msg, value end if !------------------------------------------------------------------------------- diff --git a/sources/integrals.F90 b/sources/integrals.F90 index 1668e26..b84f681 100644 --- a/sources/integrals.F90 +++ b/sources/integrals.F90 @@ -180,9 +180,10 @@ module integrals ! 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' & - , 'ener', 'ekin', 'emag', 'eint' + , 'ener', 'ekin', 'emag', 'eint' & + , 'einj', 'erat', 'arms' write(funit,"('#')") ! 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 : magnetized, gamma, csnd use evolution , only : step, time, dtn + use forcing , only : einj, rinj, arms use mpitools , only : master #ifdef MPI use mpitools , only : reduce_sum_real_array @@ -518,6 +520,11 @@ module integrals #endif /* NDIMS == 3 */ end if +! sum up the injected energy and injection rate +! + inarr( 9) = einj + inarr(10) = rinj + ! get average, minimum and maximum values of density ! #if NDIMS == 3 @@ -881,7 +888,7 @@ module integrals ! write down the integrals and statistics to appropriate files ! 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 & , avarr(1), mnarr(1), mxarr(1) & , avarr(2), mnarr(2), mxarr(2) & diff --git a/sources/io.F90 b/sources/io.F90 index 0b8447d..6ffa4bc 100644 --- a/sources/io.F90 +++ b/sources/io.F90 @@ -45,13 +45,23 @@ module io ! subroutine interfaces ! interface read_snapshot_parameter + module procedure read_snapshot_parameter_string + module procedure read_snapshot_parameter_integer + module procedure read_snapshot_parameter_double + end interface + interface write_attribute_xml + module procedure write_attribute_xml_string + module procedure write_attribute_xml_integer + module procedure write_attribute_xml_double + end interface #ifdef HDF5 + interface read_snapshot_parameter_h5 module procedure read_snapshot_parameter_string_h5 module procedure read_snapshot_parameter_integer_h5 module procedure read_snapshot_parameter_integer_vector_h5 module procedure read_snapshot_parameter_double_h5 -#endif /* HDF5 */ end interface +#endif /* HDF5 */ interface write_attribute #ifdef HDF5 module procedure write_scalar_attribute_string_h5 @@ -59,6 +69,8 @@ module io module procedure write_scalar_attribute_double_h5 module procedure write_vector_attribute_integer_h5 module procedure write_vector_attribute_double_h5 + module procedure write_array_attribute_long_h5 + module procedure write_array_attribute_complex_h5 #endif /* HDF5 */ end interface interface read_attribute @@ -67,6 +79,8 @@ module io module procedure read_scalar_attribute_double_h5 module procedure read_vector_attribute_integer_h5 module procedure read_vector_attribute_double_h5 + module procedure read_array_attribute_long_h5 + module procedure read_array_attribute_complex_h5 #endif /* HDF5 */ end interface interface write_array @@ -108,6 +122,19 @@ module io ! MODULE PARAMETERS: ! ================= ! +! snapshot formats +! + integer, parameter :: snapshot_xml = 0 +#ifdef HDF5 + integer, parameter :: snapshot_hdf5 = 1 +#endif /* HDF5 */ + +! snapshot_format - the format of snapshots; +! restart_format - the format of restart snapshots; +! + integer, save :: snapshot_format = snapshot_xml + integer, save :: restart_format = snapshot_xml + ! respath - the directory from which the restart snapshots should be read; ! ftype - the type of snapshots to write: ! 'p' -> all primitive variables (default); @@ -231,6 +258,7 @@ module io ! local variables ! + character(len=255) :: sformat = "xml" character(len=255) :: precise = "off" character(len=255) :: ghosts = "on" character(len=255) :: xdmf = "off" @@ -274,6 +302,36 @@ module io call get_parameter("include_ghosts" , ghosts ) call get_parameter("generate_xdmf" , xdmf ) +! add slash at the end of respath if not present +! + if (index(respath, '/', back = .true.) /= len(trim(respath))) then + write(respath,"(a)") trim(adjustl(respath)) // '/' + end if + +! check the snapshot format +! + call get_parameter("snapshot_format" , sformat) + select case(sformat) +#ifdef HDF5 + case('h5', 'hdf5', 'H5', 'HDF5') + snapshot_format = snapshot_hdf5 +#endif /* HDF5 */ + case default + snapshot_format = snapshot_xml + end select + +! check the restart snapshot format +! + call get_parameter("restart_format" , sformat) + select case(sformat) +#ifdef HDF5 + case('h5', 'hdf5', 'H5', 'HDF5') + restart_format = snapshot_hdf5 +#endif /* HDF5 */ + case default + restart_format = snapshot_xml + end select + ! check the snapshot type ! select case(ftype) @@ -508,6 +566,22 @@ module io ! if (verbose) then call print_section(verbose, "Snapshots") + select case(snapshot_format) +#ifdef HDF5 + case(snapshot_hdf5) + call print_parameter(verbose, "snapshot format", "HDF5") +#endif /* HDF5 */ + case default + call print_parameter(verbose, "snapshot format", "XML+binary") + end select + select case(restart_format) +#ifdef HDF5 + case(snapshot_hdf5) + call print_parameter(verbose, "restart snapshot format", "HDF5") +#endif /* HDF5 */ + case default + call print_parameter(verbose, "restart snapshot format", "XML+binary") + end select if (precise_snapshots) then call print_parameter(verbose, "precise snapshot intervals", "on" ) else @@ -614,58 +688,43 @@ module io ! ! Arguments: ! -! iret - the return flag to inform if subroutine succeeded or failed; +! status - the status flag to inform if subroutine succeeded or failed; ! !=============================================================================== ! - subroutine read_restart_snapshot(iret) + subroutine read_restart_snapshot(status) -! import external variables -! - use evolution , only : time + use evolution, only : time -! local variables are not implicit by default -! implicit none -! input and output arguments -! - integer, intent(out) :: iret + integer, intent(out) :: status ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for the data reading -! call start_timer(ios) #endif /* PROFILE */ - -! reset the return flag -! - iret = 0 - -! start accounting time for I/O -! call start_timer(iio) -#ifdef HDF5 -! read HDF5 restart file and rebuild the meta and data block structures -! - call read_restart_snapshot_h5(iret) -#endif /* HDF5 */ + status = 0 -! stop accounting time for I/O -! - call stop_timer(iio) + select case(restart_format) +#ifdef HDF5 + case(snapshot_hdf5) + call read_restart_snapshot_h5(status) +#endif /* HDF5 */ + case default + call read_restart_snapshot_xml(status) + end select ! calculate the shift of the snapshot counter, and the next snapshot time ! ishift = int(time / hsnap) - isnap + 1 tsnap = (ishift + isnap) * hsnap + call stop_timer(iio) #ifdef PROFILE -! stop accounting time for the data reading -! call stop_timer(ios) #endif /* PROFILE */ @@ -675,21 +734,21 @@ module io ! !=============================================================================== ! -! subroutine WRITE_RESTART_SNAPSHOTS: -! ---------------------------------- +! subroutine WRITE_RESTART_SNAPSHOT: +! --------------------------------- ! ! Subroutine stores current restart snapshot files. This is a wrapper ! calling specific format subroutine. ! ! Arguments: ! -! thrs - the current execution time in hours; -! nrun - the run number; -! iret - the return flag; +! thrs - the current execution time in hours; +! nrun - the run number; +! status - the status flag; ! !=============================================================================== ! - subroutine write_restart_snapshot(thrs, nrun, iret) + subroutine write_restart_snapshot(thrs, nrun, status) ! local variables are not implicit by default ! @@ -699,45 +758,36 @@ module io ! real(kind=8), intent(in) :: thrs integer , intent(in) :: nrun - integer , intent(out) :: iret + integer , intent(out) :: status ! !------------------------------------------------------------------------------- ! -! reset the return flag -! - iret = 0 + status = 0 ! check if conditions for storing the restart snapshot have been met ! if (hrest < 5.0d-02 .or. thrs < irest * hrest) return #ifdef PROFILE -! start accounting time for the data writing -! call start_timer(iow) #endif /* PROFILE */ - -! start accounting time for I/O -! call start_timer(iio) + select case(snapshot_format) #ifdef HDF5 -! store restart file -! - call write_restart_snapshot_h5(nrun, iret) + case(snapshot_hdf5) + call write_restart_snapshot_h5(nrun, status) #endif /* HDF5 */ - -! stop accounting time for I/O -! - call stop_timer(iio) + case default + call write_restart_snapshot_xml(nrun, status) + end select ! increase the restart snapshot counter ! irest = irest + 1 + call stop_timer(iio) #ifdef PROFILE -! stop accounting time for the data writing -! call stop_timer(iow) #endif /* PROFILE */ @@ -759,53 +809,44 @@ module io ! subroutine write_snapshot() -! import external variables -! use evolution , only : time use mpitools , only : master -! local variables are not implicit by default -! implicit none + +! local variables +! + integer :: status ! !------------------------------------------------------------------------------- -! -! check if conditions for storing the regular snapshot have been met ! if (hsnap <= 0.0e+00 .or. time < tsnap) return #ifdef PROFILE -! start accounting time for the data writing -! call start_timer(iow) #endif /* PROFILE */ - -! start accounting time for I/O -! call start_timer(iio) + select case(snapshot_format) #ifdef HDF5 -! store variable snapshot file -! - call write_snapshot_h5() - if (with_xdmf) then - call write_snapshot_xdmf() - if (master) call write_snapshot_xdmf_master() - end if + case(snapshot_hdf5) + call write_snapshot_h5() + if (with_xdmf) then + call write_snapshot_xdmf() + if (master) call write_snapshot_xdmf_master() + end if #endif /* HDF5 */ - -! stop accounting time for I/O -! - call stop_timer(iio) + case default + call write_snapshot_xml(status) + end select ! increase the snapshot counter and calculate the next snapshot time ! isnap = isnap + 1 tsnap = (ishift + isnap) * hsnap + call stop_timer(iio) #ifdef PROFILE -! stop accounting time for the data writing -! call stop_timer(iow) #endif /* PROFILE */ @@ -847,6 +888,1980 @@ module io !! !=============================================================================== ! +!=============================================================================== +! +! subroutine READ_SNAPSHOT_PARAMETER_STRING: +! ----------------------------------------- +! +! Subroutine reads a string parameter from the restart snapshot. +! +! Arguments: +! +! pname - the parameter name; +! pvalue - the parameter value; +! status - the status flag (the success is 0, failure otherwise); +! +!=============================================================================== +! + subroutine read_snapshot_parameter_string(pname, pvalue, status) + + use iso_fortran_env, only : error_unit + + implicit none + +! subroutine arguments +! + character(len=*), intent(in) :: pname + character(len=*), intent(inout) :: pvalue + integer , intent(inout) :: status + +! local variables +! + logical :: test + character(len=255) :: dname, fname, line + integer(kind=4) :: lun = 103 + integer :: l, u + +! local parameters +! + character(len=*), parameter :: loc = 'IO::read_snapshot_parameter_string' +! +!------------------------------------------------------------------------------- +! + status = 0 + + select case(restart_format) +#ifdef HDF5 + case(snapshot_hdf5) + call read_snapshot_parameter_h5(pname, pvalue, status) +#endif /* HDF5 */ + case default + +! check if the snapshot directory and metafile exist +! + write(dname, "(a,'restart-',i5.5)") trim(respath), nrest + + inquire(file = dname, exist = test) + if (.not. test) then + write(error_unit,"('[',a,']: ',a)") trim(loc), & + trim(dname) // " does not exists!" + status = 121 + return + end if + + write(fname,"(a,'/metadata.xml')") trim(dname) + inquire(file = fname, exist = test) + if (.not. test) then + write(error_unit,"('[',a,']: ',a)") trim(loc), & + trim(fname) // " does not exists!" + status = 121 + return + end if + +! read requested parameter from the file +! + open(newunit = lun, file = fname, status = 'old') +10 read(lun, fmt = "(a)", end = 20) line + l = index(line, trim(adjustl(pname))) + if (l > 0) then + l = index(line, '>') + 1 + u = index(line, '<', back = .true.) - 1 + pvalue = trim(adjustl(line(l:u))) + end if + go to 10 +20 close(lun) + + end select + +!------------------------------------------------------------------------------- +! + end subroutine read_snapshot_parameter_string +! +!=============================================================================== +! +! subroutine READ_SNAPSHOT_PARAMETER_INTEGER: +! ------------------------------------------ +! +! Subroutine reads an integer parameter from the restart snapshot. +! +! Arguments: +! +! pname - the parameter name; +! pvalue - the parameter value; +! status - the status flag (the success is 0, failure otherwise); +! +!=============================================================================== +! + subroutine read_snapshot_parameter_integer(pname, pvalue, status) + + use iso_fortran_env, only : error_unit + + implicit none + +! subroutine arguments +! + character(len=*), intent(in) :: pname + integer , intent(inout) :: pvalue + integer , intent(inout) :: status + +! local variables +! + logical :: test + character(len=255) :: dname, fname, line, svalue + integer(kind=4) :: lun = 103 + integer :: l, u + +! local parameters +! + character(len=*), parameter :: loc = 'IO::read_snapshot_parameter_integer' +! +!------------------------------------------------------------------------------- +! + status = 0 + + select case(restart_format) +#ifdef HDF5 + case(snapshot_hdf5) + call read_snapshot_parameter_h5(pname, pvalue, status) +#endif /* HDF5 */ + case default + +! check if the snapshot directory and metafile exist +! + write(dname, "(a,'restart-',i5.5)") trim(respath), nrest + + inquire(file = dname, exist = test) + if (.not. test) then + write(error_unit,"('[',a,']: ',a)") trim(loc), & + trim(dname) // " does not exists!" + status = 121 + return + end if + + write(fname,"(a,'/metadata.xml')") trim(dname) + inquire(file = fname, exist = test) + if (.not. test) then + write(error_unit,"('[',a,']: ',a)") trim(loc), & + trim(fname) // " does not exists!" + status = 121 + return + end if + +! read parameter from the file +! + open(newunit = lun, file = fname, status = 'old') +10 read(lun, fmt = "(a)", end = 20) line + l = index(line, trim(adjustl(pname))) + if (l > 0) then + l = index(line, '>') + 1 + u = index(line, '<', back = .true.) - 1 + svalue = trim(adjustl(line(l:u))) + read(svalue, fmt = *) pvalue + end if + go to 10 +20 close(lun) + + end select + +!------------------------------------------------------------------------------- +! + end subroutine read_snapshot_parameter_integer +! +!=============================================================================== +! +! subroutine READ_SNAPSHOT_PARAMETER_DOUBLE: +! ----------------------------------------- +! +! Subroutine reads a floating point parameter from the restart snapshot. +! +! Arguments: +! +! pname - the parameter name; +! pvalue - the parameter value; +! status - the status flag (the success is 0, failure otherwise); +! +!=============================================================================== +! + subroutine read_snapshot_parameter_double(pname, pvalue, status) + + use iso_fortran_env, only : error_unit + + implicit none + +! subroutine arguments +! + character(len=*), intent(in) :: pname + real(kind=8) , intent(inout) :: pvalue + integer , intent(inout) :: status + +! local variables +! + logical :: test + character(len=255) :: dname, fname, line, svalue + integer(kind=4) :: lun = 103 + integer :: l, u + +! local parameters +! + character(len=*), parameter :: loc = 'IO::read_snapshot_parameter_double' +! +!------------------------------------------------------------------------------- +! + status = 0 + + select case(restart_format) +#ifdef HDF5 + case(snapshot_hdf5) + call read_snapshot_parameter_h5(pname, pvalue, status) +#endif /* HDF5 */ + case default + +! check if the snapshot directory and metafile exist +! + write(dname, "(a,'restart-',i5.5)") trim(respath), nrest + + inquire(file = dname, exist = test) + if (.not. test) then + write(error_unit,"('[',a,']: ',a)") trim(loc), & + trim(dname) // " does not exists!" + status = 121 + return + end if + + write(fname,"(a,'/metadata.xml')") trim(dname) + inquire(file = fname, exist = test) + if (.not. test) then + write(error_unit,"('[',a,']: ',a)") trim(loc), & + trim(fname) // " does not exists!" + status = 121 + return + end if + +! read parameter from the file +! + open(newunit = lun, file = fname, status = 'old') +10 read(lun, fmt = "(a)", end = 20) line + l = index(line, trim(adjustl(pname))) + if (l > 0) then + l = index(line, '>') + 1 + u = index(line, '<', back = .true.) - 1 + svalue = trim(adjustl(line(l:u))) + read(svalue, fmt = *) pvalue + end if + go to 10 +20 close(lun) + + end select + +!------------------------------------------------------------------------------- +! + end subroutine read_snapshot_parameter_double +! +!=============================================================================== +! +! subroutine READ_RESTART_SNAPSHOT_XML: +! ------------------------------------ +! +! Subroutine reads restart snapshot, i.e. parameters, meta and data blocks +! stored in the XML+binary format restart files and reconstructs +! the data structure in order to resume a terminated job. +! +! Arguments: +! +! status - the return flag to inform if subroutine succeeded or failed; +! +!=============================================================================== +! + subroutine read_restart_snapshot_xml(status) + + use blocks , only : block_meta, block_data, pointer_meta, list_meta + use blocks , only : ndims, ns => nsides, nc => nchildren + use blocks , only : append_metablock, append_datablock, link_blocks + use blocks , only : get_mblocks + use blocks , only : set_last_id, get_last_id + use blocks , only : metablock_set_id, metablock_set_process + use blocks , only : metablock_set_refinement + use blocks , only : metablock_set_configuration + use blocks , only : metablock_set_level, metablock_set_position + use blocks , only : metablock_set_coordinates, metablock_set_bounds + use blocks , only : metablock_set_leaf + use blocks , only : change_blocks_process + use coordinates , only : nn => bcells, ncells, nghosts, minlev, maxlev + use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax + use evolution , only : step, time, dt, dtn + use forcing , only : nmodes, fcoefs, einj + use hash , only : xxh64_integer, xxh64_long, xxh64_double, & + xxh64_complex + use iso_fortran_env, only : error_unit +#ifdef MPI + use mesh , only : redistribute_blocks +#endif /* MPI */ + use mpitools , only : nprocs, nproc + use random , only : gentype, nseeds, set_seeds + + implicit none + + integer, intent(out) :: status + +! local variables +! + logical :: test + character(len=255) :: dname, fname, line, sname, svalue + integer :: id, il, iu, nl, nx, nm, nd, nv, i, j, k, l, n, p, nu + integer(kind=4) :: lndims, lnprocs, lnproc, lmblocks, lnleafs, llast_id + integer(kind=4) :: ldblocks, lncells, lnseeds, lnmodes + real(kind=8) :: deinj + +! local pointers +! + type(block_meta), pointer :: pmeta + type(block_data), pointer :: pdata + +! local variables +! + integer(kind=4) :: lun = 104 + integer(kind=8) :: hash + +! hash strings +! + character(len=22) :: hfield, hchild, hface, hedge, hcorner, hbound + character(len=22) :: hids, harray, hseed, hforce + +! local arrays +! + integer(kind=4), dimension(:) , allocatable :: ids + integer(kind=4), dimension(:,:) , allocatable :: fields + integer(kind=4), dimension(:,:) , allocatable :: children +#if NDIMS == 2 + integer(kind=4), dimension(:,:,:,:) , allocatable :: edges + integer(kind=4), dimension(:,:,:) , allocatable :: corners +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + integer(kind=4), dimension(:,:,:,:,:) , allocatable :: faces + integer(kind=4), dimension(:,:,:,:,:) , allocatable :: edges + integer(kind=4), dimension(:,:,:,:) , allocatable :: corners +#endif /* NDIMS == 3 */ + integer(kind=8), dimension(:,:) , allocatable :: seeds + real(kind=8) , dimension(:,:,:) , allocatable :: bounds + real(kind=8) , dimension(:,:,:,:,:,:), allocatable :: arrays + +! array of pointer used during job restart +! + type(pointer_meta), dimension(:), allocatable :: barray + +! local parameters +! + character(len=*), parameter :: loc = 'IO::read_restart_snapshot_xml()' + character(len=*), parameter :: fmt = "(a,a,'_',i6.6,'.',a)" +! +!------------------------------------------------------------------------------- +! + status = 0 + +! check if the snapshot directory exists +! + write(dname, "(a,'restart-',i5.5)") trim(respath), nrest + + inquire(file = dname, exist = test) + if (.not. test) then + write(*,*) trim(dname) // " does not exists!" + status = 121 + return + end if + dname = trim(dname) // "/" + + write(fname,"(a,'metadata.xml')") trim(dname) + inquire(file = fname, exist = test) + if (.not. test) then + write(*,*) trim(fname) // " does not exists!" + status = 121 + return + end if + +! read attributes from the metadata file +! + open(newunit = lun, file = fname, status = 'old') +10 read(lun, fmt = "(a)", end = 20) line + il = index(line, 'name="') + if (il > 0) then + iu = index(line, '"', back = .true.) + write(sname,*) line(il+6:iu-1) + il = index(line, '>') + 1 + iu = index(line, '<', back = .true.) - 1 + write(svalue,*) line(il:iu) + select case(trim(adjustl(sname))) + case('ndims') + read(svalue, fmt = *) lndims + case('nprocs') + read(svalue, fmt = *) lnprocs + case('nproc') + read(svalue, fmt = *) lnproc + case('mblocks') + read(svalue, fmt = *) lmblocks + case('dblocks') + read(svalue, fmt = *) ldblocks + case('nleafs') + read(svalue, fmt = *) lnleafs + case('last_id') + read(svalue, fmt = *) llast_id + case('ncells') + read(svalue, fmt = *) lncells + case('nseeds') + read(svalue, fmt = *) lnseeds + case('step') + read(svalue, fmt = *) step + case('isnap') + read(svalue, fmt = *) isnap + case('nvars') + read(svalue, fmt = *) nv + case('nmodes') + read(svalue, fmt = *) lnmodes + case('xmin') + read(svalue, fmt = *) xmin + case('xmax') + read(svalue, fmt = *) xmax + case('ymin') + read(svalue, fmt = *) ymin + case('ymax') + read(svalue, fmt = *) ymax + case('zmin') + read(svalue, fmt = *) zmin + case('zmax') + read(svalue, fmt = *) zmax + case('time') + read(svalue, fmt = *) time + case('dt') + read(svalue, fmt = *) dt + case('dtn') + read(svalue, fmt = *) dtn + case('fields') + read(svalue, fmt = *) hfield + case('children') + read(svalue, fmt = *) hchild + case('faces') + read(svalue, fmt = *) hface + case('edges') + read(svalue, fmt = *) hedge + case('corners') + read(svalue, fmt = *) hcorner + case('bounds') + read(svalue, fmt = *) hbound + case('forcing') + read(svalue, fmt = *) hforce + end select + end if + go to 10 +20 close(lun) + +! check the number of dimensions +! + if (lndims /= NDIMS) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "The number of dimensions does not match!" + return + end if + +! check the block dimensions +! + if (lncells /= ncells) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "The block dimensions do not match!" + return + end if + +! allocate all metablocks +! + do l = 1, lmblocks + call append_metablock(pmeta, status) + end do + +! check if the number of created metablocks is equal to lbmcloks +! + if (lmblocks /= get_mblocks()) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Number of metablocks does not match!" + end if + +! set the last_id +! + call set_last_id(llast_id) + +! get numbers of meta and data blocks +! + nx = llast_id + nm = lmblocks + +! prepare and store metablocks +! + allocate(barray(nx), fields(nm,14), children(nm,nc), bounds(nm,NDIMS,2), & +#if NDIMS == 3 + faces(nm,NDIMS,ns,ns,ns), & + edges(nm,NDIMS,ns,ns,ns), corners(nm,ns,ns,ns), & +#else /* NDIMS == 3 */ + edges(nm,NDIMS,ns,ns), corners(nm,ns,ns), & +#endif /* NDIMS == 3 */ + stat = status) + + if (status == 0) then + + fields(:,:) = -1 + children(:,:) = -1 +#if NDIMS == 3 + faces(:,:,:,:,:) = -1 + edges(:,:,:,:,:) = -1 + corners(:,:,:,:) = -1 +#else /* NDIMS == 3 */ + edges(:,:,:,:) = -1 + corners(:,:,:) = -1 +#endif /* NDIMS == 3 */ + bounds(:,:,:) = 0.0d+00 + +! read metablocks from binary files and check hashes +! + write(fname,"(a,'metablock_fields.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(fields)) + read(lun, rec = 1) fields + close(lun) + read(hfield(7:), fmt = "(1z16)") hash + if (hash /= xxh64_integer(size(fields), fields)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + write(fname,"(a,'metablock_children.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(children)) + read(lun, rec = 1) children + close(lun) + read(hchild(7:), fmt = "(1z16)") hash + if (hash /= xxh64_integer(size(children), children)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if +#if NDIMS == 3 + write(fname,"(a,'metablock_faces.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(faces)) + read(lun, rec = 1) faces + close(lun) + read(hface(7:), fmt = "(1z16)") hash + if (hash /= xxh64_integer(size(faces), faces)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if +#endif /* NDIMS == 3 */ + write(fname,"(a,'metablock_edges.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(edges)) + read(lun, rec = 1) edges + close(lun) + read(hedge(7:), fmt = "(1z16)") hash + if (hash /= xxh64_integer(size(edges), edges)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + write(fname,"(a,'metablock_corners.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(corners)) + read(lun, rec = 1) corners + close(lun) + read(hcorner(7:), fmt = "(1z16)") hash + if (hash /= xxh64_integer(size(corners), corners)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + write(fname,"(a,'metablock_bounds.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(bounds)) + read(lun, rec = 1) bounds + close(lun) + read(hbound(7:), fmt = "(1z16)") hash + if (hash /= xxh64_double(size(bounds), bounds)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + +! iterate over all meta blocks and restore their fields +! + l = 0 + pmeta => list_meta + do while(associated(pmeta)) + + l = l + 1 + + barray(fields(l,1))%ptr => pmeta + + call metablock_set_id (pmeta, fields(l, 1)) + call metablock_set_process (pmeta, fields(l, 2)) + call metablock_set_level (pmeta, fields(l, 3)) + call metablock_set_configuration(pmeta, fields(l, 4)) + call metablock_set_refinement (pmeta, fields(l, 5)) + call metablock_set_position (pmeta, fields(l, 6), fields(l, 7), & + fields(l, 8)) + call metablock_set_coordinates (pmeta, fields(l, 9), fields(l,10), & + fields(l,11)) + call metablock_set_bounds (pmeta, bounds(l,1,1), bounds(l,1,2), & + bounds(l,2,1), bounds(l,2,2), & + bounds(l,3,1), bounds(l,3,2)) + + if (fields(l,12) == 1) call metablock_set_leaf(pmeta) + + pmeta => pmeta%next + + end do ! over all meta blocks + +! iterate over all meta blocks and restore their pointers +! + l = 0 + pmeta => list_meta + do while(associated(pmeta)) + + l = l + 1 + + if (fields(l,14) > 0) pmeta%parent => barray(fields(l,14))%ptr + + do p = 1, nc + if (children(l,p) > 0) then + pmeta%child(p)%ptr => barray(children(l,p))%ptr + end if + end do ! p = 1, nc + +#if NDIMS == 2 + do j = 1, ns + do i = 1, ns + do n = 1, NDIMS + if (edges(l,n,i,j) > 0) & + pmeta%edges(i,j,n)%ptr => barray(edges(l,n,i,j))%ptr + end do ! NDIMS + if (corners(l,i,j) > 0) & + pmeta%corners(i,j)%ptr => barray(corners(l,i,j))%ptr + end do ! i = 1, ns + end do ! j = 1, ns +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + do k = 1, ns + do j = 1, ns + do i = 1, ns + do n = 1, NDIMS + if (faces(l,n,i,j,k) > 0) & + pmeta%faces(i,j,k,n)%ptr => barray(faces(l,n,i,j,k))%ptr + if (edges(l,n,i,j,k) > 0) & + pmeta%edges(i,j,k,n)%ptr => barray(edges(l,n,i,j,k))%ptr + end do ! NDIMS + if (corners(l,i,j,k) > 0) & + pmeta%corners(i,j,k)%ptr => barray(corners(l,i,j,k))%ptr + end do ! i = 1, ns + end do ! j = 1, ns + end do ! k = 1, ns +#endif /* NDIMS == 3 */ + + pmeta => pmeta%next + + end do ! over all meta blocks + + if (allocated(fields)) deallocate(fields) + if (allocated(children)) deallocate(children) + if (allocated(bounds)) deallocate(bounds) +#if NDIMS == 3 + if (allocated(faces)) deallocate(faces) +#endif /* NDIMS == 3 */ + if (allocated(edges)) deallocate(edges) + if (allocated(corners)) deallocate(corners) + + end if ! allocation + +! check the number of forcing modes +! + if (lnmodes == nmodes .and. lnmodes > 0) then + write(fname,"(a,'forcing_coefficients.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(fcoefs)) + read(lun, rec = 1) fcoefs + close(lun) + read(hforce(7:), fmt = "(1z16)") hash + if (hash /= xxh64_complex(size(fcoefs), fcoefs)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + else + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "The number of forcing modes does not match!" + end if + +! if the number of processes is bigger after the restart than before +! + if (nprocs >= lnprocs) then + + if (nproc < lnprocs) then + + write(fname,fmt) trim(dname), "datablocks", nproc, "xml" + inquire(file = fname, exist = test) + if (.not. test) then + write(*,*) trim(fname) // " does not exists!" + status = 121 + return + end if + +! read attributes from the metadata file +! + open(newunit = lun, file = fname, status = 'old') +30 read(lun, fmt = "(a)", end = 40) line + il = index(line, 'name="') + if (il > 0) then + iu = index(line, '"', back = .true.) + write(sname,*) line(il+6:iu-1) + il = index(line, '>') + 1 + iu = index(line, '<', back = .true.) - 1 + write(svalue,*) line(il:iu) + select case(trim(adjustl(sname))) + case('dblocks') + read(svalue, fmt = *) nd + case('einj') + read(svalue, fmt = *) einj + case('ids') + read(svalue, fmt = *) hids + case('arrays') + read(svalue, fmt = *) harray + case('seeds') + read(svalue, fmt = *) hseed + end select + end if + go to 30 +40 close(lun) + +#if NDIMS == 3 + allocate(ids(nd), arrays(nd,3,nv,nn,nn,nn), stat = status) +#else /* NDIMS == 3 */ + allocate(ids(nd), arrays(nd,3,nv,nn,nn, 1), stat = status) +#endif /* NDIMS == 3 */ + + if (status == 0) then + + write(fname, fmt) trim(dname), "datablock_ids", nproc, "bin" + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct' , recl = sizeof(ids)) + read(lun, rec = 1) ids + close(lun) + read(hids(7:), fmt = "(1z16)") hash + if (hash /= xxh64_integer(size(ids), ids)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + write(fname, fmt) trim(dname), "datablock_arrays", nproc, "bin" + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct' , recl = sizeof(arrays)) + read(lun, rec = 1) arrays + close(lun) + read(harray(7:), fmt = "(1z16)") hash + if (hash /= xxh64_double(size(arrays), arrays)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + + do l = 1, nd + call append_datablock(pdata, status) + call link_blocks(barray(ids(l))%ptr, pdata) + + pdata%q (:,:,:,:) = arrays(l,1,:,:,:,:) + pdata%u0(:,:,:,:) = arrays(l,2,:,:,:,:) + pdata%u1(:,:,:,:) = arrays(l,3,:,:,:,:) + end do + + if (allocated(ids)) deallocate(ids) + if (allocated(arrays)) deallocate(arrays) + + end if ! allocation + +! restore PRNG seeds +! + allocate(seeds(4,lnseeds), stat = status) + + if (status == 0) then + + write(fname, fmt) trim(dname), "random_seeds", nproc, "bin" + open(unit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(seeds)) + read(lun, rec = 1) seeds + close(lun) + read(hseed(7:), fmt = "(1z16)") hash + if (hash /= xxh64_long(size(seeds), seeds)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + call set_seeds(lnseeds, seeds(:,:), .false.) + + if (allocated(seeds)) deallocate(seeds) + + end if ! allocation + + else ! nproc < lnprocs + + write(fname,fmt) trim(dname), "datablocks", 0, "xml" + inquire(file = fname, exist = test) + if (.not. test) then + write(*,*) trim(fname) // " does not exists!" + status = 121 + return + end if + +! read attributes from the metadata file +! + open(newunit = lun, file = fname, status = 'old') +50 read(lun, fmt = "(a)", end = 60) line + il = index(line, 'name="') + if (il > 0) then + iu = index(line, '"', back = .true.) + write(sname,*) line(il+6:iu-1) + il = index(line, '>') + 1 + iu = index(line, '<', back = .true.) - 1 + write(svalue,*) line(il:iu) + select case(trim(adjustl(sname))) + case('seeds') + read(svalue, fmt = *) hseed + end select + end if + go to 50 +60 close(lun) + +! restore PRNG seeds for remaining processes +! + if (trim(gentype) == "same") then + + allocate(seeds(4,lnseeds), stat = status) + + if (status == 0) then + + write(fname, fmt) trim(dname), "random_seeds", 0, "bin" + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(seeds)) + read(lun, rec = 1) seeds + close(lun) + read(hseed(7:), fmt = "(1z16)") hash + if (hash /= xxh64_long(size(seeds), seeds)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + call set_seeds(lnseeds, seeds(:,:), .false.) + + if (allocated(seeds)) deallocate(seeds) + + end if ! allocation + + end if ! gentype == "same" + + end if ! nproc < nprocs + + else ! nprocs < lnprocs + +! divide files between processes +! + i = mod(lnprocs, nprocs) + j = lnprocs / nprocs + do p = 0, nprocs + k = 0 + do n = 0, p + nl = k + if (n < i) then + nu = k + j + else + nu = k + j - 1 + end if + k = nu + 1 + end do + do n = nl, nu + call change_blocks_process(n, p) + end do + end do + k = 0 + do n = 0, nproc + nl = k + if (n < i) then + nu = k + j + else + nu = k + j - 1 + end if + k = nu + 1 + end do + + do n = nl, nu + + write(fname,fmt) trim(dname), "datablocks", n, "xml" + inquire(file = fname, exist = test) + if (.not. test) then + write(*,*) trim(fname) // " does not exists!" + status = 121 + return + end if + +! read attributes from the metadata file +! + open(newunit = lun, file = fname, status = 'old') +70 read(lun, fmt = "(a)", end = 80) line + il = index(line, 'name="') + if (il > 0) then + iu = index(line, '"', back = .true.) + write(sname,*) line(il+6:iu-1) + il = index(line, '>') + 1 + iu = index(line, '<', back = .true.) - 1 + write(svalue,*) line(il:iu) + select case(trim(adjustl(sname))) + case('dblocks') + read(svalue, fmt = *) nd + case('einj') + read(svalue, fmt = *) deinj + case('ids') + read(svalue, fmt = *) hids + case('arrays') + read(svalue, fmt = *) harray + case('seeds') + read(svalue, fmt = *) hseed + end select + end if + go to 70 +80 close(lun) + einj = einj + deinj + +#if NDIMS == 3 + allocate(ids(nd), arrays(nd,3,nv,nn,nn,nn), stat = status) +#else /* NDIMS == 3 */ + allocate(ids(nd), arrays(nd,3,nv,nn,nn, 1), stat = status) +#endif /* NDIMS == 3 */ + + if (status == 0) then + + write(fname, fmt) trim(dname), "datablock_ids", n, "bin" + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct' , recl = sizeof(ids)) + read(lun, rec = 1) ids + close(lun) + read(hids(7:), fmt = "(1z16)") hash + if (hash /= xxh64_integer(size(ids), ids)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + write(fname, fmt) trim(dname), "datablock_arrays", n, "bin" + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct' , recl = sizeof(arrays)) + read(lun, rec = 1) arrays + close(lun) + read(harray(7:), fmt = "(1z16)") hash + if (hash /= xxh64_double(size(arrays), arrays)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + + do l = 1, nd + call append_datablock(pdata, status) + call link_blocks(barray(ids(l))%ptr, pdata) + + pdata%q (:,:,:,:) = arrays(l,1,:,:,:,:) + pdata%u0(:,:,:,:) = arrays(l,2,:,:,:,:) + pdata%u1(:,:,:,:) = arrays(l,3,:,:,:,:) + end do + + if (allocated(ids)) deallocate(ids) + if (allocated(arrays)) deallocate(arrays) + + end if ! allocation + + end do ! n = nl, nu + +! restore seeds +! + allocate(seeds(4,lnseeds), stat = status) + + if (status == 0) then + + write(fname, fmt) trim(dname), "random_seeds", nproc, "bin" + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(seeds)) + read(lun, rec = 1) seeds + close(lun) + read(hseed(7:), fmt = "(1z16)") hash + if (hash /= xxh64_long(size(seeds), seeds)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "'" // trim(fname) // "' seems to be corrupted!" + end if + call set_seeds(lnseeds, seeds(:,:), .false.) + + if (allocated(seeds)) deallocate(seeds) + + end if ! allocation + + end if ! nprocs >= lnprocs + + if (allocated(barray)) deallocate(barray) + +#ifdef MPI +! redistribute blocks between processors +! + call redistribute_blocks() +#endif /* MPI */ + +!------------------------------------------------------------------------------- +! + end subroutine read_restart_snapshot_xml +! +!=============================================================================== +! +! subroutine WRITE_RESTART_SNAPSHOT_XML: +! ------------------------------------- +! +! Subroutine saves a restart snapshot, i.e. parameters, meta and data blocks +! using the XML format for parameters and binary format for meta and data +! block fields. +! +! Arguments: +! +! nrun - the snapshot number; +! status - the status flag to inform if subroutine succeeded or failed; +! +!=============================================================================== +! + subroutine write_restart_snapshot_xml(nrun, status) + + use blocks , only : block_meta, block_data, list_meta, list_data + use blocks , only : get_mblocks, get_dblocks, get_nleafs + use blocks , only : get_last_id + use blocks , only : ns => nsides, nc => nchildren + use coordinates , only : nn => bcells, ncells, nghosts, minlev, maxlev + use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax + use coordinates , only : bdims => domain_base_dims + use equations , only : eqsys, eos, nv + use evolution , only : step, time, dt, dtn + use forcing , only : nmodes, fcoefs, einj + use hash , only : xxh64_integer, xxh64_long, xxh64_double, & + xxh64_complex + use iso_fortran_env, only : error_unit + use mpitools , only : nprocs, nproc + use parameters , only : get_parameter_file + use problems , only : problem_name + use random , only : gentype, nseeds, get_seeds + + implicit none + +! input and output arguments +! + integer, intent(in) :: nrun + integer, intent(out) :: status + +! local variables +! + logical :: test + character(len=64) :: dname, fname + integer(kind=8) :: hash + integer(kind=4) :: lun = 103 + integer :: nd, nl, nm, nx, i, j, k, l, n, p + +! hash strings +! + character(len=22) :: hfield, hchild, hface, hedge, hcorner, hbound + character(len=22) :: hids, harray, hseed, hforce + +! local pointers +! + type(block_meta), pointer :: pmeta + type(block_data), pointer :: pdata + +! local arrays +! + integer(kind=4), dimension(:) , allocatable :: ids + integer(kind=4), dimension(:,:) , allocatable :: fields + integer(kind=4), dimension(:,:) , allocatable :: children +#if NDIMS == 2 + integer(kind=4), dimension(:,:,:,:) , allocatable :: edges + integer(kind=4), dimension(:,:,:) , allocatable :: corners +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + integer(kind=4), dimension(:,:,:,:,:) , allocatable :: faces + integer(kind=4), dimension(:,:,:,:,:) , allocatable :: edges + integer(kind=4), dimension(:,:,:,:) , allocatable :: corners +#endif /* NDIMS == 3 */ + integer(kind=8), dimension(:,:) , allocatable :: seeds + real(kind=8) , dimension(:,:,:) , allocatable :: bounds + real(kind=8) , dimension(:,:,:,:,:,:), allocatable :: arrays + +! local parameters +! + character(len=*), parameter :: loc = "IO::write_restart_snapshot_xml()" + character(len=*), parameter :: fmt = "(a,a,'_',i6.6,'.',a)" +! +!------------------------------------------------------------------------------- +! + status = 0 + +! create snapshot directory, if it does not exists +! + write(dname, "('restart-',i5.5)") nrun + + inquire(file = dname, exist = test) + do while(.not. test) + if (.not. test) call system("mkdir -p " // trim(dname)) + inquire(file = dname, exist = test) + end do + dname = trim(dname) // "/" + +! get numbers of meta and data blocks +! + nx = get_last_id() + nm = get_mblocks() + nd = get_dblocks() + nl = get_nleafs() + +! only process 0 stores the metadata +! + if (nproc == 0) then + +! copy the parameter file to the restart directory +! + call get_parameter_file(fname, status) + if (status == 0) then + call system("cp -a " // trim(fname) // " " // trim(dname)) + else + write(error_unit,"('[',a,']: ',a)") trim(loc), & + "Cannot get the location of parameter file!" + return + end if + +! prepare and store metablocks +! + allocate(fields(nm,14), children(nm,nc), bounds(nm,3,2), & +#if NDIMS == 3 + faces(nm,NDIMS,ns,ns,ns), & + edges(nm,NDIMS,ns,ns,ns), corners(nm,ns,ns,ns), & +#else /* NDIMS == 3 */ + edges(nm,NDIMS,ns,ns), corners(nm,ns,ns), & +#endif /* NDIMS == 3 */ + stat = status) + + if (status == 0) then + + fields(:,:) = -1 + children(:,:) = -1 +#if NDIMS == 3 + faces(:,:,:,:,:) = -1 + edges(:,:,:,:,:) = -1 + corners(:,:,:,:) = -1 +#else /* NDIMS == 3 */ + edges(:,:,:,:) = -1 + corners(:,:,:) = -1 +#endif /* NDIMS == 3 */ + bounds(:,:,:) = 0.0d+00 + + l = 0 + pmeta => list_meta + do while(associated(pmeta)) + + l = l + 1 + + fields(l, 1) = pmeta%id + fields(l, 2) = pmeta%process + fields(l, 3) = pmeta%level + fields(l, 4) = pmeta%conf + fields(l, 5) = pmeta%refine + fields(l, 6) = pmeta%pos(1) + fields(l, 7) = pmeta%pos(2) +#if NDIMS == 3 + fields(l, 8) = pmeta%pos(3) +#endif /* NDIMS == 3 */ + fields(l, 9) = pmeta%coords(1) + fields(l,10) = pmeta%coords(2) +#if NDIMS == 3 + fields(l,11) = pmeta%coords(3) +#endif /* NDIMS == 3 */ + + if (pmeta%leaf) fields(l,12) = 1 + if (associated(pmeta%data) ) fields(l,13) = 1 + if (associated(pmeta%parent)) fields(l,14) = pmeta%parent%id + + do p = 1, nc + if (associated(pmeta%child(p)%ptr)) & + children(l,p) = pmeta%child(p)%ptr%id + end do + +#if NDIMS == 2 + do j = 1, ns + do i = 1, ns + do n = 1, NDIMS + if (associated(pmeta%edges(i,j,n)%ptr)) & + edges(l,n,i,j) = pmeta%edges(i,j,n)%ptr%id + end do ! NDIMS + if (associated(pmeta%corners(i,j)%ptr)) & + corners(l,i,j) = pmeta%corners(i,j)%ptr%id + end do ! i = 1, ns + end do ! j = 1, ns +#endif /* NDIMS == 2 */ +#if NDIMS == 3 + do k = 1, ns + do j = 1, ns + do i = 1, ns + do n = 1, NDIMS + if (associated(pmeta%faces(i,j,k,n)%ptr)) & + faces(l,n,i,j,k) = pmeta%faces(i,j,k,n)%ptr%id + if (associated(pmeta%edges(i,j,k,n)%ptr)) & + edges(l,n,i,j,k) = pmeta%edges(i,j,k,n)%ptr%id + end do ! NDIMS + if (associated(pmeta%corners(i,j,k)%ptr)) & + corners(l,i,j,k) = pmeta%corners(i,j,k)%ptr%id + end do ! i = 1, ns + end do ! j = 1, ns + end do ! k = 1, ns +#endif /* NDIMS == 3 */ + + bounds(l,1,1) = pmeta%xmin + bounds(l,1,2) = pmeta%xmax + bounds(l,2,1) = pmeta%ymin + bounds(l,2,2) = pmeta%ymax +#if NDIMS == 3 + bounds(l,3,1) = pmeta%zmin + bounds(l,3,2) = pmeta%zmax +#endif /* NDIMS == 3 */ + + pmeta => pmeta%next + + end do ! metablocks + +! get hashes for metadablock data +! + hash = xxh64_integer(size(fields), fields) + write(hfield ,"('xxh64:', 1z0.16)") hash + hash = xxh64_integer(size(children), children) + write(hchild ,"('xxh64:', 1z0.16)") hash +#if NDIMS == 3 + hash = xxh64_integer(size(faces), faces) + write(hface ,"('xxh64:', 1z0.16)") hash +#endif /* NDIMS == 3 */ + hash = xxh64_integer(size(edges), edges) + write(hedge ,"('xxh64:', 1z0.16)") hash + hash = xxh64_integer(size(corners), corners) + write(hcorner,"('xxh64:', 1z0.16)") hash + hash = xxh64_double(size(bounds), bounds) + write(hbound ,"('xxh64:', 1z0.16)") hash + if (nmodes > 0) then + hash = xxh64_complex(size(fcoefs), fcoefs) + write(hforce,"('xxh64:', 1z0.16)") hash + end if + +! store metablock data +! + write(fname,"(a,'metablock_fields.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(fields), status = 'replace') + write(lun, rec = 1) fields + close(lun) + write(fname,"(a,'metablock_children.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(children), status = 'replace') + write(lun, rec = 1) children + close(lun) + write(fname,"(a,'metablock_bounds.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(bounds), status = 'replace') + write(lun, rec = 1) bounds + close(lun) +#if NDIMS == 3 + write(fname,"(a,'metablock_faces.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(faces), status = 'replace') + write(lun, rec = 1) faces + close(lun) +#endif /* NDIMS == 3 */ + write(fname,"(a,'metablock_edges.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(edges), status = 'replace') + write(lun, rec = 1) edges + close(lun) + write(fname,"(a,'metablock_corners.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(corners), status = 'replace') + write(lun, rec = 1) corners + close(lun) + if (nmodes > 0) then + write(fname,"(a,'forcing_coefficients.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(fcoefs), status = 'replace') + write(lun, rec = 1) fcoefs + close(lun) + end if + + if (allocated(fields)) deallocate(fields) + if (allocated(children)) deallocate(children) + if (allocated(bounds)) deallocate(bounds) +#if NDIMS == 3 + if (allocated(faces)) deallocate(faces) +#endif /* NDIMS == 3 */ + if (allocated(edges)) deallocate(edges) + if (allocated(corners)) deallocate(corners) + + else + write(error_unit,"('[',a,']: ',a)") trim(loc), & + "Cannot allocate space for metablocks!" + status = 1001 + return + end if ! allocation + +! store metadata (parameters and attributes) +! + write(fname,"(a,'metadata.xml')") trim(dname) + open(newunit = lun, file = fname, status = 'replace') + write(lun,"(a)") '' + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "problem" , problem_name) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "nprocs" , nprocs) + call write_attribute_xml(lun, "nproc" , nproc) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "eqsys" , eqsys) + call write_attribute_xml(lun, "eos" , eos) + call write_attribute_xml(lun, "nvars" , nv) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "ndims" , NDIMS) + call write_attribute_xml(lun, "xblocks" , bdims(1)) + call write_attribute_xml(lun, "yblocks" , bdims(2)) +#if NDIMS == 3 + call write_attribute_xml(lun, "zblocks" , bdims(3)) +#endif /* NDIMS */ + call write_attribute_xml(lun, "xmin" , xmin) + call write_attribute_xml(lun, "xmax" , xmax) + call write_attribute_xml(lun, "ymin" , ymin) + call write_attribute_xml(lun, "ymax" , ymax) +#if NDIMS == 3 + call write_attribute_xml(lun, "zmin" , zmin) + call write_attribute_xml(lun, "zmax" , zmax) +#endif /* NDIMS */ + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "minlev" , minlev) + call write_attribute_xml(lun, "maxlev" , maxlev) + call write_attribute_xml(lun, "ncells" , ncells) + call write_attribute_xml(lun, "nghosts" , nghosts) + call write_attribute_xml(lun, "bcells" , nn) + call write_attribute_xml(lun, "nchildren", nc) + call write_attribute_xml(lun, "mblocks" , nm) + call write_attribute_xml(lun, "nleafs" , nl) + call write_attribute_xml(lun, "last_id" , nx) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "step" , step) + call write_attribute_xml(lun, "time" , time) + call write_attribute_xml(lun, "dt" , dt) + call write_attribute_xml(lun, "dtn" , dtn) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "nmodes" , nmodes) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "gentype" , gentype) + call write_attribute_xml(lun, "nseeds" , nseeds) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "isnap" , isnap) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "fields" , hfield) + call write_attribute_xml(lun, "children" , hchild) +#if NDIMS == 3 + call write_attribute_xml(lun, "faces" , hface ) +#endif /* NDIMS == 3 */ + call write_attribute_xml(lun, "edges" , hedge ) + call write_attribute_xml(lun, "corners" , hcorner) + call write_attribute_xml(lun, "bounds" , hbound) + if (nmodes > 0) then + call write_attribute_xml(lun, "forcing", hforce) + end if + write(lun,"(a)") '' + write(lun,"(a)") '' + close(lun) + + end if ! meta data file is stored only on the master process + + if (nd > 0) then + +#if NDIMS == 3 + allocate(ids(nd), arrays(nd,3,nv,nn,nn,nn), stat = status) +#else /* NDIMS == 3 */ + allocate(ids(nd), arrays(nd,3,nv,nn,nn, 1), stat = status) +#endif /* NDIMS == 3 */ + + if (status == 0) then + + arrays = 0.0d+00 + + l = 0 + pdata => list_data + do while(associated(pdata)) + + l = l + 1 + ids(l) = pdata%meta%id + arrays(l,1,:,:,:,:) = pdata%q (:,:,:,:) + arrays(l,2,:,:,:,:) = pdata%u0(:,:,:,:) + arrays(l,3,:,:,:,:) = pdata%u1(:,:,:,:) + + pdata => pdata%next + + end do ! data blocks + +! get hashes for block IDs and arrays +! + hash = xxh64_integer(size(ids), ids) + write(hids ,"('xxh64:', 1z0.16)") hash + hash = xxh64_double(size(arrays), arrays) + write(harray,"('xxh64:', 1z0.16)") hash + +! store block IDs and arrays +! + write(fname,fmt) trim(dname), "datablock_ids", nproc, "bin" + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(ids), status = 'replace') + write(lun, rec = 1) ids + close(lun) + + write(fname,fmt) trim(dname), "datablock_arrays", nproc, "bin" + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(arrays), status = 'replace') + write(lun, rec = 1) arrays + close(lun) + + if (allocated(ids)) deallocate(ids) + if (allocated(arrays)) deallocate(arrays) + + else + write(error_unit,"('[',a,']: ',a)") trim(loc), & + "Cannot allocate space for datablocks!" + status = 1001 + return + end if ! allocation + + end if + +! store PRNG seeds +! + allocate(seeds(4,nseeds), stat = status) + + if (status == 0) then + + call get_seeds(seeds(:,:)) + +! get the hash for seeds +! + hash = xxh64_long(size(seeds), seeds) + write(hseed,"('xxh64:', 1z0.16)") hash + +! store seeds +! + write(fname,fmt) trim(dname), "random_seeds", nproc, "bin" + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(seeds), status = 'replace') + write(lun, rec = 1) seeds + close(lun) + + if (allocated(seeds)) deallocate(seeds) + + else + + write(error_unit,"('[',a,']: ',a)") trim(loc), & + "Cannot allocate space for random generator seeds!" + status = 1001 + return + + end if + +! prepare and store data block info +! + write(fname,fmt) trim(dname), "datablocks", nproc, "xml" + open(newunit = lun, file = fname, status = 'replace') + write(lun,"(a)") '' + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "nprocs" , nprocs) + call write_attribute_xml(lun, "nproc" , nproc) + call write_attribute_xml(lun, "ndims" , NDIMS) + call write_attribute_xml(lun, "ncells" , ncells) + call write_attribute_xml(lun, "nghosts", nghosts) + call write_attribute_xml(lun, "bcells" , nn) + call write_attribute_xml(lun, "dblocks", nd) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "einj" , einj) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "gentype", gentype) + call write_attribute_xml(lun, "nseeds" , nseeds) + write(lun,"(a)") '' + write(lun,"(a)") '' + if (nd > 0) then + call write_attribute_xml(lun, "ids" , hids) + call write_attribute_xml(lun, "arrays", harray) + end if + if (nseeds > 0) then + call write_attribute_xml(lun, "seeds", hseed) + end if + write(lun,"(a)") '' + write(lun,"(a)") '' + close(lun) + +!------------------------------------------------------------------------------- +! + end subroutine write_restart_snapshot_xml +! +!=============================================================================== +! +! subroutine WRITE_SNAPSHOT_XML: +! ----------------------------- +! +! Subroutine saves a regular snapshot, i.e. parameters, leafs and data blocks +! using the XML format for metadata and binary format for meta and data +! block fields. +! +! Arguments: +! +! status - the status flag to inform if subroutine succeeded or failed; +! +!=============================================================================== +! + subroutine write_snapshot_xml(status) + + use blocks , only : block_meta, block_data, list_meta, list_data + use blocks , only : get_dblocks, get_nleafs + use blocks , only : ns => nsides, nc => nchildren + use coordinates , only : nn => bcells, ncells, nghosts, minlev, maxlev + use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax + use coordinates , only : bdims => domain_base_dims + use equations , only : eqsys, eos, nv, pvars, gamma, csnd + use evolution , only : step, time, dt, dtn + use hash , only : xxh64_integer, xxh64_double + use iso_fortran_env, only : error_unit + use mpitools , only : nprocs, nproc + use parameters , only : get_parameter_file + use problems , only : problem_name + use sources , only : viscosity, resistivity + + implicit none + +! input and output arguments +! + integer, intent(out) :: status + +! local variables +! + logical :: test + character(len=64) :: dname, fname + character(len=256) :: vars + integer(kind=4) :: lun = 103, dun = 104 + integer(kind=8) :: hash + integer :: nd, nl, l, p + +! hash strings +! + character(len=22) :: hfield, hbound, hvar + +! local pointers +! + type(block_meta), pointer :: pmeta + type(block_data), pointer :: pdata + +! local arrays +! + integer(kind=4), dimension(:) , allocatable :: ids + integer(kind=4), dimension(:,:) , allocatable :: fields + real(kind=8) , dimension(:,:,:) , allocatable :: bounds + real(kind=8) , dimension(:,:,:,:,:), allocatable :: arrays + +! local parameters +! + character(len=*), parameter :: loc = "IO::write_snapshot_xml()" + character(len=*), parameter :: fmt = "(a,a,'_',i6.6,'.',a)" +! +!------------------------------------------------------------------------------- +! + status = 0 + +! create snapshot directory, if it does not exists +! + write(dname, "('snapshot-',i9.9)") isnap + + inquire(file = dname, exist = test) + do while(.not. test) + if (.not. test) call system("mkdir -p " // trim(dname)) + inquire(file = dname, exist = test) + end do + dname = trim(dname) // "/" + +! get numbers of meta and data blocks, leafs, and prepare stored variables +! + nd = get_dblocks() + nl = get_nleafs() + vars = "" + do l = 1, nv + vars = trim(vars) // " " // trim(pvars(l)) + end do + +! only process 0 stores the metadata +! + if (nproc == 0) then + +! copy the parameter file to the restart directory +! + call get_parameter_file(fname, status) + if (status == 0) then + call system("cp -a " // trim(fname) // " " // trim(dname)) + else + write(error_unit,"('[',a,']: ',a)") trim(loc), & + "Cannot get the location of parameter file!" + return + end if + +! prepare and store metablocks +! + allocate(fields(nl,5), bounds(nl,3,2), stat = status) + + if (status == 0) then + + fields(:,:) = -1 + bounds(:,:,:) = 0.0d+00 + + l = 0 + pmeta => list_meta + do while(associated(pmeta)) + + if (pmeta%leaf) then + + l = l + 1 + + fields(l,1) = pmeta%id + fields(l,2) = pmeta%level + fields(l,3) = pmeta%coords(1) + fields(l,4) = pmeta%coords(2) +#if NDIMS == 3 + fields(l,5) = pmeta%coords(3) +#endif /* NDIMS == 3 */ + + bounds(l,1,1) = pmeta%xmin + bounds(l,1,2) = pmeta%xmax + bounds(l,2,1) = pmeta%ymin + bounds(l,2,2) = pmeta%ymax +#if NDIMS == 3 + bounds(l,3,1) = pmeta%zmin + bounds(l,3,2) = pmeta%zmax +#endif /* NDIMS == 3 */ + + end if ! leaf + + pmeta => pmeta%next + + end do ! metablocks + +! get hashes for metadablock data +! + hash = xxh64_integer(size(fields), fields) + write(hfield ,"('xxh64:', 1z0.16)") hash + hash = xxh64_double(size(bounds), bounds) + write(hbound ,"('xxh64:', 1z0.16)") hash + +! store metablock data +! + write(fname,"(a,'metablock_fields.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(fields), status = 'replace') + write(lun, rec = 1) fields + close(lun) + write(fname,"(a,'metablock_bounds.bin')") trim(dname) + open(newunit = lun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(bounds), status = 'replace') + write(lun, rec = 1) bounds + close(lun) + + if (allocated(fields)) deallocate(fields) + if (allocated(bounds)) deallocate(bounds) + + else + write(error_unit,"('[',a,']: ',a)") trim(loc), & + "Cannot allocate space for metablocks!" + status = 1001 + return + end if ! allocation + +! store metadata (parameters and attributes) +! + write(fname,"(a,'metadata.xml')") trim(dname) + open(newunit = lun, file = fname, status = 'replace') + write(lun,"(a)") '' + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "problem" , problem_name) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "nprocs" , nprocs) + call write_attribute_xml(lun, "nproc" , nproc) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "eqsys" , eqsys) + call write_attribute_xml(lun, "eos" , eos) + call write_attribute_xml(lun, "nvars" , nv) + call write_attribute_xml(lun, "gamma" , gamma) + call write_attribute_xml(lun, "csnd" , csnd) + call write_attribute_xml(lun, "viscosity" , viscosity) + call write_attribute_xml(lun, "resistivity", resistivity) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "ndims" , NDIMS) + call write_attribute_xml(lun, "xblocks" , bdims(1)) + call write_attribute_xml(lun, "yblocks" , bdims(2)) +#if NDIMS == 3 + call write_attribute_xml(lun, "zblocks" , bdims(3)) +#endif /* NDIMS */ + call write_attribute_xml(lun, "xmin" , xmin) + call write_attribute_xml(lun, "xmax" , xmax) + call write_attribute_xml(lun, "ymin" , ymin) + call write_attribute_xml(lun, "ymax" , ymax) +#if NDIMS == 3 + call write_attribute_xml(lun, "zmin" , zmin) + call write_attribute_xml(lun, "zmax" , zmax) +#endif /* NDIMS */ + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "minlev" , minlev) + call write_attribute_xml(lun, "maxlev" , maxlev) + call write_attribute_xml(lun, "ncells" , ncells) + call write_attribute_xml(lun, "nghosts" , nghosts) + call write_attribute_xml(lun, "bcells" , nn) + call write_attribute_xml(lun, "nleafs" , nl) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "step" , step) + call write_attribute_xml(lun, "time" , time) + call write_attribute_xml(lun, "dt" , dt) + call write_attribute_xml(lun, "dtn" , dtn) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "isnap" , isnap) + call write_attribute_xml(lun, "variables", trim(vars)) + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "fields" , hfield) + call write_attribute_xml(lun, "bounds" , hbound) + write(lun,"(a)") '' + write(lun,"(a)") '' + close(lun) + + end if ! meta data file is stored only on the master process + +! prepare and store data block info +! + write(fname,fmt) trim(dname), "datablocks", nproc, "xml" + open(newunit = lun, file = fname, status = 'replace') + write(lun,"(a)") '' + write(lun,"(a)") '' + write(lun,"(a)") '' + call write_attribute_xml(lun, "nprocs" , nprocs) + call write_attribute_xml(lun, "nproc" , nproc) + call write_attribute_xml(lun, "ndims" , NDIMS) + call write_attribute_xml(lun, "ncells" , ncells) + call write_attribute_xml(lun, "nghosts" , nghosts) + call write_attribute_xml(lun, "bcells" , nn) + call write_attribute_xml(lun, "nvars" , nv) + call write_attribute_xml(lun, "dblocks" , nd) + call write_attribute_xml(lun, "variables", trim(vars)) + write(lun,"(a)") '' + write(lun,"(a)") '' + + if (nd > 0) then + +#if NDIMS == 3 + allocate(ids(nd), arrays(nv,nd,nn,nn,nn), stat = status) +#else /* NDIMS == 3 */ + allocate(ids(nd), arrays(nv,nd,nn,nn, 1), stat = status) +#endif /* NDIMS == 3 */ + + if (status == 0) then + + l = 0 + pdata => list_data + do while(associated(pdata)) + + l = l + 1 + ids(l) = pdata%meta%id + arrays(:,l,:,:,:) = pdata%q(:,:,:,:) + + pdata => pdata%next + + end do ! data blocks + + write(fname,fmt) trim(dname), "datablock_ids", nproc, "bin" + open(newunit = dun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(ids), status = 'replace') + write(dun, rec = 1) ids + close(dun) + hash = xxh64_integer(size(ids), ids) + write(hvar,"('xxh64:', 1z0.16)") hash + call write_attribute_xml(lun, "ids", hvar) + + do p = 1, nv + write(fname,fmt) trim(dname), "datablock_" // trim(pvars(p)), & + nproc, "bin" + open(newunit = dun, file = fname, form = 'unformatted', & + access = 'direct', recl = sizeof(arrays(p,:,:,:,:)), & + status = 'replace') + write(dun, rec = 1) arrays(p,:,:,:,:) + close(dun) + hash = xxh64_double(size(arrays(p,:,:,:,:)), arrays(p,:,:,:,:)) + write(hvar,"('xxh64:', 1z0.16)") hash + call write_attribute_xml(lun, trim(pvars(p)), hvar) + end do + + if (allocated(ids)) deallocate(ids) + if (allocated(arrays)) deallocate(arrays) + + else + write(error_unit,"('[',a,']: ',a)") trim(loc), & + "Cannot allocate space for datablocks!" + status = 1001 + return + end if ! allocation + + end if + + write(lun,"(a)") '' + write(lun,"(a)") '' + close(lun) + +!------------------------------------------------------------------------------- +! + end subroutine write_snapshot_xml +! +!=============================================================================== +! +! subroutine WRITE_ATTRIBUTE_XML_STRING: +! ------------------------------------- +! +! Subroutine writes a string attribute in XML format to specified +! file handler. +! +! Arguments: +! +! lun - the file handler to write to; +! aname - the name of attribute; +! avalue - the value of attribute; +! +!=============================================================================== +! + subroutine write_attribute_xml_string(lun, aname, avalue) + + implicit none + +! input and output arguments +! + integer , intent(in) :: lun + character(len=*), intent(in) :: aname + character(len=*), intent(in) :: avalue + +! local parameters +! + character(len=*), parameter :: afmt = "('',a,'')" +! +!------------------------------------------------------------------------------- +! + write(lun,afmt) "string", trim(adjustl(aname)), trim(adjustl(avalue)) + +!------------------------------------------------------------------------------- +! + end subroutine write_attribute_xml_string +! +!=============================================================================== +! +! subroutine WRITE_ATTRIBUTE_XML_INTEGER: +! -------------------------------------- +! +! Subroutine writes an integer attribute in XML format to specified +! file handler. +! +! Arguments: +! +! lun - the file handler to write to; +! aname - the name of attribute; +! avalue - the value of attribute; +! +!=============================================================================== +! + subroutine write_attribute_xml_integer(lun, aname, avalue) + + implicit none + +! input and output arguments +! + integer , intent(in) :: lun + character(len=*), intent(in) :: aname + integer(kind=4) , intent(in) :: avalue + +! local variables +! + character(len=32) :: svalue + +! local parameters +! + character(len=*), parameter :: afmt = "('',a,'')" +! +!------------------------------------------------------------------------------- +! + write(svalue,"(1i32)") avalue + write(lun,afmt) "integer", trim(adjustl(aname)), trim(adjustl(svalue)) + +!------------------------------------------------------------------------------- +! + end subroutine write_attribute_xml_integer +! +!=============================================================================== +! +! subroutine WRITE_ATTRIBUTE_XML_DOUBLE: +! -------------------------------------- +! +! Subroutine writes a double precision attribute in XML format to specified +! file handler. +! +! Arguments: +! +! lun - the file handler to write to; +! aname - the name of attribute; +! avalue - the value of attribute; +! +!=============================================================================== +! + subroutine write_attribute_xml_double(lun, aname, avalue) + + implicit none + +! input and output arguments +! + integer , intent(in) :: lun + character(len=*), intent(in) :: aname + real(kind=8) , intent(in) :: avalue + +! local variables +! + character(len=32) :: svalue + +! local parameters +! + character(len=*), parameter :: afmt = "('',a,'')" +! +!------------------------------------------------------------------------------- +! + write(svalue,"(1es32.20)") avalue + write(lun,afmt) "double", trim(adjustl(aname)), trim(adjustl(svalue)) + +!------------------------------------------------------------------------------- +! + end subroutine write_attribute_xml_double #ifdef HDF5 ! !=============================================================================== @@ -1275,9 +3290,11 @@ module io ! import external procedures and variables ! use blocks , only : change_blocks_process + use forcing , only : einj use hdf5 , only : hid_t use hdf5 , only : H5F_ACC_RDONLY_F use hdf5 , only : h5fis_hdf5_f, h5fopen_f, h5fclose_f + use hdf5 , only : h5gopen_f, h5gclose_f use iso_fortran_env, only : error_unit #ifdef MPI use mesh , only : redistribute_blocks @@ -1295,9 +3312,10 @@ module io ! local variables ! character(len=255) :: fl, msg - integer(hid_t) :: fid + integer(hid_t) :: fid, gid integer :: err, lfile logical :: info + real(kind=8) :: deinj ! local parameters ! @@ -1315,13 +3333,8 @@ module io ! not exist decrease it until the file corresponding to lower process number ! is found; ! - info = .false. - lfile = nproc + 1 - do while (.not. info .and. lfile > 0) - lfile = lfile - 1 - write (fl, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, lfile - inquire(file = fl, exist = info) - end do + write (fl, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, 0 + inquire(file = fl, exist = info) ! quit, if file does not exist ! @@ -1435,6 +3448,21 @@ module io else ! file is opened +! restore injected energy +! + call h5gopen_f(fid, 'attributes', gid, err) + if (err >= 0) then + call read_attribute(gid, 'einj', einj) + call h5gclose_f(gid, err) + if (err /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot close the group!" + end if + else + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot open the group!" + end if + ! read data blocks ! call read_datablocks_h5(fid) @@ -1517,6 +3545,22 @@ module io else ! file is opened +! restore injected energy +! + call h5gopen_f(fid, 'attributes', gid, err) + if (err >= 0) then + call read_attribute(gid, 'einj', deinj) + einj = einj + deinj + call h5gclose_f(gid, err) + if (err /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot close the group!" + end if + else + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot open the group!" + end if + ! read data blocks ! call read_datablocks_h5(fid) @@ -1797,6 +3841,7 @@ module io use coordinates , only : periodic use equations , only : eqsys, eos, gamma, csnd use evolution , only : step, time, dt, dtn + use forcing , only : nmodes, einj, fcoefs use hdf5 , only : hid_t use hdf5 , only : h5gcreate_f, h5gclose_f use iso_fortran_env, only : error_unit @@ -1824,7 +3869,7 @@ module io ! local allocatable arrays ! - integer(kind=4), dimension(:), allocatable :: seeds + integer(kind=8), dimension(:,:), allocatable :: seeds ! local parameters ! @@ -1905,8 +3950,19 @@ module io ! store the vector attributes ! dims(1:NDIMS) = ncells - call write_attribute(gid, 'dims' , dims) - call write_attribute(gid, 'bdims', bdims) + call write_attribute(gid, 'dims' , dims) + call write_attribute(gid, 'bdims' , bdims) + call write_attribute(gid, 'xblocks', bdims(1)) + call write_attribute(gid, 'yblocks', bdims(2)) + call write_attribute(gid, 'zblocks', bdims(3)) + +! store forcing parameters +! + call write_attribute(gid, 'nmodes', nmodes) + call write_attribute(gid, 'einj' , einj) + if (nmodes > 0) then + call write_attribute(gid, 'fcoefs', fcoefs) + end if ! store random number generator seed values ! @@ -1914,15 +3970,15 @@ module io ! allocate space for seeds ! - allocate(seeds(nseeds)) + allocate(seeds(4,nseeds)) ! get the seed values ! - call get_seeds(seeds) + call get_seeds(seeds(:,:)) ! store them in the current group ! - call write_attribute(gid, 'seeds', seeds(:)) + call write_attribute(gid, 'seeds', seeds(:,:)) ! deallocate seed array ! @@ -1974,11 +4030,12 @@ module io use coordinates , only : ncells, nghosts use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax use evolution , only : step, time, dt, dtn + use forcing , only : nmodes, fcoefs use hdf5 , only : hid_t, hsize_t use hdf5 , only : h5gopen_f, h5gclose_f use iso_fortran_env, only : error_unit use mpitools , only : nprocs, nproc - use random , only : nseeds, set_seeds + use random , only : nseeds, set_seeds, gentype ! local variables are not implicit by default ! @@ -1993,7 +4050,7 @@ module io integer(hid_t) :: gid integer :: ierr, l integer :: lndims, lmblocks, lnleafs, llast_id - integer :: lncells, lnproc, lnseeds + integer :: lncells, lnproc, lnseeds, lnmodes integer :: status ! local pointers @@ -2002,7 +4059,7 @@ module io ! allocatable arrays ! - integer(kind=4), dimension(:), allocatable :: seeds + integer(kind=8), dimension(:,:), allocatable :: seeds ! local parameters ! @@ -2062,21 +4119,24 @@ module io , "The block dimensions do not match!" end if -! allocate space for seeds +! restore forcing coefficients ! - allocate(seeds(lnseeds)) + call read_attribute(gid, 'nmodes', lnmodes) + if (lnmodes == nmodes .and. lnmodes > 0) then + call read_attribute(gid, 'fcoefs', fcoefs) + else + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "The number of driving modes does not match!" + end if -! store them in the current group +! restore seeds ! - call read_attribute(gid, 'seeds', seeds(:)) - -! set the seed values -! - call set_seeds(lnseeds, seeds(:), nproc /= lnproc) - -! deallocate seed array -! - deallocate(seeds) + if (trim(gentype) == "same") then + allocate(seeds(4,0:lnseeds-1)) + call read_attribute(gid, 'seeds', seeds(:,:)) + call set_seeds(lnseeds, seeds(:,:), .false.) + deallocate(seeds) + end if ! allocate all metablocks ! @@ -4201,6 +6261,219 @@ module io !------------------------------------------------------------------------------- ! end subroutine write_vector_attribute_double_h5 +! +!=============================================================================== +! +! subroutine WRITE_ARRAY_ATTRIBUTE_LONG_H5: +! ---------------------------------------- +! +! Subroutine stores a 2D array of the integer attribute in the group provided +! by an identifier and the attribute name. +! +! Arguments: +! +! gid - the group identifier to which the attribute should be linked; +! aname - the attribute name; +! avalue - the attribute values; +! +!=============================================================================== +! + subroutine write_array_attribute_long_h5(gid, aname, avalue) + +! import procedures and variables from other modules +! + use hdf5 , only : hid_t, hsize_t, H5T_STD_I64LE + use hdf5 , only : h5screate_simple_f, h5sclose_f + use hdf5 , only : h5acreate_f, h5awrite_f, h5aclose_f + use iso_fortran_env, only : error_unit + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*) , intent(in) :: aname + integer(kind=8) , dimension(:,:), intent(in) :: avalue + +! local variables +! + integer(hid_t) :: sid, aid + integer(hsize_t), dimension(2) :: am = (/ 1, 1 /) + integer :: ierr + +! local parameters +! + character(len=*), parameter :: loc = 'IO::write_array_attribute_long_h5()' +! +!------------------------------------------------------------------------------- +! +! set the proper attribute length +! + am(1) = size(avalue, 1) + am(2) = size(avalue, 2) + +! create space for the attribute value +! + call h5screate_simple_f(2, am, sid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot create space for attribute :" // trim(aname) + return + end if + +! create the attribute in the given group +! + call h5acreate_f(gid, aname, H5T_STD_I64LE, sid, aid, ierr) + if (ierr == 0) then + +! write the attribute data +! + call h5awrite_f(aid, H5T_STD_I64LE, avalue, am, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot write the attribute data in :" // trim(aname) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot close attribute :" // trim(aname) + end if + + else + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot create attribute :" // trim(aname) + end if + +! release the space +! + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot close space for attribute :" // trim(aname) + end if + +!------------------------------------------------------------------------------- +! + end subroutine write_array_attribute_long_h5 +! +!=============================================================================== +! +! subroutine WRITE_ARRAY_ATTRIBUTE_COMPLEX_H5: +! ------------------------------------------- +! +! Subroutine stores a 2D array of the double precision complex attribute +! in the group provided by an identifier and the attribute name. +! +! Arguments: +! +! gid - the group identifier to which the attribute should be linked; +! aname - the attribute name; +! avalue - the attribute values; +! +!=============================================================================== +! + subroutine write_array_attribute_complex_h5(gid, aname, avalue) + +! import procedures and variables from other modules +! + use hdf5 , only : hid_t, hsize_t, H5T_NATIVE_DOUBLE + use hdf5 , only : h5screate_simple_f, h5sclose_f + use hdf5 , only : h5acreate_f, h5awrite_f, h5aclose_f + use iso_fortran_env, only : error_unit + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*) , intent(in) :: aname + complex(kind=8) , dimension(:,:), intent(in) :: avalue + +! local variables +! + integer(hid_t) :: sid, aid + integer(hsize_t), dimension(3) :: am = (/ 1, 1, 1 /) + integer :: ierr + +! allocatable arrays +! + real(kind=8), dimension(:,:,:), allocatable :: tvalue + +! local parameters +! + character(len=*), parameter :: loc = 'IO::write_array_attribute_complex_h5()' +! +!------------------------------------------------------------------------------- +! +! set the proper attribute length +! + am(1) = 2 + am(2) = size(avalue, 1) + am(3) = size(avalue, 2) + +! create space for the attribute value +! + call h5screate_simple_f(3, am, sid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot create space for attribute :" // trim(aname) + return + end if + +! allocate space for temporary array and copy the input array values +! + allocate(tvalue(am(1),am(2),am(3))) + tvalue(1,:,:) = dreal(avalue(:,:)) + tvalue(2,:,:) = dimag(avalue(:,:)) + +! create the attribute in the given group +! + call h5acreate_f(gid, aname, H5T_NATIVE_DOUBLE, sid, aid, ierr) + if (ierr == 0) then + +! write the attribute data +! + call h5awrite_f(aid, H5T_NATIVE_DOUBLE, tvalue, am, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot write the attribute data in :" // trim(aname) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot close attribute :" // trim(aname) + end if + + else + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot create attribute :" // trim(aname) + end if + +! deallocate temporary array +! + deallocate(tvalue) + +! release the space +! + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot close space for attribute :" // trim(aname) + end if + +!------------------------------------------------------------------------------- +! + end subroutine write_array_attribute_complex_h5 !=============================================================================== ! @@ -4640,6 +6913,268 @@ module io ! !=============================================================================== ! +! subroutine READ_ARRAY_ATTRIBUTE_LONG_H5: +! --------------------------------------- +! +! Subroutine reads a 2D array of the integer attribute provided by the group +! identifier to which it is linked and its name. +! +! Arguments: +! +! gid - the group identifier to which the attribute is linked; +! aname - the attribute name; +! avalue - the attribute value; +! +!=============================================================================== +! + subroutine read_array_attribute_long_h5(gid, aname, avalue) + +! import procedures and variables from other modules +! + use hdf5 , only : H5T_STD_I64LE + use hdf5 , only : hid_t, hsize_t + use hdf5 , only : h5aexists_by_name_f, h5aget_space_f + use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f + use hdf5 , only : h5sclose_f, h5sget_simple_extent_dims_f + use iso_fortran_env, only : error_unit + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*) , intent(in) :: aname + integer(kind=8) , dimension(:,:), intent(inout) :: avalue + +! local variables +! + logical :: exists = .false. + integer(hid_t) :: aid, sid + integer(hsize_t), dimension(2) :: am, bm + integer(hsize_t) :: alen + integer :: ierr + +! local parameters +! + character(len=*), parameter :: loc = 'IO::read_array_attribute_long_h5()' +! +!------------------------------------------------------------------------------- +! +! check if the attribute exists in the group provided by gid +! + call h5aexists_by_name_f(gid, '.', aname, exists, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot check if attribute exists :" // trim(aname) + return + end if + if (.not. exists) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Attribute does not exist :" // trim(aname) + return + end if + +! open the attribute +! + call h5aopen_by_name_f(gid, '.', aname, aid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot open attribute :" // trim(aname) + return + end if + +! get the attribute space +! + call h5aget_space_f(aid, sid, ierr) + if (ierr == 0) then + call h5sget_simple_extent_dims_f(sid, am, bm, ierr) + if (ierr /= 2) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot get attribute dimensions :" // trim(aname) + end if + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot close the attribute space :" // trim(aname) + end if + else + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot get the attribute space :" // trim(aname) + return + end if + +! check if the output array is large enough +! + if (am(1) * am(2) > size(avalue)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Attribute too large for output argument :" // trim(aname) + return + end if + +! read attribute value +! + call h5aread_f(aid, H5T_STD_I64LE, avalue, am(:), ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot read attribute :" // trim(aname) + end if + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot close attribute :" // trim(aname) + return + end if + +!------------------------------------------------------------------------------- +! + end subroutine read_array_attribute_long_h5 +! +!=============================================================================== +! +! subroutine READ_ARRAY_ATTRIBUTE_COMPLEX_H5: +! ------------------------------------------ +! +! Subroutine reads a 2D array of the double precision complex attribute +! provided by the group identifier to which it is linked and its name. +! +! Arguments: +! +! gid - the group identifier to which the attribute is linked; +! aname - the attribute name; +! avalue - the attribute value; +! +!=============================================================================== +! + subroutine read_array_attribute_complex_h5(gid, aname, avalue) + +! import procedures and variables from other modules +! + use hdf5 , only : H5T_NATIVE_DOUBLE + use hdf5 , only : hid_t, hsize_t + use hdf5 , only : h5aexists_by_name_f, h5aget_space_f + use hdf5 , only : h5aopen_by_name_f, h5aread_f, h5aclose_f + use hdf5 , only : h5sclose_f, h5sget_simple_extent_dims_f + use iso_fortran_env, only : error_unit + +! local variables are not implicit by default +! + implicit none + +! attribute arguments +! + integer(hid_t) , intent(in) :: gid + character(len=*) , intent(in) :: aname + complex(kind=8) , dimension(:,:), intent(inout) :: avalue + +! local variables +! + logical :: exists = .false. + integer(hid_t) :: aid, sid + integer(hsize_t), dimension(3) :: am, bm + integer(hsize_t) :: alen + integer :: ierr + +! local allocatable arrays +! + real(kind=8), dimension(:,:,:), allocatable :: tvalue + +! local parameters +! + character(len=*), parameter :: loc = 'IO::read_array_attribute_complex_h5()' +! +!------------------------------------------------------------------------------- +! +! check if the attribute exists in the group provided by gid +! + call h5aexists_by_name_f(gid, '.', aname, exists, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot check if attribute exists :" // trim(aname) + return + end if + if (.not. exists) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Attribute does not exist :" // trim(aname) + return + end if + +! open the attribute +! + call h5aopen_by_name_f(gid, '.', aname, aid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot open attribute :" // trim(aname) + return + end if + +! get the attribute space +! + call h5aget_space_f(aid, sid, ierr) + if (ierr == 0) then + call h5sget_simple_extent_dims_f(sid, am, bm, ierr) + if (ierr /= 3) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot get attribute dimensions :" // trim(aname) + end if + call h5sclose_f(sid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot close the attribute space :" // trim(aname) + end if + else + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot get the attribute space :" // trim(aname) + return + end if + +! allocate temporary array for attribute +! + allocate(tvalue(am(1),am(2),am(3))) + +! check if the output array is large enough +! + if (am(1) * am(2) * am(3) > size(tvalue)) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Attribute too large for output argument :" // trim(aname) + return + end if + +! read attribute value +! + call h5aread_f(aid, H5T_NATIVE_DOUBLE, tvalue, am(:), ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot read attribute :" // trim(aname) + end if + +! copy array to complex one +! + avalue(:,:) = dcmplx(tvalue(1,:,:), tvalue(2,:,:)) + +! deallocate temporary array +! + deallocate(tvalue) + +! close the attribute +! + call h5aclose_f(aid, ierr) + if (ierr /= 0) then + write(error_unit,"('[',a,']: ',a)") trim(loc) & + , "Cannot close attribute :" // trim(aname) + return + end if + +!------------------------------------------------------------------------------- +! + end subroutine read_array_attribute_complex_h5 +! +!=============================================================================== +! ! WRITE_ARRAY SUBROUTINES ! !=============================================================================== diff --git a/sources/mpitools.F90 b/sources/mpitools.F90 index caf28e5..fd30174 100644 --- a/sources/mpitools.F90 +++ b/sources/mpitools.F90 @@ -1345,7 +1345,7 @@ module mpitools ! include external procedures and variables ! 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 ! @@ -1353,13 +1353,13 @@ module mpitools ! subroutine arguments ! - integer , intent(in) :: n - complex, dimension(n), intent(inout) :: cbuf - integer , intent(out) :: iret + integer , intent(in) :: n + complex(kind=8), dimension(n), intent(inout) :: cbuf + integer , intent(out) :: iret ! local variables ! - real(kind=8), dimension(n) :: rbuf, ibuf, tbuf + complex(kind=8), dimension(n) :: tbuf ! local parameters ! @@ -1377,10 +1377,7 @@ module mpitools call start_timer(imm) #endif /* PROFILE */ - tbuf(:) = real(cbuf(:)) - 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) + call mpi_allreduce(cbuf, tbuf, n, mpi_double_complex, mpi_sum, comm, iret) #ifdef PROFILE ! stop time accounting for the MPI reduce @@ -1390,7 +1387,7 @@ module mpitools ! substitute the result ! - cbuf(1:n) = cmplx(rbuf(1:n), ibuf(1:n)) + cbuf(:) = tbuf(:) ! check if the operation was successful ! diff --git a/sources/parameters.F90 b/sources/parameters.F90 index 4817866..e03af2c 100644 --- a/sources/parameters.F90 +++ b/sources/parameters.F90 @@ -75,7 +75,7 @@ module parameters ! declare public subroutines ! public :: read_parameters, finalize_parameters - public :: get_parameter + public :: get_parameter_file, get_parameter #ifdef MPI public :: redistribute_parameters #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: ! -------------------------------- ! diff --git a/sources/problems.F90 b/sources/problems.F90 index 8c4a5c9..5c7bae6 100644 --- a/sources/problems.F90 +++ b/sources/problems.F90 @@ -1515,7 +1515,7 @@ module problems use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl use parameters , only : get_parameter - use random , only : randomn + use random , only : randsym ! local variables are not implicit by default ! @@ -1658,10 +1658,10 @@ module problems ! add a random seed velocity component ! do i = 1, nn - q(ivx,i) = q(ivx,i) + vper * randomn() - q(ivy,i) = q(ivy,i) + vper * randomn() + q(ivx,i) = q(ivx,i) + vper * randsym() + q(ivy,i) = q(ivy,i) + vper * randsym() #if NDIMS == 3 - q(ivz,i) = q(ivz,i) + vper * randomn() + q(ivz,i) = q(ivz,i) + vper * randsym() #endif /* NDIMS == 3 */ end do @@ -1729,7 +1729,7 @@ module problems use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp, isl use equations , only : csnd2 use parameters , only : get_parameter - use random , only : randomn + use random , only : randsym ! local variables are not implicit by default ! @@ -1882,7 +1882,7 @@ module problems ! if (vper /= 0.0d+00) then do i = 1, nn - q(ivy,i) = q(ivy,i) + vper * randomn() + q(ivy,i) = q(ivy,i) + vper * randsym() end do end if diff --git a/sources/random.F90 b/sources/random.F90 index cb8be11..c175f6d 100644 --- a/sources/random.F90 +++ b/sources/random.F90 @@ -27,8 +27,10 @@ !! !! References: !! -!! [1] Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating -!! random variables', J. Statist. Software, v5(8) +!! [1] David Blackman and Sebastiano Vigna, +!! "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 ! timer indices ! - integer , save :: iri, irc + integer, save :: iri, irg #endif /* PROFILE */ ! random generator type @@ -56,9 +58,15 @@ module random ! variables to store seeds for random generator ! - integer , save :: kp = 0 - integer , save :: nseeds, lseed - integer(kind=4), dimension(:), allocatable, save :: seeds + integer , save :: kp = 0, np = 0 + integer , save :: nseeds, lseed + 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 ! @@ -67,7 +75,8 @@ module random ! declare public subroutines ! 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 ! @@ -76,77 +85,89 @@ module 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 -! declare all variables as implicit -! implicit none ! subroutine arguments ! - integer, intent(in) :: nthreads, nthread + integer, intent(in) :: nthreads, nthread, nproc + integer, intent(out) :: status ! local variables ! - integer :: i - real(kind=4) :: r + integer :: i + +! local arrays +! + integer(kind=8), dimension(4) :: s ! !------------------------------------------------------------------------------- ! -#ifdef PROFILE -! set timer descriptions -! - call set_timer('random:: initialization' , iri) - call set_timer('random:: number generation', irc) + status = 0 + +#ifdef PROFILE + call set_timer('random:: initialization' , iri) + call set_timer('random:: number generation', irg) -! start accounting time for the random number generator -! call start_timer(iri) #endif /* PROFILE */ -! set the thread number -! - kp = nthread - -! obtain the generator type -! - call get_parameter("gentype", gentype) - -! calculate the number of seeds +! set local copies of the thread and process numbers, +! determine the number of seeds ! + kp = nthread + np = nproc nseeds = nthreads 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 ! - select case(gentype) - case('random') - do i = 0, lseed - call random_number(r) - seeds(i) = 123456789 * r - end do - case default - r = 0.1234567890123456789 - do i = 0, lseed - seeds(i) = 123456789 * r - end do - end select + call get_parameter("gentype", gentype) + select case(gentype) + case('random') + state = 1234567890 + s(1) = splitmix64() + 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 + case default + state = 1234567890 + seeds(1,:) = splitmix64() + seeds(2,:) = splitmix64() + seeds(3,:) = splitmix64() + seeds(4,:) = splitmix64() + end select + + end if #ifdef PROFILE -! stop accounting time for the random number generator -! call stop_timer(iri) #endif /* PROFILE */ @@ -159,31 +180,23 @@ module 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() -! declare all variables as implicit -! implicit none ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for the random number generator -! call start_timer(iri) #endif /* PROFILE */ -! deallocate seeds if they are allocated -! if (allocated(seeds)) deallocate(seeds) #ifdef PROFILE -! stop accounting time for the random number generator -! call stop_timer(iri) #endif /* PROFILE */ @@ -196,59 +209,70 @@ module random ! 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 ! input arguments ! - integer , intent(in) :: np - integer(kind=4), dimension(np), intent(in) :: nseeds - logical , intent(in) :: generate + integer , intent(in) :: n + integer(kind=8), dimension(4,n), intent(in) :: sds + logical , intent(in) :: gen ! local variables ! - integer :: i, l - real(kind=4) :: r + integer :: i + +! local arrays +! + integer(kind=8), dimension(4) :: s ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for the random number generator -! call start_timer(iri) #endif /* PROFILE */ ! set the seeds only if the input array and seeds have the same sizes ! - l = min(lseed, np - 1) - seeds(0:l) = nseeds(1:l+1) - select case(gentype) - case('random') - if (generate) then - do i = 0, lseed - call random_number(r) - seeds(i) = 123456789 * r + if (gen .or. n /= nseeds) then + select case(gentype) + case('random') + state = 1234567890 + s(1) = splitmix64() + s(2) = splitmix64() + s(3) = splitmix64() + s(4) = splitmix64() + do i = 1, nseeds * np + call jump(s(:)) end do - else - do i = l + 1, lseed - call random_number(r) - seeds(i) = 123456789 * r + seeds(:,0) = s(:) + do i = 1, lseed + call jump(s(:)) + seeds(:,i) = s(:) end do - end if - case default - seeds(l+1:lseed) = seeds(0) - end select + case default + state = 1234567890 + seeds(1,:) = splitmix64() + seeds(2,:) = splitmix64() + seeds(3,:) = splitmix64() + seeds(4,:) = splitmix64() + end select + else + seeds(:,:) = sds(:,:) + end if #ifdef PROFILE -! stop accounting time for the random number generator -! call stop_timer(iri) #endif /* PROFILE */ @@ -261,33 +285,33 @@ module random ! 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 ! output arguments ! - integer(kind=4), dimension(nseeds), intent(out) :: rseeds + integer(kind=8), dimension(4,nseeds), intent(out) :: sds ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for the random number generator -! call start_timer(iri) #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 -! stop accounting time for the random number generator -! call stop_timer(iri) #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 -! 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 ! - 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 -! start accounting time for the random number generation -! - call start_timer(irc) + call start_timer(irg) #endif /* PROFILE */ - jsr = seeds(kp) - 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) + randuni = i64tor8 * xoshiro256p() #ifdef PROFILE -! stop accounting time for the random number generation -! - call stop_timer(irc) + call stop_timer(irg) #endif /* PROFILE */ 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 - -! output variables -! - real(kind=8) :: val - -! local variables -! - integer(kind=4) :: jz, jsr ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for the random number generation -! - call start_timer(irc) + call start_timer(irg) #endif /* PROFILE */ - jsr = seeds(kp) - 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) + randsym = 2.0d+00 * (i64tor8 * xoshiro256p() - 5.0d-01) #ifdef PROFILE -! stop accounting time for the random number generation -! - call stop_timer(irc) + call stop_timer(irg) #endif /* PROFILE */ 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 -! output variables -! - real(kind=8) :: val - -! local variables -! - integer(kind=4) :: jz, jsr + logical :: c + real(kind=8) :: u, v, s ! !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for the random number generation -! - call start_timer(irc) + call start_timer(irg) #endif /* PROFILE */ - jsr = seeds(kp) - jz = jsr + c = .true. + 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 ) ) - jsr = ieor( jsr, ishft( jsr, -17 ) ) - jsr = ieor( jsr, ishft( jsr, 5 ) ) - - seeds(kp) = jsr - - val = 0.46566128730d-09 * (jz + jsr) + randnorz = dsqrt(- 2.0d+00 * dlog(s) / s) * dcmplx(u, v) #ifdef PROFILE -! stop accounting time for the random number generation -! - call stop_timer(irc) + call stop_timer(irg) #endif /* PROFILE */ return !------------------------------------------------------------------------------- ! - end function randomn + end function randnorz !=============================================================================== !