fluxes_base.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: fluxes_base.f90 #
5!# #
6!# Copyright (C) 2007-2017 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Manuel Jung <mjung@astrophysik.uni-kiel.de> #
9!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
10!# #
11!# This program is free software; you can redistribute it and/or modify #
12!# it under the terms of the GNU General Public License as published by #
13!# the Free Software Foundation; either version 2 of the License, or (at #
14!# your option) any later version. #
15!# #
16!# This program is distributed in the hope that it will be useful, but #
17!# WITHOUT ANY WARRANTY; without even the implied warranty of #
18!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
19!# NON INFRINGEMENT. See the GNU General Public License for more #
20!# details. #
21!# #
22!# You should have received a copy of the GNU General Public License #
23!# along with this program; if not, write to the Free Software #
24!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
25!# #
26!#############################################################################
27
28!----------------------------------------------------------------------------!
38!----------------------------------------------------------------------------!
47 USE common_dict
48#ifdef PARALLEL
49#ifdef HAVE_MPI_MOD
50 USE mpi
51#endif
52#endif
53 IMPLICIT NONE
54#ifdef PARALLEL
55#ifdef HAVE_MPIF_H
56 include 'mpif.h'
57#endif
58#endif
59 !--------------------------------------------------------------------------!
60 PRIVATE
61 TYPE, ABSTRACT, EXTENDS (logging_base) :: fluxes_base
63 CLASS(reconstruction_base), ALLOCATABLE &
64 :: reconstruction
65 CLASS(marray_base), ALLOCATABLE &
66 :: minwav,maxwav
67 CLASS(marray_compound), POINTER &
68 :: prim =>null(),&
69 cons =>null(),&
70 pfluxes =>null()
73 REAL, DIMENSION(:,:,:,:), POINTER, CONTIGUOUS &
74 :: bxflux,byflux, &
75 bzflux,bxfold, &
76 byfold,bzfold
77 CONTAINS
78 PROCEDURE :: initfluxes
79 procedure(calculatefluxes), DEFERRED :: calculatefluxes
80 PROCEDURE :: calculatefacedata
81 PROCEDURE :: getboundaryflux
82 procedure(finalize), DEFERRED :: finalize
83 PROCEDURE :: finalize_base
84 END TYPE fluxes_base
85
86 abstract INTERFACE
87 SUBROUTINE calculatefluxes(this,Mesh,Physics,pvar,cvar, &
88 xfluxdydz,yfluxdzdx,zfluxdxdy)
90 IMPLICIT NONE
91 CLASS(fluxes_base), INTENT(INOUT) :: this
92 CLASS(mesh_base), INTENT(IN) :: Mesh
93 CLASS(physics_base), INTENT(INOUT) :: Physics
94 CLASS(marray_compound),INTENT(INOUT):: pvar,cvar
95 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%VNUM), &
96 INTENT(OUT) :: xfluxdydz,yfluxdzdx,zfluxdxdy
97 END SUBROUTINE
98 SUBROUTINE finalize(this)
99 IMPORT fluxes_base
100 IMPLICIT NONE
101 CLASS(fluxes_base), INTENT(INOUT) :: this
102 END SUBROUTINE
103 END INTERFACE
104 INTEGER, PARAMETER :: kt = 1
105! INTEGER, PARAMETER :: HLL = 2
106! INTEGER, PARAMETER :: HLLC = 3
107! INTEGER, PARAMETER :: EXACT = 4
108 !--------------------------------------------------------------------------!
109 PUBLIC :: &
110 ! types
111 fluxes_base, &
112 ! constants
113 kt!, HLL, HLLC, EXACT
114 !--------------------------------------------------------------------------!
115
116CONTAINS
117
119 SUBROUTINE initfluxes(this,Mesh,Physics,config,IO,ftype,fname)
120 IMPLICIT NONE
121 !------------------------------------------------------------------------!
122 CLASS(fluxes_base), INTENT(INOUT) :: this
123 CLASS(mesh_base), INTENT(IN) :: Mesh
124 CLASS(physics_base), INTENT(IN) :: Physics
125 TYPE(dict_typ),POINTER :: config,IO
126 INTEGER :: err
127 INTEGER :: ftype
128 CHARACTER(LEN=*) :: fname
129 !------------------------------------------------------------------------!
130 TYPE(dict_typ),POINTER :: IOrec => null()
131 INTEGER :: valwrite,i
132 CHARACTER(LEN=60) :: key
133 !------------------------------------------------------------------------!
134 CALL this%InitLogging(ftype,fname)
135
136 ! check initialization of Mesh and Physics
137 IF (.NOT.mesh%Initialized().OR..NOT.physics%Initialized()) &
138 CALL this%Error("InitFluxes","mesh and/or physics module uninitialized")
139
140 ! allocate memory for all arrays used in fluxes
141 ALLOCATE(this%minwav,this%maxwav, &
142 this%bxflux(mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,2,physics%VNUM), &
143 this%byflux(mesh%KGMIN:mesh%KGMAX,mesh%IGMIN:mesh%IGMAX,2,physics%VNUM), &
144 this%bzflux(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,2,physics%VNUM), &
145 this%bxfold(mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,2,physics%VNUM), &
146 this%byfold(mesh%KGMIN:mesh%KGMAX,mesh%IGMIN:mesh%IGMAX,2,physics%VNUM), &
147 this%bzfold(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,2,physics%VNUM), &
148 stat = err)
149 IF (err.NE.0) THEN
150 CALL this%Error("InitFluxes", "Unable to allocate memory.")
151 END IF
152
153 ! create RANK 1 mesh arrays
154 this%minwav = marray_base(mesh%NDIMS)
155 this%maxwav = marray_base(mesh%NDIMS)
156
157 ! initialize state vectors on cell interfaces
158 CALL physics%new_statevector(this%prim,primitive,mesh%NFACES)
159 CALL physics%new_statevector(this%cons,conservative,mesh%NFACES)
160
161 ! initialize state vector for physical fluxes
162 CALL physics%new_statevector(this%pfluxes,conservative,mesh%NFACES)
163
164 ! print some information
165 CALL this%Info(" FLUXES---> fluxes type " // trim(this%GetName()))
166
167 ! enable output of reconstructed cell face data if requested
168 CALL getattr(config, "output/pstates", valwrite,0)
169 IF(valwrite.EQ.1) THEN
170 DO i=1, physics%VNUM
171 key = trim(physics%pvarname(i)) // "_pstates"
172 CALL setattr(io, trim(key), &
173 this%prim%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,:,i))
174 END DO
175 END IF
176 CALL getattr(config, "output/cstates", valwrite, 0)
177 IF(valwrite.EQ.1) THEN
178 DO i=1, physics%VNUM
179 key = trim(physics%cvarname(i)) // "_cstates"
180 CALL setattr(io, trim(key), &
181 this%cons%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,:,i))
182 END DO
183 END IF
184 CALL getattr(config, "output/pfluxes", valwrite, 0)
185 IF(valwrite.EQ.1) THEN
186 DO i=1, physics%VNUM
187 key = trim(physics%pvarname(i)) // "_pfluxes"
188 CALL setattr(io, trim(key), &
189 this%pfluxes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,:,i))
190 END DO
191 END IF
192
193 CALL getattr(config, "output/wave_speeds", valwrite, 0)
194 IF(valwrite.EQ.1) THEN
195 CALL setattr(io, "minwav", &
196 this%minwav%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1:mesh%NDIMS))
197 CALL setattr(io, "maxwav", &
198 this%maxwav%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1:mesh%NDIMS))
199 END IF
200
201 ! initialize reconstruction modules
202 CALL new_reconstruction(this%Reconstruction,mesh,physics,config,iorec)
203
204 ! add output data of reconstruction modules to IO dictionary
205 IF(ASSOCIATED(iorec)) CALL setattr(io,"reconstruction",iorec)
206
207
208 ! initialize boundary fluxes
209 this%bxflux(:,:,:,:) = 0.
210 this%byflux(:,:,:,:) = 0.
211 this%bzflux(:,:,:,:) = 0.
212 this%bxfold(:,:,:,:) = 0.
213 this%byfold(:,:,:,:) = 0.
214 this%bzfold(:,:,:,:) = 0.
215 END SUBROUTINE initfluxes
216
220 FUNCTION getboundaryflux(this,Mesh,Physics,direction,comm) RESULT(bflux)
221 IMPLICIT NONE
222 !------------------------------------------------------------------------!
223 CLASS(fluxes_base), INTENT(IN) :: this
224 CLASS(mesh_base), INTENT(IN) :: mesh
225 CLASS(physics_base), INTENT(IN) :: physics
226 INTEGER :: direction
227 REAL, DIMENSION(Physics%VNUM) :: bflux
228 INTEGER, OPTIONAL :: comm ! communicator for MPI
229 !------------------------------------------------------------------------!
230#ifdef PARALLEL
231 REAL, DIMENSION(Physics%VNUM) :: bflux_local
232 INTEGER :: sender_rank(1),dest_ranks(1),rank0(1)
233 INTEGER :: dest_comm,union_comm
234 INTEGER :: world_group,dest_group,union_group,sender_group
235 INTEGER :: ierror
236#endif
237 !------------------------------------------------------------------------!
238 INTENT(IN) :: direction
239 !------------------------------------------------------------------------!
240#ifdef PARALLEL
241 IF (PRESENT(comm)) THEN
242 dest_comm = comm
243 ELSE
244 dest_comm = mpi_comm_null
245 END IF
246#endif
247 SELECT CASE(direction)
248 CASE(west) ! western boundary flux
249#ifdef PARALLEL
250 IF (mesh%mycoords(1).EQ.0) THEN
251 bflux_local(:) = sum(sum(this%bxflux(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1,:),1),1)
252 ELSE
253 bflux_local(:) = 0.0
254 END IF
255#else
256 bflux = sum(sum(this%bxflux(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1,:),1),1)
257#endif
258 CASE(east) ! eastern boundary flux
259#ifdef PARALLEL
260 IF (mesh%mycoords(1).EQ.mesh%dims(1)-1) THEN
261 bflux_local(:) = sum(sum(this%bxflux(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2,:),1),1)
262 ELSE
263 bflux_local(:) = 0.0
264 END IF
265#else
266 bflux = sum(sum(this%bxflux(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2,:),1),1)
267#endif
268 CASE(south) ! southern boundary flux
269#ifdef PARALLEL
270 IF (mesh%mycoords(2).EQ.0) THEN
271 bflux_local(:) = sum(sum(this%byflux(mesh%KMIN:mesh%KMAX,mesh%IMIN:mesh%IMAX,1,:),1),1)
272 ELSE
273 bflux_local(:) = 0.0
274 END IF
275#else
276 bflux = sum(sum(this%byflux(mesh%KMIN:mesh%KMAX,mesh%IMIN:mesh%IMAX,1,:),1),1)
277#endif
278 CASE (north) ! northern boundary flux
279#ifdef PARALLEL
280 IF (mesh%mycoords(2).EQ.mesh%dims(2)-1) THEN
281 bflux_local(:) = sum(sum(this%byflux(mesh%KMIN:mesh%KMAX,mesh%IMIN:mesh%IMAX,2,:),1),1)
282 ELSE
283 bflux_local(:) = 0.0
284 END IF
285#else
286 bflux = sum(sum(this%byflux(mesh%KMIN:mesh%KMAX,mesh%IMIN:mesh%IMAX,2,:),1),1)
287#endif
288 CASE (bottom) ! bottom boundary flux
289#ifdef PARALLEL
290 IF (mesh%mycoords(3).EQ.0) THEN
291 bflux_local(:) = sum(sum(this%bzflux(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1,:),1),1)
292 ELSE
293 bflux_local(:) = 0.0
294 END IF
295#else
296 bflux = sum(sum(this%bzflux(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1,:),1),1)
297#endif
298 CASE (top) ! topper boundary flux
299#ifdef PARALLEL
300 IF (mesh%mycoords(3).EQ.mesh%dims(3)-1) THEN
301 bflux_local(:) = sum(sum(this%bzflux(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2,:),1),1)
302 ELSE
303 bflux_local(:) = 0.0
304 END IF
305#else
306 bflux = sum(sum(this%bzflux(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2,:),1),1)
307#endif
308 CASE DEFAULT
309 CALL this%Error("GetBoundaryFlux","wrong direction")
310 END SELECT
311
312#ifdef PARALLEL
313 ! if dest_comm is the world comm, use simpler and faster AllReduce
314 IF(dest_comm.EQ.mpi_comm_world) THEN
315 CALL mpi_allreduce(bflux_local,bflux,physics%VNUM,default_mpi_real, &
316 mpi_sum,mpi_comm_world,ierror)
317 ELSE
318 ! create world group
319 CALL mpi_comm_group(mpi_comm_world,world_group,ierror)
320 ! determine the destination for the result
321 IF(dest_comm.NE.mpi_comm_null) THEN
322 dest_comm = comm
323 ! set destination group
324 CALL mpi_comm_group(dest_comm,dest_group,ierror)
325 ELSE
326 ! if no communicator is given, send the result
327 ! to the rank0 process
328 dest_ranks(1) = 0
329 CALL mpi_group_incl(world_group,1,dest_ranks,dest_group,ierror)
330 END IF
331 ! collect and sum up the result in process with rank 0 with respect to the
332 ! subset of MPI processes at the boundary under consideration
333 IF (mesh%comm_boundaries(direction).NE.mpi_comm_null) THEN
334 CALL mpi_reduce(bflux_local,bflux,physics%VNUM,default_mpi_real, &
335 mpi_sum,0,mesh%comm_boundaries(direction),ierror)
336 ELSE
337 bflux(:) = 0.0
338 END IF
339 ! get sender group
340 rank0(1) = mesh%rank0_boundaries(direction)
341 CALL mpi_group_incl(world_group,1,rank0,sender_group,ierror)
342 ! merge sender with destination
343 CALL mpi_group_union(sender_group,dest_group,union_group,ierror)
344 ! create a communicator for the union group
345 CALL mpi_comm_create(mpi_comm_world,union_group,union_comm,ierror)
346 IF (union_comm.NE.mpi_comm_null) THEN
347 ! get rank of sender in union group
348 rank0(1) = 0
349 CALL mpi_group_translate_ranks(sender_group,1,rank0,union_group,sender_rank,ierror)
350 IF (sender_rank(1).EQ.mpi_undefined) &
351 CALL this%Error("GetBoundaryFlux","sender rank undefined")
352 ! send result to all processes in communicator 'union_comm'
353 CALL mpi_bcast(bflux,physics%VNUM,default_mpi_real,sender_rank(1),union_comm,ierror)
354 ! free union communicator
355 CALL mpi_comm_free(union_comm,ierror)
356 END IF
357 ! free all groups
358 CALL mpi_group_free(union_group,ierror)
359 CALL mpi_group_free(sender_group,ierror)
360 CALL mpi_group_free(world_group,ierror)
361 CALL mpi_group_free(dest_group,ierror)
362 END IF
363
364#endif
365 END FUNCTION getboundaryflux
366
368 SUBROUTINE calculatefacedata(this,Mesh,Physics,pvar,cvar)
370 IMPLICIT NONE
371 !------------------------------------------------------------------------!
372 CLASS(fluxes_base), INTENT(INOUT) :: this
373 CLASS(mesh_base), INTENT(IN) :: Mesh
374 CLASS(physics_base), INTENT(INOUT) :: Physics
375 CLASS(marray_compound),INTENT(INOUT):: pvar,cvar
376 !------------------------------------------------------------------------!
377 ! reconstruct data on cell faces
378 IF (this%Reconstruction%PrimRecon()) THEN
379 CALL this%Reconstruction%CalculateStates(mesh,physics,pvar,this%prim)
380 CALL physics%Convert2Conservative(this%prim,this%cons)
381 ELSE
382 CALL this%Reconstruction%CalculateStates(mesh,physics,cvar,this%cons)
383 CALL physics%Convert2Primitive(this%cons,this%prim)
384 END IF
385
386 ! update the speed of sound on cell faces (non-isotherml physics only)
387 SELECT TYPE(phys => physics)
388 CLASS IS(physics_euler)
389 SELECT TYPE(prim => this%prim)
390 CLASS IS(statevector_euler)
391 CALL phys%UpdateSoundSpeed(prim)
392 END SELECT
393 END SELECT
394
395 ! get minimal & maximal wave speeds on cell interfaces
396 CALL physics%CalculateWaveSpeeds(mesh,this%prim%data5d,this%cons%data5d,this%minwav,this%maxwav)
397 END SUBROUTINE calculatefacedata
398
400 SUBROUTINE finalize_base(this)
401 IMPLICIT NONE
402 !------------------------------------------------------------------------!
403 CLASS(fluxes_base), INTENT(INOUT) :: this
404 !------------------------------------------------------------------------!
405 IF (.NOT.this%Initialized()) &
406 CALL this%Error("CloseFluxes","not initialized")
407 DEALLOCATE(this%minwav,this%maxwav,this%prim,this%cons,this%pfluxes, &
408 this%bxflux,this%byflux,this%bzflux,this%bxfold,this%byfold,this%bzfold)
409 CALL this%Reconstruction%Finalize()
410 DEALLOCATE(this%Reconstruction)
411 END SUBROUTINE finalize_base
412
413END MODULE fluxes_base_mod
Dictionary for generic data types.
Definition: common_dict.f90:61
type(logging_base), save this
base module for numerical flux functions
Definition: fluxes_base.f90:39
subroutine calculatefacedata(this, Mesh, Physics, pvar, cvar)
Calcualtes face data with reconstruction methods (e. g. limiters)
subroutine initfluxes(this, Mesh, Physics, config, IO, ftype, fname)
Initialize Fluxes.
subroutine finalize_base(this)
Destructor.
integer, parameter, public kt
real function, dimension(physics%vnum) getboundaryflux(this, Mesh, Physics, direction, comm)
Get fluxes at boundaries.
Basic fosite module.
subroutine finalize(this)
Destructor of logging_base class.
integer, save default_mpi_real
default real type for MPI
base class for mesh arrays
Definition: marray_base.f90:36
derived class for compound of mesh arrays
basic mesh module
Definition: mesh_base.f90:72
integer, parameter east
named constant for eastern boundary
Definition: mesh_base.f90:101
integer, parameter bottom
named constant for bottom boundary
Definition: mesh_base.f90:101
integer, parameter south
named constant for southern boundary
Definition: mesh_base.f90:101
integer, parameter top
named constant for top boundary
Definition: mesh_base.f90:101
integer, parameter north
named constant for northern boundary
Definition: mesh_base.f90:101
integer, parameter west
named constant for western boundary
Definition: mesh_base.f90:101
Basic physics module.
@, public primitive
@, public conservative
physics module for 1D,2D and 3D non-isothermal Euler equations
constructor for physics class
constructor for reconstruction class
subroutine new_reconstruction(Reconstruction, Mesh, Physics, config, IO)
common data structure
basic mesh array class
Definition: marray_base.f90:69
mesh data structure
Definition: mesh_base.f90:122