boundary_shearing.f90
Go to the documentation of this file.
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: boundary_shearing.f90 #
5 !# #
6 !# Copyright (C) 2006-2024 #
7 !# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
8 !# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
9 !# #
10 !# This program is free software; you can redistribute it and/or modify #
11 !# it under the terms of the GNU General Public License as published by #
12 !# the Free Software Foundation; either version 2 of the License, or (at #
13 !# your option) any later version. #
14 !# #
15 !# This program is distributed in the hope that it will be useful, but #
16 !# WITHOUT ANY WARRANTY; without even the implied warranty of #
17 !# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
18 !# NON INFRINGEMENT. See the GNU General Public License for more #
19 !# details. #
20 !# #
21 !# You should have received a copy of the GNU General Public License #
22 !# along with this program; if not, write to the Free Software #
23 !# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
24 !# #
25 !#############################################################################
26 !----------------------------------------------------------------------------!
57 !----------------------------------------------------------------------------!
64 USE common_dict
65#ifdef PARALLEL
66#ifdef HAVE_MPI_MOD
67 USE mpi
68#endif
69#endif
70 IMPLICIT NONE
71#ifdef PARALLEL
72#ifdef HAVE_MPIF_H
73 include 'mpif.h'
74#endif
75#endif
76 !--------------------------------------------------------------------------!
77 PRIVATE
78 CHARACTER(LEN=32), PARAMETER :: boundcond_name = "shearing"
79
81 REAL :: velocity_offset
82 REAL, DIMENSION(:), POINTER :: velocity_shift
83 CONTAINS
85 PROCEDURE :: setboundarydata
86 PROCEDURE :: finalize
87 END TYPE
88 !--------------------------------------------------------------------------!
89 PUBLIC :: boundary_shearing
90 !--------------------------------------------------------------------------!
91
92CONTAINS
93
95 SUBROUTINE initboundary_shearing(this,Mesh,Physics,dir,config)
96 IMPLICIT NONE
97 !------------------------------------------------------------------------!
98 CLASS(boundary_shearing), INTENT(INOUT) :: this
99 CLASS(mesh_base), INTENT(IN) :: Mesh
100 CLASS(physics_base), INTENT(IN) :: Physics
101 TYPE(dict_typ), POINTER :: config
102 INTEGER, INTENT(IN) :: dir
103 !------------------------------------------------------------------------!
104 INTEGER :: err, l
105 !------------------------------------------------------------------------!
106 CALL this%InitBoundary(mesh,physics,shearing,boundcond_name,dir,config)
107
108 ALLOCATE( &
109 this%velocity_shift(physics%VNUM+physics%PNUM), &
110 stat = err)
111 IF (err.NE.0) THEN
112 CALL this%Error("boundary_shearing::InitBoundaryg", "Unable to allocate memory.")
113 END IF
114
115 DO l=1,physics%VNUM + physics%PNUM
116 ! this part for WEST-EAST shear
117 IF (mesh%shear_dir.EQ.2) THEN
118 IF (l.EQ.physics%YVELOCITY) THEN
119 SELECT CASE(mesh%fargo%GetType())
120 CASE(0) ! fargo disabled
121 this%velocity_shift(l) = mesh%Q*mesh%OMEGA*(mesh%xmax-mesh%xmin)
122 CASE(3) ! shearing sheet/box fargo enabled
123 this%velocity_shift(l) = 0.0
124 CASE DEFAULT
125 CALL this%Error("boundary_shearing::InitBoundary", &
126 "Shearing boundaries are only compatible without Fargo or with Fargo type 3 (shearing box).")
127 END SELECT
128 ELSE
129 this%velocity_shift(l) = 0.0
130 END IF
131 ! this part for SOUTH-NORTH shear
132 ELSE IF (mesh%shear_dir.EQ.1) THEN
133 IF (l.EQ.physics%XVELOCITY) THEN
134 SELECT CASE(mesh%fargo%GetType())
135 CASE(0) ! fargo disabled
136 this%velocity_shift(l) = mesh%Q*mesh%OMEGA*(mesh%ymax-mesh%ymin)
137 CASE(3) ! shearing sheet/box fargo enabled
138 this%velocity_shift(l) = 0.0
139 CASE DEFAULT
140 CALL this%Error("boundary_shearing::InitBoundary", &
141 "Shearing boundaries are only compatible without Fargo or with Fargo type 3 (shearing box).")
142 END SELECT
143 ELSE
144 this%velocity_shift(l) = 0.0
145 END IF
146 ELSE
147 CALL this%Error("boundary_shearing::InitBoundary", &
148 "Shearing boundaries in top/south direction not allowed, yet.")
149 END IF
150 END DO
151 this%velocity_offset = mesh%Q*mesh%OMEGA*(mesh%xmax-mesh%xmin)/mesh%dy
152
153 END SUBROUTINE initboundary_shearing
154
192 SUBROUTINE setboundarydata(this,Mesh,Physics,time,pvar)
193 IMPLICIT NONE
194 !------------------------------------------------------------------------!
195 CLASS(boundary_shearing), INTENT(INOUT) :: this
196 CLASS(mesh_base), INTENT(IN) :: Mesh
197 CLASS(physics_base), INTENT(IN) :: Physics
198 REAL, INTENT(IN) :: time
199 CLASS(marray_compound), INTENT(INOUT) :: pvar
200 !------------------------------------------------------------------------!
201 INTEGER :: i,j,k,l,intshift
202 REAL :: offremain,offset,offset_tmp
203#ifdef PARALLEL
204 INTEGER :: status(MPI_STATUS_SIZE)
205 INTEGER :: ierror,req(4)
206 CHARACTER(LEN=80) :: str
207 REAL :: mpi_buf(2*Mesh%GNUM)
208#endif
209 !------------------------------------------------------------------------!
210 ! the routine below expect that periodic boundaries are already applied
211#ifndef PARALLEL
212 CALL this%boundary_periodic%SetBoundaryData(mesh,physics,time,pvar)
213#endif
214
215 ! make sure all MPI processes use the same step if domain is decomposed
216 ! along the y-direction (can be different due to round-off errors)
217 offset = -this%velocity_offset*time
218 DO WHILE(offset<0)
219 offset = offset + mesh%JNUM
220 END DO
221
222 DO l=1,physics%VNUM+physics%PNUM
223 ! chose velocity offset for fargo and considered pvar
224 SELECT CASE(this%GetDirection())
225 CASE(west)
226 offset_tmp = offset
227 intshift = floor(offset_tmp)
228 offremain = offset_tmp - intshift
229 DO i=1,mesh%GNUM
230 !------- integral shift ---------------------------------------------!
231 pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,l) = &
232 cshift(pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,l),intshift)
233 !------- residual shift ---------------------------------------------!
234 pvar%data4d(mesh%IMIN-i,mesh%JMAX+1,:,l) = pvar%data4d(mesh%IMIN-i,mesh%JMIN,:,l)
235 DO k=mesh%KMIN,mesh%KMAX
236 DO j=mesh%JMIN,mesh%JMAX
237 pvar%data4d(mesh%IMIN-i,j,k,l) = (1.0 - offremain)*pvar%data4d(mesh%IMIN-i,j,k,l) + &
238 offremain*pvar%data4d(mesh%IMIN-i,j+1,k,l) + this%velocity_shift(l)
239 END DO
240 END DO
241 END DO
242 CASE(east)
243 offset_tmp = -offset
244 intshift = floor(offset_tmp)
245 offremain = offset_tmp - intshift
246 DO i=1,mesh%GNUM
247 !------- integral shift ---------------------------------------------!
248 pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,l) = &
249 cshift(pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,l),intshift)
250 !------- residual shift ---------------------------------------------!
251 pvar%data4d(mesh%IMAX+i,mesh%JMAX+1,:,l) = pvar%data4d(mesh%IMAX+i,mesh%JMIN,:,l)
252 DO k=mesh%KMIN,mesh%KMAX
253 DO j=mesh%JMIN,mesh%JMAX
254 pvar%data4d(mesh%IMAX+i,j,k,l) = (1.0 - offremain)*pvar%data4d(mesh%IMAX+i,j,k,l) + &
255 offremain*pvar%data4d(mesh%IMAX+i,j+1,k,l) - this%velocity_shift(l)
256 END DO
257 END DO
258 END DO
259 CASE(south)
260 offset_tmp = -offset
261 intshift = floor(offset_tmp)
262 offremain = offset_tmp - intshift
263 DO j=1,mesh%GNUM
264 !------- integral shift ---------------------------------------------!
265 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,l) = &
266 cshift(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,l),intshift)
267 !------- residual shift ---------------------------------------------!
268 pvar%data4d(mesh%IMAX+1,mesh%JMIN-j,:,l) = pvar%data4d(mesh%IMIN,mesh%JMIN-j,:,l)
269 DO k=mesh%KMIN,mesh%KMAX
270 DO i=mesh%IMIN,mesh%IMAX
271 pvar%data4d(i,mesh%JMIN-j,k,l) = (1.0 - offremain)*pvar%data4d(i,mesh%JMIN-j,k,l) + &
272 offremain*pvar%data4d(i+1,mesh%JMIN-j,k,l) - this%velocity_shift(l)
273 END DO
274 END DO
275 END DO
276 CASE(north)
277 offset_tmp = offset
278 intshift = floor(offset_tmp)
279 offremain = offset_tmp - intshift
280 DO j=1,mesh%GNUM
281 !------- integral shift ---------------------------------------------!
282 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,l) = &
283 cshift(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,l),intshift)
284 !------- residual shift ---------------------------------------------!
285 pvar%data4d(mesh%IMAX+1,mesh%JMAX+j,:,l) = pvar%data4d(mesh%IMIN,mesh%JMAX+j,:,l)
286 DO k=mesh%KMIN,mesh%KMAX
287 DO i=mesh%IMIN,mesh%IMAX
288 pvar%data4d(i,mesh%JMAX+j,k,l) = (1.0 - offremain)*pvar%data4d(i,mesh%JMAX+j,k,l) + &
289 offremain*pvar%data4d(i+1,mesh%JMAX+j,k,l) + this%velocity_shift(l)
290 END DO
291 END DO
292 END DO
293 CASE(bottom)
294 !CALL this%Error("SetBoundary", "Shearing not supported in BOTTOM direction.")
295 CASE(top)
296 !CALL this%Error("SetBoundary", "Shearing not supported in TOP direction.")
297 END SELECT
298 END DO
299 END SUBROUTINE setboundarydata
300
301
303 SUBROUTINE finalize(this)
304 IMPLICIT NONE
305 !------------------------------------------------------------------------!
306 CLASS(boundary_shearing), INTENT(INOUT) :: this
307 !------------------------------------------------------------------------!
308 DEALLOCATE(this%velocity_shift)
309 CALL this%Finalize_base()
310 END SUBROUTINE finalize
311
312
313
314END MODULE boundary_shearing_mod
integer, parameter shearing
periodic with shear for shearing sheet/box
Boundary module for periodic boundary conditions.
subroutine finalize(this)
Destructor for periodic boundary conditions.
character(len=32), parameter boundcond_name
subroutine setboundarydata(this, Mesh, Physics, time, pvar)
Applies the periodic boundary condition.
Boundary module for a shearingsheet/shearingbox. (see e.g. the standard run )
subroutine initboundary_shearing(this, Mesh, Physics, dir, config)
Constructor for shearing boundary conditions.
Dictionary for generic data types.
Definition: common_dict.f90:61
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.
mesh data structure
Definition: mesh_base.f90:122