From c35d5ccf64cbf4b928cb5935cfe2abf5fd4078f7 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 31 Mar 2018 11:33:06 -0300 Subject: [PATCH 01/10] Improve slightly makefile for makedeps generation and cleaning. File 'makedeps' is now generated automatically if it does not exist or any of source files has been modified. Improve the cleaning of compilation files. Only files created during the compilation are removed. Also if the OBJDIR is empty, it is removed as well. Signed-off-by: Grzegorz Kowal --- build/makefile | 30 ++++++++++++++---------------- 1 file changed, 14 insertions(+), 16 deletions(-) diff --git a/build/makefile b/build/makefile index b1ed814..0ce77fd 100644 --- a/build/makefile +++ b/build/makefile @@ -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 @@ -89,17 +79,18 @@ 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 @@ -113,8 +104,15 @@ $(OBJSDIR): $(DESTDIR): mkdir -p $(DESTDIR) -clean: - rm -rf $(OBJSDIR) $(DESTDIR)/$(name).x +makedeps: $(sources) + ./mkdeps.sh $(SRCSDIR) $(OBJSDIR) > makedeps + +.PHONY: clean + +clean: $(OBJSDIR) + rm -f $(objects) $(modules) $(DESTDIR)/$(name).x + rmdir $(OBJSDIR) + rm -f makedeps #------------------------------------------------------------------------------- # From 461260ed38415c0393196da165691c55864a149e Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 31 Mar 2018 11:44:37 -0300 Subject: [PATCH 02/10] Improve 'make clean' even more. Signed-off-by: Grzegorz Kowal --- build/makefile | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/build/makefile b/build/makefile index 0ce77fd..ce7fbc7 100644 --- a/build/makefile +++ b/build/makefile @@ -109,10 +109,9 @@ makedeps: $(sources) .PHONY: clean -clean: $(OBJSDIR) - rm -f $(objects) $(modules) $(DESTDIR)/$(name).x - rmdir $(OBJSDIR) - rm -f makedeps +clean: + rm -f $(objects) $(modules) $(DESTDIR)/$(name).x makedeps + if [ -d $(OBJSDIR) ]; then rmdir $(OBJSDIR); fi #------------------------------------------------------------------------------- # From 444b78450eaa15e962cecdf59e374cafc8c41516 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 31 Mar 2018 12:01:40 -0300 Subject: [PATCH 03/10] PYTHON: Add two subroutines to read integral (.dat) files. Signed-off-by: Grzegorz Kowal --- python/amun.py | 97 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 97 insertions(+) diff --git a/python/amun.py b/python/amun.py index 863e064..76261e9 100644 --- a/python/amun.py +++ b/python/amun.py @@ -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' From cd4474ad23cb3dc902151f22d0d174e0aa574c48 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 3 May 2018 16:57:46 -0300 Subject: [PATCH 04/10] Makefile: Fix module path for non-GNU compilers. Signed-off-by: Grzegorz Kowal --- build/makefile | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/build/makefile b/build/makefile index ce7fbc7..4793f3d 100644 --- a/build/makefile +++ b/build/makefile @@ -70,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: @@ -94,7 +102,7 @@ $(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) From 5964d507202b5534364a7f02ceec4e81739b0cfa Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 7 Aug 2018 14:13:34 -0300 Subject: [PATCH 05/10] EVOLUTION: Fix unphysical cells after mesh update. Unphysical cells can appear due to various reasons, e.g conservative to primitive variable conversion in recently created block after mesh refinements. Such unphysical cells can be immediately propagated to neighbor blocks through boundary update, blocks which have not been marked as recently updated. In such situation, the unphysical cells are not corrected in these neighbors. This change fixes this by marking all neighbors of blocks created in the mesh update to be verified for unphysical cells too. Signed-off-by: Grzegorz Kowal --- src/blocks.F90 | 2 +- src/evolution.F90 | 21 +++++++++++++++++++-- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/src/blocks.F90 b/src/blocks.F90 index 181fd64..f529178 100644 --- a/src/blocks.F90 +++ b/src/blocks.F90 @@ -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 diff --git a/src/evolution.F90 b/src/evolution.F90 index c38b3f4..dc86c5c 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -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 From 8c0c89d0af7a8a50d71bbe1e1de5b0a921d6bda9 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 7 Aug 2018 14:59:42 -0300 Subject: [PATCH 06/10] INTERPOLATIONS: Avoid infinite loop in interfaces_tvd() and interfaces_mgp(). In case an interpolated cell, which should be positive, is zero or negative, and the sum of derivatives is zero as well, the correction of derivatives can enter into an infinite loop. It shouldn't normally happen, since such a situation is unphysical. Avoid the infinite loop and print a warning about encountered problems. Signed-off-by: Grzegorz Kowal --- src/interpolations.F90 | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/src/interpolations.F90 b/src/interpolations.F90 index 617d486..40af9e2 100644 --- a/src/interpolations.F90 +++ b/src/interpolations.F90 @@ -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 From a3383488ce5fd0ed90917fc3761a36064e2a1132 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 7 Aug 2018 15:15:08 -0300 Subject: [PATCH 07/10] BOUNDARIES: Avoid infinite loop in prolong subroutines. Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 39 ++++++++++++++++++++++++++++++--------- 1 file changed, 30 insertions(+), 9 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index 32986d2..fc49aca 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -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(:) From ae06f4fb733a8ff6d09491c55f78ccdb62fba20f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 7 Aug 2018 15:19:56 -0300 Subject: [PATCH 08/10] MESH: Avoid infinite loop in prolong_block(). Signed-off-by: Grzegorz Kowal --- src/mesh.F90 | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/mesh.F90 b/src/mesh.F90 index 5327e65..f11343e 100644 --- a/src/mesh.F90 +++ b/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(:) From 808c1583ba8fd6372ea4d72d968ee46564f5060d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 21 Aug 2018 22:45:49 -0300 Subject: [PATCH 09/10] EQUATIONS: Improve unphysical cell correction. If there are no enough physical neighbors, the cell cannot be corrected and the code simply stops. Allow the code to continue the execution in such a case by simply correcting the unphysical cells by replacing the negative values of positive variables with a lower bounds, dmin and pmin, for density and pressure, respectively. Signed-off-by: Grzegorz Kowal --- src/equations.F90 | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/equations.F90 b/src/equations.F90 index 6c8c3c4..5c536df 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -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|² ! @@ -1218,13 +1218,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 From 1cefe88de49728eafb0cab4481f80f9f8b93873e Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 21 Aug 2018 23:02:51 -0300 Subject: [PATCH 10/10] EQUATIONS: More control over the unphysical cell correction. Give more control over the unphysical cell correction by allowing to specify the minimum number of physical cells to perform averaging (parameter "npavg") and the maximum distance to be considered from the corrected cell (parameter "ngavg"). These two parameters allow to avoid averaging of extended unphysical regions from two physical points. Signed-off-by: Grzegorz Kowal --- src/equations.F90 | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/src/equations.F90 b/src/equations.F90 index 5c536df..380677e 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -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), &