Merge branch 'master' into reconnection
This commit is contained in:
commit
da7037c517
@ -49,16 +49,6 @@ $(info )
|
||||
$(shell sleep 15)
|
||||
endif
|
||||
|
||||
#-------------------------------------------------------------------------------
|
||||
#
|
||||
# generate object file dependencies
|
||||
#
|
||||
ifneq ($(wildcard makedeps),makedeps)
|
||||
$(info Generating dependencies.)
|
||||
$(shell ./mkdeps.sh $(SRCSDIR) $(OBJSDIR) > makedeps)
|
||||
$(info )
|
||||
endif
|
||||
|
||||
#-------------------------------------------------------------------------------
|
||||
#
|
||||
# check flag conditions
|
||||
@ -80,6 +70,14 @@ FFLAGS += ${CPPPREFIX}-DNDIMS=${NDIMS}
|
||||
#
|
||||
FFLAGS += ${CPPPREFIX}-D${OUTPUT}
|
||||
|
||||
# add module path to compiler options
|
||||
#
|
||||
ifeq ($(COMPILER),GNU)
|
||||
FFLAGS += -J $(OBJSDIR)
|
||||
else
|
||||
FFLAGS += -module $(OBJSDIR)
|
||||
endif
|
||||
|
||||
#-------------------------------------------------------------------------------
|
||||
|
||||
.SUFFIXES:
|
||||
@ -89,21 +87,22 @@ FFLAGS += ${CPPPREFIX}-D${OUTPUT}
|
||||
|
||||
name = amun
|
||||
|
||||
modules = algebra blocks boundaries constants coordinates domains driver \
|
||||
files = algebra blocks boundaries constants coordinates domains driver \
|
||||
equations error evolution gravity integrals interpolations io mesh \
|
||||
mpitools operators parameters problems random refinement schemes \
|
||||
shapes sources timers user_problem utils
|
||||
|
||||
sources := $(addprefix $(SRCSDIR)/,$(addsuffix .F90, $(modules)))
|
||||
objects := $(addprefix $(OBJSDIR)/,$(addsuffix .o, $(modules)))
|
||||
sources := $(addprefix $(SRCSDIR)/,$(addsuffix .F90, $(files)))
|
||||
objects := $(addprefix $(OBJSDIR)/,$(addsuffix .o, $(files)))
|
||||
modules := $(addprefix $(OBJSDIR)/,$(addsuffix .mod, $(files)))
|
||||
|
||||
all: $(name).x
|
||||
|
||||
$(name).x: $(deps) $(objects) | $(DESTDIR)
|
||||
$(name).x: $(objects) | $(DESTDIR)
|
||||
$(LD) $(LDFLAGS) $(objects) $(LIBS) -o $(DESTDIR)/$(name).x
|
||||
|
||||
$(OBJSDIR)/%.o : $(SRCSDIR)/%.F90
|
||||
$(FC) -c $(FFLAGS) -J $(OBJSDIR) $< -o $@
|
||||
$(FC) -c $(FFLAGS) $< -o $@
|
||||
|
||||
$(objects): | $(OBJSDIR)
|
||||
|
||||
@ -113,8 +112,14 @@ $(OBJSDIR):
|
||||
$(DESTDIR):
|
||||
mkdir -p $(DESTDIR)
|
||||
|
||||
makedeps: $(sources)
|
||||
./mkdeps.sh $(SRCSDIR) $(OBJSDIR) > makedeps
|
||||
|
||||
.PHONY: clean
|
||||
|
||||
clean:
|
||||
rm -rf $(OBJSDIR) $(DESTDIR)/$(name).x
|
||||
rm -f $(objects) $(modules) $(DESTDIR)/$(name).x makedeps
|
||||
if [ -d $(OBJSDIR) ]; then rmdir $(OBJSDIR); fi
|
||||
|
||||
#-------------------------------------------------------------------------------
|
||||
#
|
||||
|
||||
@ -418,6 +418,103 @@ def rebin(a, newshape):
|
||||
a = np.repeat(a, newshape[n] / a.shape[n], axis = n)
|
||||
return(a)
|
||||
|
||||
|
||||
def amun_integrals(field, filename, pathlist):
|
||||
'''
|
||||
get_integral: iterate over pathlist and read and merge field values from filename files in the provided paths
|
||||
'''
|
||||
# Initiate the return values with empty array and file number.
|
||||
#
|
||||
vals = array([])
|
||||
num = 1
|
||||
|
||||
# Iterate over all paths provided in the list 'pathlist'.
|
||||
#
|
||||
for path in pathlist:
|
||||
|
||||
# Iterate over all files in the current path.
|
||||
#
|
||||
while True:
|
||||
|
||||
# Generate file name.
|
||||
#
|
||||
dfile = path + '/' + filename + '_' + str(num).zfill(2) + '.dat'
|
||||
|
||||
# Check if the file exists.
|
||||
#
|
||||
if isfile(dfile):
|
||||
|
||||
# Read values from the current integrals file.
|
||||
#
|
||||
lvals = read_integrals(dfile, field)
|
||||
|
||||
# Append to the return array.
|
||||
#
|
||||
vals = append(vals, lvals)
|
||||
|
||||
# Increase the number file.
|
||||
#
|
||||
num = num + 1
|
||||
|
||||
else:
|
||||
|
||||
# File does not exists, so go to the next path.
|
||||
#
|
||||
break
|
||||
|
||||
# Return appended values.
|
||||
#
|
||||
return vals
|
||||
|
||||
def read_integrals(filename, column):
|
||||
'''
|
||||
read_integrals: reads a given column from an integral file.
|
||||
'''
|
||||
# Open the given file and check if it is text file.
|
||||
#
|
||||
f = open(filename, 'r')
|
||||
|
||||
# Read fist line and store it in h, since it will be used to obtain the
|
||||
# column headers.
|
||||
#
|
||||
l = f.readline()
|
||||
h = l
|
||||
|
||||
# Read first line which is not comment in order to determine the number of
|
||||
# columns and store the number of columns in nc. Calculate the column width
|
||||
# and store it in wc.
|
||||
#
|
||||
while l.startswith('#'):
|
||||
l = f.readline()
|
||||
nc = len(l.rsplit())
|
||||
wc = int((len(l) - 9) / (nc - 1))
|
||||
|
||||
# Split header line into a list.
|
||||
#
|
||||
lh = [h[1:9].strip()]
|
||||
for i in range(nc - 1):
|
||||
ib = i * wc + 10
|
||||
ie = ib + wc - 1
|
||||
lh.append(h[ib:ie].strip())
|
||||
|
||||
|
||||
ic = lh.index(column)
|
||||
|
||||
# Read given column.
|
||||
#
|
||||
if (ic > -1):
|
||||
lc = [float(l.split()[ic])]
|
||||
for l in f:
|
||||
lc.append(float(l.split()[ic]))
|
||||
|
||||
# Close the file.
|
||||
#
|
||||
f.close()
|
||||
|
||||
# Return values.
|
||||
#
|
||||
return(array(lc))
|
||||
|
||||
if __name__ == "__main__":
|
||||
fname = './p000030_00000.h5'
|
||||
|
||||
|
||||
@ -343,7 +343,7 @@ module blocks
|
||||
public :: set_last_id, get_last_id, get_mblocks, get_dblocks, get_nleafs
|
||||
public :: set_blocks_update
|
||||
public :: change_blocks_process
|
||||
public :: set_neighbors_refine
|
||||
public :: set_neighbors_refine, set_neighbors_update
|
||||
public :: metablock_set_id, metablock_set_process, metablock_set_level
|
||||
public :: metablock_set_configuration, metablock_set_refinement
|
||||
public :: metablock_set_position, metablock_set_coordinates
|
||||
|
||||
@ -5765,6 +5765,7 @@ module boundaries
|
||||
use coordinates , only : im , jm , km
|
||||
use coordinates , only : faces_dp
|
||||
use equations , only : nv, idn, ipr
|
||||
use error , only : print_warning
|
||||
use interpolations , only : limiter_prol
|
||||
|
||||
! local variables are not implicit by default
|
||||
@ -5842,9 +5843,15 @@ module boundaries
|
||||
dq(3) = limiter_prol(0.5d+00, dql, dqr)
|
||||
|
||||
if (p == idn .or. p == ipr) then
|
||||
do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
if (qn(p,i,j,k) > 0.0d+00) then
|
||||
do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
else
|
||||
call print_warning("boundaries::block_face_prolong" &
|
||||
, "Positive variable is not positive!")
|
||||
dq(:) = 0.0d+00
|
||||
end if
|
||||
end if
|
||||
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
@ -6063,6 +6070,7 @@ module boundaries
|
||||
use coordinates , only : im , jm , km
|
||||
use coordinates , only : edges_dp
|
||||
use equations , only : nv, idn, ipr
|
||||
use error , only : print_warning
|
||||
use interpolations , only : limiter_prol
|
||||
|
||||
! local variables are not implicit by default
|
||||
@ -6155,9 +6163,15 @@ module boundaries
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
if (p == idn .or. p == ipr) then
|
||||
do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
if (qn(p,i,j,k) > 0.0d+00) then
|
||||
do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
else
|
||||
call print_warning("boundaries::block_edge_prolong" &
|
||||
, "Positive variable is not positive!")
|
||||
dq(:) = 0.0d+00
|
||||
end if
|
||||
end if
|
||||
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
@ -6330,6 +6344,7 @@ module boundaries
|
||||
use coordinates , only : im , jm , km
|
||||
use coordinates , only : corners_dp
|
||||
use equations , only : nv, idn, ipr
|
||||
use error , only : print_warning
|
||||
use interpolations , only : limiter_prol
|
||||
|
||||
! local variables are not implicit by default
|
||||
@ -6427,9 +6442,15 @@ module boundaries
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
if (p == idn .or. p == ipr) then
|
||||
do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
if (qn(p,i,j,k) > 0.0d+00) then
|
||||
do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
else
|
||||
call print_warning("boundaries::block_corner_prolong" &
|
||||
, "Positive variable is not positive!")
|
||||
dq(:) = 0.0d+00
|
||||
end if
|
||||
end if
|
||||
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
|
||||
@ -139,7 +139,7 @@ module equations
|
||||
|
||||
! the lower limits for density and pressure to be treated as physical
|
||||
!
|
||||
real(kind=8) , save :: dmin = 1.0d-08, pmin = 1.0d-08
|
||||
real(kind=8) , save :: dmin = 1.0d-16, pmin = 1.0d-16
|
||||
|
||||
! the upper limits for the Lorentz factor and corresponding |v|²
|
||||
!
|
||||
@ -158,9 +158,12 @@ module equations
|
||||
integer , save :: nrmax = 100
|
||||
integer , save :: nrext = 2
|
||||
|
||||
! flags for corrections
|
||||
! flag for unphysical cells correction, the maximum distance of neighbors for
|
||||
! averaging region, and the minimum number of cells for averaging
|
||||
!
|
||||
logical , save :: fix_unphysical_cells = .false.
|
||||
integer , save :: ngavg = 2
|
||||
integer , save :: npavg = 4
|
||||
|
||||
! by default everything is private
|
||||
!
|
||||
@ -868,6 +871,16 @@ module equations
|
||||
fix_unphysical_cells = .false.
|
||||
end select
|
||||
|
||||
! get parameters for unphysical cells correction
|
||||
!
|
||||
call get_parameter_integer("ngavg", ngavg)
|
||||
call get_parameter_integer("npavg", npavg)
|
||||
|
||||
! correct the above parameters to reasonable values
|
||||
!
|
||||
ngavg = max(1, ngavg)
|
||||
npavg = max(2, npavg)
|
||||
|
||||
! print information about the equation module
|
||||
!
|
||||
if (verbose) then
|
||||
@ -879,6 +892,10 @@ module equations
|
||||
end if
|
||||
write (*,"(4x,a20, 3x,'=',1x,a)") "fix unphysical cells" &
|
||||
, trim(unphysical_fix)
|
||||
if (fix_unphysical_cells) then
|
||||
write (*,"(4x,a20, 3x,'=',1x,i4)") "ngavg ", ngavg
|
||||
write (*,"(4x,a20, 3x,'=',1x,i4)") "npavg ", npavg
|
||||
end if
|
||||
|
||||
end if
|
||||
|
||||
@ -1193,7 +1210,7 @@ module equations
|
||||
!
|
||||
np = 0
|
||||
p = 1
|
||||
do while (np <= 2 .and. p <= 4)
|
||||
do while (np <= npavg .and. p <= ngavg)
|
||||
il = max( 1, i - p)
|
||||
iu = min(im, i + p)
|
||||
jl = max( 1, j - p)
|
||||
@ -1208,7 +1225,7 @@ module equations
|
||||
|
||||
! average primitive variables
|
||||
!
|
||||
if (np > 2) then
|
||||
if (np >= npavg) then
|
||||
|
||||
do p = 1, nv
|
||||
q(p,n) = sum(qq(p,il:iu,jl:ju,kl:ku), &
|
||||
@ -1218,13 +1235,17 @@ module equations
|
||||
|
||||
else
|
||||
|
||||
! print error, since no physical cells found for averaging
|
||||
! limit density or pressure to minimum value, since the averaging over
|
||||
! neighbours failed
|
||||
!
|
||||
write(msg,'(a,1x,a)') &
|
||||
"Cannot correct the unphysical cell." &
|
||||
, "Not sufficient number of physical neighbors!"
|
||||
call print_error(loc, trim(msg))
|
||||
stop
|
||||
"Not sufficient number of physical neighbors!" &
|
||||
, "Applying lower bounds for positive variables."
|
||||
call print_warning(loc, trim(msg))
|
||||
|
||||
q(1:nv,n) = qq(1:nv,i,j,k)
|
||||
q(idn ,n) = max(dmin, qq(idn,i,j,k))
|
||||
if (ipr > 0) q(ipr,n) = max(pmin, qq(ipr,i,j,k))
|
||||
|
||||
end if ! not sufficient number of physical cells for averaging
|
||||
|
||||
|
||||
@ -1988,6 +1988,7 @@ module evolution
|
||||
|
||||
! include external procedures
|
||||
!
|
||||
use blocks , only : set_neighbors_update
|
||||
use boundaries , only : boundary_variables
|
||||
use equations , only : update_primitive_variables
|
||||
use equations , only : fix_unphysical_cells, correct_unphysical_states
|
||||
@ -2037,12 +2038,28 @@ module evolution
|
||||
! correct unphysical states if detected
|
||||
!
|
||||
if (fix_unphysical_cells) then
|
||||
|
||||
! if an unphysical cell appeared in a block while updating its primitive
|
||||
! variables it could be propagated to its neighbors through boundary update;
|
||||
! mark all neighbors of such a block to be verified and corrected for
|
||||
! unphysical cells too
|
||||
!
|
||||
pmeta => list_meta
|
||||
do while (associated(pmeta))
|
||||
|
||||
if (pmeta%update) call set_neighbors_update(pmeta)
|
||||
|
||||
pmeta => pmeta%next
|
||||
end do
|
||||
|
||||
! verify and correct, if necessary, unphysical cells in recently updated blocks
|
||||
!
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
pmeta => pdata%meta
|
||||
|
||||
if (pmeta%update) call correct_unphysical_states(pmeta%id &
|
||||
, pdata%q, pdata%u)
|
||||
if (pmeta%update) &
|
||||
call correct_unphysical_states(pmeta%id, pdata%q, pdata%u)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
@ -695,6 +695,7 @@ module interpolations
|
||||
use coordinates , only : im , jm , km
|
||||
use coordinates , only : ib , jb , kb , ie , je , ke
|
||||
use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu
|
||||
use error , only : print_warning
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
@ -768,9 +769,15 @@ module interpolations
|
||||
! variables
|
||||
!
|
||||
if (positive) then
|
||||
do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
if (q(i,j,k) > 0.0d+00) then
|
||||
do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
else
|
||||
call print_warning("interpolations::interfaces_tvd" &
|
||||
, "Positive variable is not positive!")
|
||||
dq(:) = 0.0d+00
|
||||
end if
|
||||
end if
|
||||
|
||||
! interpolate states
|
||||
@ -933,6 +940,7 @@ module interpolations
|
||||
use coordinates , only : im , jm , km
|
||||
use coordinates , only : ib , jb , kb , ie , je , ke
|
||||
use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu
|
||||
use error , only : print_warning
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
@ -1051,9 +1059,15 @@ module interpolations
|
||||
! variables
|
||||
!
|
||||
if (positive) then
|
||||
do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
if (q(i,j,k) > 0.0d+00) then
|
||||
do while (q(i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
else
|
||||
call print_warning("interpolations::interfaces_mgp" &
|
||||
, "Positive variable is not positive!")
|
||||
dq(:) = 0.0d+00
|
||||
end if
|
||||
end if
|
||||
|
||||
! interpolate states
|
||||
|
||||
13
src/mesh.F90
13
src/mesh.F90
@ -1037,6 +1037,7 @@ module mesh
|
||||
use coordinates , only : ng, nh, in, jn, kn, im, jm, km
|
||||
use coordinates , only : ib, ie, jb, je, kb, ke
|
||||
use equations , only : nv, idn, ien
|
||||
use error , only : print_warning
|
||||
use interpolations , only : limiter_prol
|
||||
|
||||
! local variables are not implicit by default
|
||||
@ -1137,9 +1138,15 @@ module mesh
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
if (p == idn .or. p == ien) then
|
||||
do while (pdata%u(p,i,j,k) <= sum(abs(du(1:NDIMS))))
|
||||
du(:) = 0.5d+00 * du(:)
|
||||
end do
|
||||
if (pdata%u(p,i,j,k) > 0.0d+00) then
|
||||
do while (pdata%u(p,i,j,k) <= sum(abs(du(1:NDIMS))))
|
||||
du(:) = 0.5d+00 * du(:)
|
||||
end do
|
||||
else
|
||||
call print_warning("mesh::prolong_block" &
|
||||
, "Positive variable is not positive!")
|
||||
du(:) = 0.0d+00
|
||||
end if
|
||||
end if
|
||||
|
||||
du(:) = 0.5d+00 * du(:)
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user