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

Merge branch 'master' into gem-reconnection-challenge

Browse Source
This commit is contained in:
Grzegorz Kowal 2023-07-27 18:45:25 -03:00
commit 385ca2eb5d
8 changed files with 453 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
README.md
View File

@ -2,7 +2,7 @@
# **The AMUN Code** # **The AMUN Code**
## Copyright (C) 2008-2023 Grzegorz Kowal ## Copyright (C) 2008-2023 Grzegorz Kowal
[![Build Status](https://ampere.clarokowal.com/api/badges/gkowal/amun-code/status.svg)](https://ampere.clarokowal.com/gkowal/amun-code) [![Build Status](https://drone.amuncode.org/api/badges/gkowal/amun-code/status.svg)](https://drone.amuncode.org/gkowal/amun-code)
AMUN is a parallel code to perform numerical simulations in fluid approximation AMUN is a parallel code to perform numerical simulations in fluid approximation
on uniform or non-uniform (adaptive) meshes. The goal in developing this code is on uniform or non-uniform (adaptive) meshes. The goal in developing this code is

2
python/amunpy/setup.py
View File

@ -5,7 +5,7 @@ with open("README.md", "r", encoding="utf-8") as fh:
setuptools.setup( setuptools.setup(
name="amunpy", name="amunpy",
version="0.9.9", version="0.9.10",
author="Grzegorz Kowal", author="Grzegorz Kowal",
author_email="grzegorz@amuncode.org", author_email="grzegorz@amuncode.org",
description="Python Interface for the AMUN code's snapshots", description="Python Interface for the AMUN code's snapshots",

2
python/amunpy/src/amunpy/__init__.py
View File

@ -21,6 +21,6 @@ __all__ = [ 'AmunXML', 'AmunH5', 'WriteVTK', \
__author__ = "Grzegorz Kowal" __author__ = "Grzegorz Kowal"
__copyright__ = "Copyright 2018-2023 Grzegorz Kowal <grzegorz@amuncode.org>" __copyright__ = "Copyright 2018-2023 Grzegorz Kowal <grzegorz@amuncode.org>"
__version__ = "0.9.9" __version__ = "0.9.10"
__maintainer__ = "Grzegorz Kowal" __maintainer__ = "Grzegorz Kowal"
__email__ = "grzegorz@amuncode.org" __email__ = "grzegorz@amuncode.org"

36
python/amunpy/src/amunpy/amunxml.py
View File

@ -224,6 +224,24 @@ class AmunXML(Amun):
print("File '{}' seems to be corrupted! Proceeding anyway...".format(filename)) print("File '{}' seems to be corrupted! Proceeding anyway...".format(filename))
def __shuffle_decode(self, a, dtype='int64'):
import numpy
s = numpy.dtype(dtype).itemsize
d = [s, len(a) // s]
return numpy.frombuffer(a, dtype="int8").reshape(d).T.tobytes()
def __bytedelta_decode(self, a, dtype='int64'):
import numpy
s = numpy.dtype(dtype).itemsize
d = [s, len(a) // s]
return numpy.cumsum(numpy.frombuffer(a, dtype="int8").reshape(d), axis=-1, dtype='int8').T.tobytes()
def __read_binary_meta(self, dataset, dtype='int32'): def __read_binary_meta(self, dataset, dtype='int32'):
''' '''
Reads binary data of metadata. Reads binary data of metadata.
@ -254,6 +272,15 @@ class AmunXML(Amun):
else: else:
raise Exception("Binary file '{}' compressed in unsupported format {}!".format(fname, comp)) raise Exception("Binary file '{}' compressed in unsupported format {}!".format(fname, comp))
if 'data_filter' in self.binaries[dataset]:
data_filter = self.binaries[dataset]['data_filter']
if data_filter == 'bytedelta':
data = self.__bytedelta_decode(data, dtype=dtype)
elif data_filter == 'shuffle':
data = self.__shuffle_decode(data, dtype=dtype)
else:
raise Exception("Binary file '{}' processed using unsupported filter {}!".format(fname, data_filter))
if 'digest' in self.binaries[dataset]: if 'digest' in self.binaries[dataset]:
htype = self.binaries[dataset]['digest_type'] htype = self.binaries[dataset]['digest_type']
dhash = self.binaries[dataset]['digest'] dhash = self.binaries[dataset]['digest']
@ -302,6 +329,15 @@ class AmunXML(Amun):
else: else:
raise Exception("Binary file '{}' compressed in unsupported format {}!".format(fname, comp)) raise Exception("Binary file '{}' compressed in unsupported format {}!".format(fname, comp))
if 'data_filter' in self.chunks[chunk_number][dataset_name]:
data_filter = self.chunks[chunk_number][dataset_name]['data_filter']
if data_filter == 'bytedelta':
data = self.__bytedelta_decode(data, dtype=dtype)
elif data_filter == 'shuffle':
data = self.__shuffle_decode(data, dtype=dtype)
else:
raise Exception("Binary file '{}' processed using unsupported filter {}!".format(fname, data_filter))
if 'digest' in self.chunks[chunk_number][dataset_name]: if 'digest' in self.chunks[chunk_number][dataset_name]:
htype = self.chunks[chunk_number][dataset_name]['digest_type'] htype = self.chunks[chunk_number][dataset_name]['digest_type']
dhash = self.chunks[chunk_number][dataset_name]['digest'] dhash = self.chunks[chunk_number][dataset_name]['digest']

81
sources/compression.F90
View File

@ -123,8 +123,9 @@ module compression
private private
public :: set_compression, get_compression, compression_bound, compress public :: set_compression, get_compression, compression_bound
public :: compression_suffix public :: compression_none, compression_suffix
public :: encode, compress
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! !
@ -270,7 +271,7 @@ module compression
! subroutine COMPRESS: ! subroutine COMPRESS:
! ------------------- ! -------------------
! !
! Subroutine compressed input buffer using ZSTD compression. ! Subroutine compresses the input buffer.
! !
! Arguments: ! Arguments:
! !
@ -338,6 +339,80 @@ module compression
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine compress end subroutine compress
!
!===============================================================================
!
! subroutine ENCODE:
! -----------------
!
! Subroutine encodes the input buffer for better data compression.
!
! Arguments:
!
! filter_type - the filter type to apply, 1 - shuffle, 2 - bytedelta;
! item_size - the size of each item in bytes;
! bytes - the length of the input data;
! input_ptr - the pointer to the input data;
! output - the output where the encoded data are written;
!
!===============================================================================
!
subroutine encode(filter_type, item_size, bytes, input_ptr, output)
use iso_c_binding, only : c_ptr, c_f_pointer
implicit none
integer(kind=1) , intent(in ) :: filter_type
integer(kind=4) , intent(in ) :: item_size
integer(kind=8) , intent(in ) :: bytes
type(c_ptr) , intent(in ) :: input_ptr
integer(kind=1), dimension(bytes), intent(inout) :: output
integer(kind=1), dimension(:), pointer :: input
integer(kind=8) :: i, j, m, n, item_count
!-------------------------------------------------------------------------------
!
call c_f_pointer(input_ptr, input, [ bytes ])
select case(filter_type)
case(1)
item_count = bytes / item_size
i = 1
j = item_count
do m = 1, item_size
output(i:j) = input(m:bytes:item_size)
i = j + 1
j = j + item_count
end do
case(2)
item_count = bytes / item_size
i = 1
j = item_count
do m = 1, item_size
output(i:j) = input(m:bytes:item_size)
do n = j, i + 1, -1
output(n) = output(n) - output(n-1)
end do
i = j + 1
j = j + item_count
end do
case default
output(:) = input(:)
end select
!-------------------------------------------------------------------------------
!
end subroutine encode
!=============================================================================== !===============================================================================
! !

116
sources/io.F90
View File

@ -118,6 +118,11 @@ module io
character(len=255), save :: cformat = "none" ! compression format character(len=255), save :: cformat = "none" ! compression format
integer , save :: clevel = 3 ! compression level integer , save :: clevel = 3 ! compression level
! data filter applied before the compression:
! (0 - no filter, 1 - shuffle, 2 - bytedelta)
!
integer(kind=1) :: data_filter = 0
! the type of digest to use and its length ! the type of digest to use and its length
! !
integer, save :: hash_type = 0 integer, save :: hash_type = 0
@ -185,7 +190,7 @@ module io
! !
subroutine initialize_io(verbose, status) subroutine initialize_io(verbose, status)
use compression, only : set_compression, get_compression use compression, only : set_compression, get_compression, compression_none
use hash , only : hash_info use hash , only : hash_info
use helpers , only : print_message use helpers , only : print_message
#ifdef HDF5 #ifdef HDF5
@ -199,7 +204,7 @@ module io
integer, intent(out) :: status integer, intent(out) :: status
logical :: test logical :: test
character(len=255) :: dname character(len=255) :: string, dname
character(len=255) :: sformat = "xml" character(len=255) :: sformat = "xml"
character(len=255) :: precise = "off" character(len=255) :: precise = "off"
character(len=255) :: ghosts = "on" character(len=255) :: ghosts = "on"
@ -285,6 +290,18 @@ module io
call get_parameter("compression_level" , clevel) call get_parameter("compression_level" , clevel)
call set_compression(cformat, clevel) call set_compression(cformat, clevel)
string = "none"
call get_parameter("data_filter", string)
select case(string)
case('shuffle', 'SHUFFLE')
data_filter = 1
case('bytedelta', 'BYTEDELTA')
data_filter = 2
case default
data_filter = 0
end select
if (get_compression() == compression_none) data_filter = 0
call get_parameter("digest_type", dtype) call get_parameter("digest_type", dtype)
call hash_info(dtype, hash_type, hash_length) call hash_info(dtype, hash_type, hash_length)
@ -493,6 +510,14 @@ module io
call print_parameter(verbose, "digest type", hash_name(hash_type)) call print_parameter(verbose, "digest type", hash_name(hash_type))
call print_parameter(verbose, "compression format", cformat) call print_parameter(verbose, "compression format", cformat)
call print_parameter(verbose, "compression level", clevel) call print_parameter(verbose, "compression level", clevel)
select case(data_filter)
case(1)
call print_parameter(verbose, "data filter", "shuffle")
case(2)
call print_parameter(verbose, "data filter", "bytedelta")
case default
call print_parameter(verbose, "data filter", "none")
end select
end select end select
call print_parameter(verbose, "precise snapshot intervals", & call print_parameter(verbose, "precise snapshot intervals", &
precise_snapshots, "on") precise_snapshots, "on")
@ -2189,49 +2214,53 @@ module io
write(fname,"(a,'.bin')") "metablock_fields" write(fname,"(a,'.bin')") "metablock_fields"
bytes = size(fields, kind=8) * kind(fields) bytes = size(fields, kind=8) * kind(fields)
call write_binary_xml(dname, fname, c_loc(fields), bytes, dtype, digest) call write_binary_xml(dname, fname, c_loc(fields), bytes, &
kind(fields), dtype, digest)
call write_attribute_xml(lun, "fields", fname, 'int32', & call write_attribute_xml(lun, "fields", fname, 'int32', &
shape(fields), bytes, dtype, digest) shape(fields), bytes, dtype, digest)
write(fname,"(a,'.bin')") "metablock_children" write(fname,"(a,'.bin')") "metablock_children"
bytes = size(children, kind=8) * kind(children) bytes = size(children, kind=8) * kind(children)
call write_binary_xml(dname, fname, c_loc(children), & call write_binary_xml(dname, fname, c_loc(children), bytes, &
bytes, dtype, digest) kind(children), dtype, digest)
call write_attribute_xml(lun, "children", fname, 'int32', & call write_attribute_xml(lun, "children", fname, 'int32', &
shape(children), bytes, dtype, digest) shape(children), bytes, dtype, digest)
#if NDIMS == 3 #if NDIMS == 3
write(fname,"(a,'.bin')") "metablock_faces" write(fname,"(a,'.bin')") "metablock_faces"
bytes = size(faces, kind=8) * kind(faces) bytes = size(faces, kind=8) * kind(faces)
call write_binary_xml(dname, fname, c_loc(faces), bytes, dtype, digest) call write_binary_xml(dname, fname, c_loc(faces), bytes, &
kind(faces), dtype, digest)
call write_attribute_xml(lun, "faces", fname, 'int32', & call write_attribute_xml(lun, "faces", fname, 'int32', &
shape(faces), bytes, dtype, digest) shape(faces), bytes, dtype, digest)
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
write(fname,"(a,'.bin')") "metablock_edges" write(fname,"(a,'.bin')") "metablock_edges"
bytes = size(edges, kind=8) * kind(edges) bytes = size(edges, kind=8) * kind(edges)
call write_binary_xml(dname, fname, c_loc(edges), bytes, dtype, digest) call write_binary_xml(dname, fname, c_loc(edges), bytes, &
kind(edges), dtype, digest)
call write_attribute_xml(lun, "edges", fname, 'int32', & call write_attribute_xml(lun, "edges", fname, 'int32', &
shape(edges), bytes, dtype, digest) shape(edges), bytes, dtype, digest)
write(fname,"(a,'.bin')") "metablock_corners" write(fname,"(a,'.bin')") "metablock_corners"
bytes = size(corners, kind=8) * kind(corners) bytes = size(corners, kind=8) * kind(corners)
call write_binary_xml(dname, fname, c_loc(corners), & call write_binary_xml(dname, fname, c_loc(corners), bytes, &
bytes, dtype, digest) kind(corners), dtype, digest)
call write_attribute_xml(lun, "corners", fname, 'int32', & call write_attribute_xml(lun, "corners", fname, 'int32', &
shape(corners), bytes, dtype, digest) shape(corners), bytes, dtype, digest)
write(fname,"(a,'.bin')") "metablock_bounds" write(fname,"(a,'.bin')") "metablock_bounds"
bytes = size(bounds, kind=8) * kind(bounds) bytes = size(bounds, kind=8) * kind(bounds)
call write_binary_xml(dname, fname, c_loc(bounds), bytes, dtype, digest) call write_binary_xml(dname, fname, c_loc(bounds), bytes, &
kind(bounds), dtype, digest)
call write_attribute_xml(lun, "bounds", fname, 'float64', & call write_attribute_xml(lun, "bounds", fname, 'float64', &
shape(bounds), bytes, dtype, digest) shape(bounds), bytes, dtype, digest)
if (nmodes > 0) then if (nmodes > 0) then
write(fname,"(a,'.bin')") "forcing_coefficients" write(fname,"(a,'.bin')") "forcing_coefficients"
bytes = size(fcoefs, kind=8) * kind(fcoefs) * 2 bytes = size(fcoefs, kind=8) * kind(fcoefs) * 2
call write_binary_xml(dname, fname, c_loc(fcoefs), & call write_binary_xml(dname, fname, c_loc(fcoefs), bytes, &
bytes, dtype, digest) kind(fcoefs), dtype, digest)
call write_attribute_xml(lun, "forcing", fname, 'complex64', & call write_attribute_xml(lun, "forcing", fname, 'complex64', &
shape(fcoefs), bytes, dtype, digest) shape(fcoefs), bytes, dtype, digest)
end if end if
@ -2295,15 +2324,15 @@ module io
write(aname,"('_',i6.6)") l write(aname,"('_',i6.6)") l
write(fname,"('datablock_prim_',i6.6,a,'.bin')") nproc, trim(aname) write(fname,"('datablock_prim_',i6.6,a,'.bin')") nproc, trim(aname)
bytes = size(pdata%q, kind=8) * kind(pdata%q) bytes = size(pdata%q, kind=8) * kind(pdata%q)
call write_binary_xml(dname, fname, c_loc(pdata%q), & call write_binary_xml(dname, fname, c_loc(pdata%q), bytes, &
bytes, dtype, digest) kind(pdata%q), dtype, digest)
call write_attribute_xml(lun, "prim" // trim(aname), fname, & call write_attribute_xml(lun, "prim" // trim(aname), fname, &
'float64', shape(pdata%q), bytes, dtype, digest) 'float64', shape(pdata%q), bytes, dtype, digest)
write(fname,"('datablock_cons_',i6.6,a,'.bin')") nproc, trim(aname) write(fname,"('datablock_cons_',i6.6,a,'.bin')") nproc, trim(aname)
bytes = size(pdata%uu, kind=8) * kind(pdata%uu) bytes = size(pdata%uu, kind=8) * kind(pdata%uu)
call write_binary_xml(dname, fname, c_loc(pdata%uu), & call write_binary_xml(dname, fname, c_loc(pdata%uu), bytes, &
bytes, dtype, digest) kind(pdata%uu), dtype, digest)
call write_attribute_xml(lun, "cons" // trim(aname), fname, & call write_attribute_xml(lun, "cons" // trim(aname), fname, &
'float64', shape(pdata%uu), bytes, dtype, digest) 'float64', shape(pdata%uu), bytes, dtype, digest)
@ -2313,7 +2342,8 @@ module io
write(fname,"(a,'_',a,'_',i6.6,'.bin')") "datablock", "ids", nproc write(fname,"(a,'_',a,'_',i6.6,'.bin')") "datablock", "ids", nproc
bytes = size(ids, kind=8) * kind(ids) bytes = size(ids, kind=8) * kind(ids)
call write_binary_xml(dname, fname, c_loc(ids), bytes, dtype, digest) call write_binary_xml(dname, fname, c_loc(ids), bytes, &
kind(ids), dtype, digest)
call write_attribute_xml(lun, "ids", fname, 'int32', & call write_attribute_xml(lun, "ids", fname, 'int32', &
shape(ids), bytes, dtype, digest) shape(ids), bytes, dtype, digest)
@ -2335,7 +2365,8 @@ module io
write(fname,"(a,'_',a,'_',i6.6,'.bin')") "random", "seeds", nproc write(fname,"(a,'_',a,'_',i6.6,'.bin')") "random", "seeds", nproc
bytes = size(seeds, kind=8) * kind(seeds) bytes = size(seeds, kind=8) * kind(seeds)
call write_binary_xml(dname, fname, c_loc(seeds), bytes, dtype, digest) call write_binary_xml(dname, fname, c_loc(seeds), bytes, &
kind(seeds), dtype, digest)
call write_attribute_xml(lun, "seeds", fname, 'int64', & call write_attribute_xml(lun, "seeds", fname, 'int64', &
shape(seeds), bytes, dtype, digest) shape(seeds), bytes, dtype, digest)
@ -2543,15 +2574,15 @@ module io
write(fname,"(a,'.bin')") "metablock_fields" write(fname,"(a,'.bin')") "metablock_fields"
dbytes = size(fields, kind=8) * kind(fields) dbytes = size(fields, kind=8) * kind(fields)
call write_binary_xml(dname, fname, c_loc(fields), & call write_binary_xml(dname, fname, c_loc(fields), dbytes, &
dbytes, hash_type, ddigest, cbytes, cdigest) kind(fields), hash_type, ddigest, cbytes, cdigest)
call write_attribute_xml(lun, "fields", fname, 'int32', & call write_attribute_xml(lun, "fields", fname, 'int32', &
shape(fields), dbytes, hash_type, ddigest, cbytes, cdigest) shape(fields), dbytes, hash_type, ddigest, cbytes, cdigest)
write(fname,"(a,'.bin')") "metablock_bounds" write(fname,"(a,'.bin')") "metablock_bounds"
dbytes = size(bounds, kind=8) * kind(bounds) dbytes = size(bounds, kind=8) * kind(bounds)
call write_binary_xml(dname, fname, c_loc(bounds), & call write_binary_xml(dname, fname, c_loc(bounds), dbytes, &
dbytes, hash_type, ddigest, cbytes, cdigest) kind(bounds), hash_type, ddigest, cbytes, cdigest)
call write_attribute_xml(lun, "bounds", fname, 'float64', & call write_attribute_xml(lun, "bounds", fname, 'float64', &
shape(bounds), dbytes, hash_type, ddigest, cbytes, cdigest) shape(bounds), dbytes, hash_type, ddigest, cbytes, cdigest)
@ -2610,8 +2641,8 @@ module io
write(fname,"(a,'_',a,'_',i6.6,'.bin')") "datablock", "ids", nproc write(fname,"(a,'_',a,'_',i6.6,'.bin')") "datablock", "ids", nproc
dbytes = size(ids, kind=8) * kind(ids) dbytes = size(ids, kind=8) * kind(ids)
call write_binary_xml(dname, fname, c_loc(ids), & call write_binary_xml(dname, fname, c_loc(ids), dbytes, &
dbytes, hash_type, ddigest, cbytes, cdigest) kind(ids), hash_type, ddigest, cbytes, cdigest)
call write_attribute_xml(lun, "ids", fname, 'int32', shape(ids), & call write_attribute_xml(lun, "ids", fname, 'int32', shape(ids), &
dbytes, hash_type, ddigest, cbytes, cdigest) dbytes, hash_type, ddigest, cbytes, cdigest)
@ -2630,8 +2661,8 @@ module io
write(fname,"(a,'_',a,'_',i6.6,'.bin')") "datablock", & write(fname,"(a,'_',a,'_',i6.6,'.bin')") "datablock", &
trim(pvars(p)), nproc trim(pvars(p)), nproc
call write_binary_xml(dname, fname, c_loc(array), & call write_binary_xml(dname, fname, c_loc(array), dbytes, &
dbytes, hash_type, ddigest, cbytes, cdigest) kind(array), hash_type, ddigest, cbytes, cdigest)
call write_attribute_xml(lun, pvars(p), fname, 'float64', & call write_attribute_xml(lun, pvars(p), fname, 'float64', &
shape(array), dbytes, hash_type, ddigest, cbytes, cdigest) shape(array), dbytes, hash_type, ddigest, cbytes, cdigest)
end do end do
@ -2860,6 +2891,12 @@ module io
' compressed_digest="' // trim(adjustl(str)) // '"' ' compressed_digest="' // trim(adjustl(str)) // '"'
end if end if
end if end if
select case(data_filter)
case(1)
string = trim(string) // ' data_filter="shuffle"'
case(2)
string = trim(string) // ' data_filter="bytedelta"'
end select
end if end if
end if end if
string = trim(string) // '>' // trim(adjustl(fname)) // '</Attribute>' string = trim(string) // '>' // trim(adjustl(fname)) // '</Attribute>'
@ -2882,6 +2919,7 @@ module io
! path, name - the path and name where the array should be written to; ! path, name - the path and name where the array should be written to;
! array_ptr - the pointer to the array to store; ! array_ptr - the pointer to the array to store;
! array_bytes - the size of the array in bytes; ! array_bytes - the size of the array in bytes;
! item_size - the size of element;
! digest_type - the type of digest to hash the data; ! digest_type - the type of digest to hash the data;
! array_digest - the digest of the original array; ! array_digest - the digest of the original array;
! compressed_bytes - the size of the compressed array in bytes; ! compressed_bytes - the size of the compressed array in bytes;
@ -2889,11 +2927,12 @@ module io
! !
!=============================================================================== !===============================================================================
! !
subroutine write_binary_xml(path, name, array_ptr, array_bytes, digest_type, & subroutine write_binary_xml(path, name, array_ptr, array_bytes, item_size, &
array_digest, compressed_bytes, compressed_digest) digest_type, array_digest, &
compressed_bytes, compressed_digest)
use compression , only : get_compression, compression_bound, compress use compression , only : get_compression, compression_bound, compress
use compression , only : compression_suffix use compression , only : compression_suffix, encode
use hash , only : digest use hash , only : digest
use iso_c_binding, only : c_ptr, c_loc, c_f_pointer use iso_c_binding, only : c_ptr, c_loc, c_f_pointer
@ -2901,6 +2940,7 @@ module io
character(len=*), intent(in) :: path, name character(len=*), intent(in) :: path, name
type(c_ptr) , intent(in) :: array_ptr type(c_ptr) , intent(in) :: array_ptr
integer(kind=4) , intent(in) :: item_size
integer(kind=8) , intent(in) :: array_bytes integer(kind=8) , intent(in) :: array_bytes
integer , intent(in) :: digest_type integer , intent(in) :: digest_type
@ -2915,8 +2955,8 @@ module io
integer(kind=1), dimension(:), pointer :: array integer(kind=1), dimension(:), pointer :: array
integer(kind=1), dimension(:), allocatable, target :: buffer integer(kind=1), dimension(:), allocatable, target :: input, buffer
type(c_ptr) :: buffer_ptr type(c_ptr) :: input_ptr, buffer_ptr
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -2929,10 +2969,12 @@ module io
! !
if (present(compressed_bytes) .and. get_compression() > 0) then if (present(compressed_bytes) .and. get_compression() > 0) then
compressed_bytes = compression_bound(array_bytes) compressed_bytes = compression_bound(array_bytes)
allocate(buffer(compressed_bytes), stat = status) allocate(input(array_bytes), buffer(compressed_bytes), stat = status)
buffer_ptr = c_loc(buffer)
if (status == 0) then if (status == 0) then
call compress(array_ptr, array_bytes, buffer_ptr, compressed_bytes) input_ptr = c_loc(input)
buffer_ptr = c_loc(buffer)
call encode(data_filter, item_size, array_bytes, array_ptr, input)
call compress(input_ptr, array_bytes, buffer_ptr, compressed_bytes)
if (compressed_bytes > 0 .and. compressed_bytes < array_bytes) then if (compressed_bytes > 0 .and. compressed_bytes < array_bytes) then
write(fname,"(a,'/',a,a)") trim(path), trim(name), & write(fname,"(a,'/',a,a)") trim(path), trim(name), &
trim(compression_suffix) trim(compression_suffix)
@ -2948,7 +2990,7 @@ module io
compressed_bytes = 0 compressed_bytes = 0
compressed_digest = 0 compressed_digest = 0
end if end if
deallocate(buffer) deallocate(input, buffer)
end if end if
end if end if

128
sources/problems.F90
View File

@ -116,6 +116,9 @@ module problems
case("turbulence") case("turbulence")
setup_problem => setup_problem_turbulence setup_problem => setup_problem_turbulence
case("stellar_wind", "stellar-wind", "wind")
setup_problem => setup_problem_wind
case default case default
end select end select
@ -2349,6 +2352,131 @@ module problems
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine setup_problem_turbulence end subroutine setup_problem_turbulence
!
!===============================================================================
!
! subroutine SETUP_PROBLEM_WIND:
! -----------------------------
!
! Subroutine sets the initial conditions for the stellar wind problem.
!
! Arguments:
!
! pdata - pointer to the datablock structure of the currently initialized
! block;
!
!===============================================================================
!
subroutine setup_problem_wind(pdata)
use blocks , only : block_data
use coordinates, only : nn => bcells
use coordinates, only : ax, ay
#if NDIMS == 3
use coordinates, only : az
#endif /* NDIMS == 3 */
use equations , only : nv, magnetized, adiabatic_index, csnd, csnd2
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
use equations , only : prim2cons
use parameters , only : get_parameter
implicit none
type(block_data), pointer, intent(inout) :: pdata
real(kind=8), save :: radius = 1.00d-01
real(kind=8), save :: dens = 1.00d+06
real(kind=8), save :: sonic = 3.00d+00
real(kind=8), save :: wind = 5.00d-01
real(kind=8), save :: pres = 5.00d-02
real(kind=8), save :: bamp = 1.00d-03
logical, save :: first = .true.
integer :: i, j, k
real(kind=8) :: r, r2, s, w
real(kind=8), dimension(nv,nn) :: q, u
real(kind=8), dimension(nn) :: x, y
#if NDIMS == 3
real(kind=8), dimension(nn) :: z
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
if (first) then
call get_parameter("radius" , radius)
call get_parameter("density" , dens )
call get_parameter("sonic" , sonic )
call get_parameter("wind_speed", wind )
if (ipr > 0) then
pres = (wind / sonic)**2 * dens / adiabatic_index
else
csnd = wind / sonic
csnd2 = csnd * csnd
end if
if (magnetized) &
call get_parameter("bamp", bamp)
first = .false.
end if
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
#if NDIMS == 3
z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
#endif /* NDIMS == 3 */
#if NDIMS == 3
do k = 1, nn
#else /* NDIMS == 3 */
k = 1
#endif /* NDIMS == 3 */
do j = 1, nn
do i = 1, nn
#if NDIMS == 3
r2 = x(i)**2 + y(j)**2 + z(k)**2
#else /* NDIMS == 3 */
r2 = x(i)**2 + y(j)**2
#endif /* NDIMS == 3 */
r = max(sqrt(r2), 1.0d-16)
s = max(1.0d+00, r / radius)**(NDIMS - 1)
q(idn,i) = dens / s
if (ipr > 0) q(ipr,i) = pres / s**adiabatic_index
q(ivx,i) = wind * x(i) / r
q(ivy,i) = wind * y(j) / r
#if NDIMS == 3
q(ivz,i) = wind * z(k) / r
#else /* NDIMS == 3 */
q(ivz,i) = 0.0d+00
#endif /* NDIMS == 3 */
if (magnetized) then
w = max(1.0d+00, r / radius)**(5.0d-1 * NDIMS)
q(ibx,i) = 0.0d+00
q(iby,i) = 0.0d+00
q(ibz,i) = bamp / w
q(ibp,i) = 0.0d+00
end if
end do
call prim2cons(q(:,:), u(:,:))
pdata%q(:,:,j,k) = q(:,:)
pdata%u(:,:,j,k) = u(:,:)
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
end subroutine setup_problem_wind
!=============================================================================== !===============================================================================
! !

129
sources/shapes.F90
View File

@ -118,6 +118,9 @@ module shapes
case("blast") case("blast")
update_shapes => update_shapes_blast update_shapes => update_shapes_blast
case("stellar_wind", "stellar-wind", "wind")
update_shapes => update_shapes_wind
end select end select
case default case default
@ -490,6 +493,132 @@ module shapes
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine update_shapes_blast end subroutine update_shapes_blast
!
!===============================================================================
!
! subroutine UPDATE_SHAPES_WIND:
! -----------------------------
!
! Subroutine resets the primitive and conserved variables within a defined
! shape for the wind problem.
!
! Arguments:
!
! pdata - pointer to the data block structure of the currently initialized
! block;
! time - time at the moment of update;
! dt - time step since the last update;
!
!===============================================================================
!
subroutine update_shapes_wind(pdata, time, dt)
use blocks , only : block_data
use coordinates , only : nn => bcells
use coordinates , only : ax, ay
#if NDIMS == 3
use coordinates , only : az
#endif /* NDIMS == 3 */
use equations , only : prim2cons
use equations , only : adiabatic_index, csnd, csnd2
use equations , only : nv, magnetized
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
use parameters , only : get_parameter
implicit none
type(block_data), pointer, intent(inout) :: pdata
real(kind=8) , intent(in) :: time, dt
real(kind=8), save :: radius = 1.00d-01
real(kind=8), save :: dens = 1.00d+06
real(kind=8), save :: sonic = 3.00d+00
real(kind=8), save :: wind = 5.00d-01
real(kind=8), save :: pres = 5.00d-02
real(kind=8), save :: bamp = 1.00d-03
logical , save :: first = .true.
!$omp threadprivate(first)
real(kind=8), save :: r2, r
integer :: i, j, k
real(kind=8), dimension(nv,nn) :: q, u
real(kind=8), dimension(nn) :: x, y
#if NDIMS == 3
real(kind=8), dimension(nn) :: z
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
if (first) then
call get_parameter("radius" , radius)
call get_parameter("density" , dens)
call get_parameter("wind_speed", wind)
if (ipr > 0) then
pres = (wind / sonic)**2 * dens / adiabatic_index
else
csnd = wind / sonic
csnd2 = csnd * csnd
end if
if (magnetized) &
call get_parameter("bamp", bamp)
first = .false.
end if
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
#if NDIMS == 3
z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
#endif /* NDIMS == 3 */
#if NDIMS == 3
do k = 1, nn
#else /* NDIMS == 3 */
k = 1
#endif /* NDIMS == 3 */
do j = 1, nn
do i = 1, nn
#if NDIMS == 3
r2 = x(i)**2 + y(j)**2 + z(k)**2
#else /* NDIMS == 3 */
r2 = x(i)**2 + y(j)**2
#endif /* NDIMS == 3 */
r = max(sqrt(r2), 1.0d-16)
if (r <= radius) then
pdata%q(idn,i,j,k) = dens
if (ipr > 0) pdata%q(ipr,i,j,k) = pres
pdata%q(ivx,i,j,k) = wind * x(i) / r
pdata%q(ivy,i,j,k) = wind * y(j) / r
#if NDIMS == 3
pdata%q(ivz,i,j,k) = wind * z(k) / r
#else /* NDIMS == 3 */
pdata%q(ivz,i,j,k) = 0.0d+00
#endif /* NDIMS == 3 */
if (magnetized) then
pdata%q(ibx,i,j,k) = 0.0d+00
pdata%q(iby,i,j,k) = 0.0d+00
pdata%q(ibz,i,j,k) = bamp
pdata%q(ibp,i,j,k) = 0.0d+00
end if
call prim2cons(pdata%q(:,i,j,k), pdata%u(:,i,j,k))
end if
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
end subroutine update_shapes_wind
!=============================================================================== !===============================================================================
! !