boundary_custom.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: boundary_custom.f90 #
5!# #
6!# Copyright (C) 2006-2018 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Jannes Klee <jklee@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
27!----------------------------------------------------------------------------!
35!----------------------------------------------------------------------------!
36MODULE boundary_custom_mod
37 USE boundary_base_mod
38 USE boundary_fixed_mod
39 USE marray_compound_mod
40 USE mesh_base_mod
41 USE physics_base_mod
42 USE common_dict
43 IMPLICIT NONE
44 !--------------------------------------------------------------------------!
45 PRIVATE
46 TYPE, EXTENDS(boundary_fixed) :: boundary_custom
47 INTEGER, DIMENSION(:), ALLOCATABLE :: cbtype
48 REAL, DIMENSION(:,:,:), ALLOCATABLE :: rscale, & !< radial scaling constants
49 invrscale
50 REAL, DIMENSION(:,:,:), POINTER :: radius
51 CONTAINS
52 PROCEDURE :: initboundary_custom
53 PROCEDURE :: finalize
54 PROCEDURE :: setboundarydata
55 PROCEDURE :: setcustomboundaries
56 END TYPE
57 CHARACTER(LEN=32), PARAMETER :: boundcond_name = "custom"
58 !--------------------------------------------------------------------------!
59 ! create bit masks for custom boundary conditions
60 ENUM, BIND(C)
61 ENUMERATOR :: custom_undefined = 0, &
62 custom_nograd = 1, & ! no gradients (default)
63 custom_period = 2, & ! periodic
64 custom_reflect = 3, & ! reflecting
65 custom_reflneg = 4, & ! reflect and change sign
66 custom_extrapol = 5, & ! linear extrapolation
67 custom_fixed = 6, & ! set fixed boundary data
68 custom_logexpol = 7, & ! extrapolation of log values
69 custom_outflow = 8, & ! nograd/reflect depending on flow direction
70 custom_kepler = 9, & ! extrapolate according to Keplers law
71 custom_angkepler = 10 ! same but for angular momentum instead of
72 END ENUM
73 !--------------------------------------------------------------------------!
74 PUBLIC :: &
75 boundary_custom, &
76 ! constants
77 custom_nograd, custom_period, custom_reflect, custom_reflneg, &
78 custom_extrapol, custom_fixed, custom_logexpol, &
79 custom_outflow, custom_kepler, custom_angkepler
80 !--------------------------------------------------------------------------!
81
82CONTAINS
83
85 SUBROUTINE initboundary_custom(this,Mesh,Physics,dir,config)
86 IMPLICIT NONE
87 !------------------------------------------------------------------------!
88 CLASS(boundary_custom), INTENT(INOUT) :: this
89 CLASS(mesh_base), INTENT(IN) :: Mesh
90 CLASS(physics_base), INTENT(IN) :: Physics
91 TYPE(dict_typ), POINTER :: config
92 INTEGER, INTENT(IN) :: dir
93 !------------------------------------------------------------------------!
94 INTEGER :: err = 0
95 INTEGER :: i,j,k
96 !------------------------------------------------------------------------!
97 CALL this%InitBoundary(mesh,physics,custom,boundcond_name,dir,config)
98
99 ! allocate memory for boundary data and mask
100 ALLOCATE(this%cbtype(physics%VNUM),stat=err)
101 IF (err.NE.0) &
102 CALL this%Error("InitBoundary_custom", "Unable to allocate memory.")
103
104 ! allocate memory for boundary data and mask
105 ! this array contains the boundary condition for each primitive variable;
106 ! the user must call SetCustomBoundaries() after initialization to specifiy
107 ! the actual boundary condition for each variable and supply user defined
108 ! data arrays if necessary
109 this%cbtype(:) = custom_nograd
110 END SUBROUTINE initboundary_custom
111
112
114 SUBROUTINE setboundarydata(this,Mesh,Physics,time,pvar)
115 IMPLICIT NONE
116 !------------------------------------------------------------------------!
117 CLASS(boundary_custom), INTENT(INOUT) :: this
118 CLASS(mesh_base), INTENT(IN) :: Mesh
119 CLASS(physics_base), INTENT(IN) :: Physics
120 REAL, INTENT(IN) :: time
121 CLASS(marray_compound), INTENT(INOUT) :: pvar
122 !------------------------------------------------------------------------!
123 INTEGER :: i,j,k,m
124 !------------------------------------------------------------------------!
125 SELECT CASE(this%GetDirection())
126 CASE(west)
127 ! REMARK: vector performance is poor, because the first dimension is short;
128 ! nested do loops seem to be the best, at least in terms run time;
129 ! anything involving array assignments is worse
130 DO m=1,physics%VNUM
131 SELECT CASE(this%cbtype(m))
132 CASE(custom_nograd)
133 DO k=mesh%KMIN,mesh%KMAX
134 DO j=mesh%JMIN,mesh%JMAX
135!NEC$ SHORTLOOP
136 DO i=1,mesh%GINUM
137 pvar%data4d(mesh%IMIN-i,j,k,m) = pvar%data4d(mesh%IMIN,j,k,m)
138 END DO
139 END DO
140 END DO
141 CASE(custom_period)
142 DO k=mesh%KMIN,mesh%KMAX
143 DO j=mesh%JMIN,mesh%JMAX
144!NEC$ SHORTLOOP
145 DO i=1,mesh%GINUM
146 pvar%data4d(mesh%IMIN-i,j,k,m) = pvar%data4d(mesh%IMAX-i+1,j,k,m)
147 END DO
148 END DO
149 END DO
150 CASE(custom_reflect)
151 DO k=mesh%KMIN,mesh%KMAX
152 DO j=mesh%JMIN,mesh%JMAX
153!NEC$ SHORTLOOP
154 DO i=1,mesh%GINUM
155 pvar%data4d(mesh%IMIN-i,j,k,m) = pvar%data4d(mesh%IMIN+i-1,j,k,m)
156 END DO
157 END DO
158 END DO
159 CASE(custom_reflneg)
160 DO k=mesh%KMIN,mesh%KMAX
161 DO j=mesh%JMIN,mesh%JMAX
162!NEC$ SHORTLOOP
163 DO i=1,mesh%GINUM
164 pvar%data4d(mesh%IMIN-i,j,k,m) = -pvar%data4d(mesh%IMIN+i-1,j,k,m)
165 END DO
166 END DO
167 END DO
168 CASE(custom_extrapol)
169 DO k=mesh%KMIN,mesh%KMAX
170 DO j=mesh%JMIN,mesh%JMAX
171!NEC$ SHORTLOOP
172 DO i=1,mesh%GINUM
173 pvar%data4d(mesh%IMIN-i,j,k,m) = (i+1)*pvar%data4d(mesh%IMIN,j,k,m) &
174 - i*pvar%data4d(mesh%IMIN+1,j,k,m)
175 END DO
176 END DO
177 END DO
178 CASE(custom_fixed)
179 DO k=mesh%KMIN,mesh%KMAX
180 DO j=mesh%JMIN,mesh%JMAX
181!NEC$ SHORTLOOP
182 DO i=1,mesh%GINUM
183 pvar%data4d(mesh%IMIN-i,j,k,m) = this%data(i,j,k,m)
184 END DO
185 END DO
186 END DO
187 CASE(custom_logexpol)
188 DO k=mesh%KMIN,mesh%KMAX
189 DO j=mesh%JMIN,mesh%JMAX
190!NEC$ SHORTLOOP
191 DO i=1,mesh%GINUM
192 pvar%data4d(mesh%IMIN-i,j,k,m) = pvar%data4d(mesh%IMIN,j,k,m) &
193 * abs(pvar%data4d(mesh%IMIN,j,k,m) / pvar%data4d(mesh%IMIN+1,j,k,m))**i
194 END DO
195 END DO
196 END DO
197 CASE(custom_outflow)
198 DO k=mesh%KMIN,mesh%KMAX
199 DO j=mesh%JMIN,mesh%JMAX
200 IF (pvar%data4d(mesh%IMIN,j,k,m).GE.0.0 ) THEN
201!NEC$ SHORTLOOP
202 DO i=1,mesh%GINUM
203 ! reflect and flip sign (REFLNEG)
204 pvar%data4d(mesh%IMIN-i,j,k,m) = -pvar%data4d(mesh%IMIN+i-1,j,k,m)
205 END DO
206 ELSE
207!NEC$ SHORTLOOP
208 DO i=1,mesh%GINUM
209 ! no gradient
210 pvar%data4d(mesh%IMIN-i,j,k,m) = pvar%data4d(mesh%IMIN,j,k,m)
211 END DO
212 END IF
213 END DO
214 END DO
215 CASE(custom_kepler)
216 DO k=mesh%KMIN,mesh%KMAX
217 DO j=mesh%JMIN,mesh%JMAX
218!NEC$ SHORTLOOP
219 DO i=1,mesh%GINUM
220 pvar%data4d(mesh%IMIN-i,j,k,m) = (pvar%data4d(mesh%IMIN,j,k,m) &
221 + mesh%radius%bcenter(mesh%IMIN,j,k)*mesh%OMEGA) * this%invRscale(i,j,k) &
222 - mesh%radius%bcenter(mesh%IMIN-i,j,k)*mesh%OMEGA
223 END DO
224 END DO
225 END DO
226 CASE(custom_angkepler)
227 DO k=mesh%KMIN,mesh%KMAX
228 DO j=mesh%JMIN,mesh%JMAX
229!NEC$ SHORTLOOP
230 DO i=1,mesh%GINUM
231 pvar%data4d(mesh%IMIN-i,j,k,m) = (pvar%data4d(mesh%IMIN,j,k,m) &
232 + mesh%radius%bcenter(mesh%IMIN,j,k)**2*mesh%OMEGA) * this%Rscale(i,j,k) &
233 - mesh%radius%bcenter(mesh%IMIN-i,j,k)**2*mesh%OMEGA
234 END DO
235 END DO
236 END DO
237!NEC$ SHORTLOOP
238 DO i=1,mesh%GINUM
239 END DO
240 CASE DEFAULT
241 this%err = ior(custom,int(z'0100'))
242 END SELECT
243 END DO
244 CASE(east)
245 ! REMARK: vector performance is poor, because the first dimension is short;
246 ! nested do loops seem to be the best, at least in terms run time;
247 ! anything involving array assignments is worse
248 DO m=1,physics%VNUM
249 SELECT CASE(this%cbtype(m))
250 CASE(custom_nograd)
251 DO k=mesh%KMIN,mesh%KMAX
252 DO j=mesh%JMIN,mesh%JMAX
253!NEC$ SHORTLOOP
254 DO i=1,mesh%GINUM
255 pvar%data4d(mesh%IMAX+i,j,k,m) = pvar%data4d(mesh%IMAX,j,k,m)
256 END DO
257 END DO
258 END DO
259 CASE(custom_period)
260 DO k=mesh%KMIN,mesh%KMAX
261 DO j=mesh%JMIN,mesh%JMAX
262!NEC$ SHORTLOOP
263 DO i=1,mesh%GINUM
264 pvar%data4d(mesh%IMAX+i,j,k,m) = pvar%data4d(mesh%IMIN+i-1,j,k,m)
265 END DO
266 END DO
267 END DO
268 CASE(custom_reflect)
269 DO k=mesh%KMIN,mesh%KMAX
270 DO j=mesh%JMIN,mesh%JMAX
271!NEC$ SHORTLOOP
272 DO i=1,mesh%GINUM
273 pvar%data4d(mesh%IMAX+i,j,k,m) = pvar%data4d(mesh%IMAX-i+1,j,k,m)
274 END DO
275 END DO
276 END DO
277 CASE(custom_reflneg)
278 DO k=mesh%KMIN,mesh%KMAX
279 DO j=mesh%JMIN,mesh%JMAX
280!NEC$ SHORTLOOP
281 DO i=1,mesh%GINUM
282 pvar%data4d(mesh%IMAX+i,j,k,m) = -pvar%data4d(mesh%IMAX-i+1,j,k,m)
283 END DO
284 END DO
285 END DO
286 CASE(custom_extrapol)
287 DO k=mesh%KMIN,mesh%KMAX
288 DO j=mesh%JMIN,mesh%JMAX
289!NEC$ SHORTLOOP
290 DO i=1,mesh%GINUM
291 pvar%data4d(mesh%IMAX+i,j,k,m) = (i+1)*pvar%data4d(mesh%IMAX,j,k,m) &
292 - i*pvar%data4d(mesh%IMAX-1,j,k,m)
293 END DO
294 END DO
295 END DO
296 CASE(custom_fixed)
297 DO k=mesh%KMIN,mesh%KMAX
298 DO j=mesh%JMIN,mesh%JMAX
299!NEC$ SHORTLOOP
300 DO i=1,mesh%GINUM
301 pvar%data4d(mesh%IMAX+i,j,k,m) = this%data(i,j,k,m)
302 END DO
303 END DO
304 END DO
305 CASE(custom_logexpol)
306 DO k=mesh%KMIN,mesh%KMAX
307 DO j=mesh%JMIN,mesh%JMAX
308!NEC$ SHORTLOOP
309 DO i=1,mesh%GINUM
310 pvar%data4d(mesh%IMAX+i,j,k,m) = pvar%data4d(mesh%IMAX,j,k,m) &
311 * abs(pvar%data4d(mesh%IMAX,j,k,m) / pvar%data4d(mesh%IMAX-1,j,k,m))**i
312 END DO
313 END DO
314 END DO
315 CASE(custom_outflow)
316 DO k=mesh%KMIN,mesh%KMAX
317 DO j=mesh%JMIN,mesh%JMAX
318!NEC$ SHORTLOOP
319 IF (pvar%data4d(mesh%IMAX,j,k,m).LE.0.0) THEN
320 DO i=1,mesh%GINUM
321 pvar%data4d(mesh%IMAX+i,j,k,m) = -pvar%data4d(mesh%IMAX-i+1,j,k,m)
322 END DO
323 ELSE
324 DO i=1,mesh%GINUM
325 pvar%data4d(mesh%IMAX+i,j,k,m) = pvar%data4d(mesh%IMAX,j,k,m)
326 END DO
327 END IF
328 END DO
329 END DO
330 CASE(custom_kepler)
331 DO k=mesh%KMIN,mesh%KMAX
332 DO j=mesh%JMIN,mesh%JMAX
333!NEC$ SHORTLOOP
334 DO i=1,mesh%GINUM
335 pvar%data4d(mesh%IMAX+i,j,k,m) = (pvar%data4d(mesh%IMAX,j,k,m) &
336 + mesh%radius%bcenter(mesh%IMAX,j,k)*mesh%OMEGA) * this%invRscale(i,j,k) &
337 - mesh%radius%bcenter(mesh%IMAX+i,j,k)*mesh%OMEGA
338 END DO
339 END DO
340 END DO
341 CASE(custom_angkepler)
342 DO k=mesh%KMIN,mesh%KMAX
343 DO j=mesh%JMIN,mesh%JMAX
344!NEC$ SHORTLOOP
345 DO i=1,mesh%GINUM
346 pvar%data4d(mesh%IMAX+i,j,k,m) = (pvar%data4d(mesh%IMAX,j,k,m) &
347 + mesh%radius%bcenter(mesh%IMAX,j,k)**2*mesh%OMEGA) * this%Rscale(i,j,k) &
348 - mesh%radius%bcenter(mesh%IMAX+i,j,k)**2*mesh%OMEGA
349 END DO
350 END DO
351 END DO
352 CASE DEFAULT
353 this%err = ior(custom,int(z'0200'))
354 END SELECT
355 END DO
356 CASE(south)
357 DO m=1,physics%VNUM
358 SELECT CASE(this%cbtype(m))
359 CASE(custom_nograd)
360!NEC$ SHORTLOOP
361 DO j=1,mesh%GJNUM
362 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
363 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m)
364 END DO
365 CASE(custom_period)
366!NEC$ SHORTLOOP
367 DO j=1,mesh%GJNUM
368 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
369 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX,m)
370 END DO
371 CASE(custom_reflect)
372!NEC$ SHORTLOOP
373 DO j=1,mesh%GJNUM
374 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
375 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX,m)
376 END DO
377 CASE(custom_reflneg)
378!NEC$ SHORTLOOP
379 DO j=1,mesh%GJNUM
380 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
381 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX,m)
382 END DO
383 CASE(custom_extrapol)
384!NEC$ SHORTLOOP
385 DO j=1,mesh%GJNUM
386 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
387 = (j+1)*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m) &
388 - j*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+1,mesh%KMIN:mesh%KMAX,m)
389 END DO
390 CASE(custom_fixed)
391!NEC$ SHORTLOOP
392 DO j=1,mesh%GJNUM
393 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
394 = this%data(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,m)
395 END DO
396 CASE(custom_logexpol)
397!NEC$ SHORTLOOP
398 DO j=1,mesh%GJNUM
399 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
400 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m) &
401 * abs(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m) &
402 / pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+1,mesh%KMIN:mesh%KMAX,m))**j
403 END DO
404 CASE(custom_outflow)
405! REMARK: this would work for any GJNUM, but with lower performance
406! !NEC$ SHORTLOOP
407! DO j=1,Mesh%GJNUM
408! WHERE (pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN,Mesh%KMIN:Mesh%KMAX,m).GE.0.0 )
409! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN-j,Mesh%KMIN:Mesh%KMAX,m) &
410! = -pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN+j-1,Mesh%KMIN:Mesh%KMAX,m)
411! ELSEWHERE
412! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN-j,Mesh%KMIN:Mesh%KMAX,m) &
413! = pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN,Mesh%KMIN:Mesh%KMAX,m)
414! END WHERE
415! END DO
416 WHERE (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m).GE.0.0 )
417 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMIN,mesh%KMIN:mesh%KMAX,m) &
418 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+1,mesh%KMIN:mesh%KMAX,m)
419 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMIN+1,mesh%KMIN:mesh%KMAX,m) &
420 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m)
421 ELSEWHERE
422 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMIN,mesh%KMIN:mesh%KMAX,m) &
423 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m)
424 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMIN+1,mesh%KMIN:mesh%KMAX,m) &
425 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m)
426 END WHERE
427 CASE(custom_kepler)
428!NEC$ SHORTLOOP
429 DO j=1,mesh%GJNUM
430 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
431 = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m) &
432 + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX)*mesh%OMEGA)&
433 * this%invRscale(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX) &
434 - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX)*mesh%OMEGA
435 END DO
436 CASE(custom_angkepler)
437!NEC$ SHORTLOOP
438 DO j=1,mesh%GJNUM
439 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
440 = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m) &
441 + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA)&
442 * this%Rscale(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX) &
443 - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA
444 END DO
445 CASE DEFAULT
446 this%err = ior(custom,int(z'0300'))
447 END SELECT
448 END DO
449 CASE(north)
450 DO m=1,physics%VNUM
451 SELECT CASE(this%cbtype(m))
452 CASE(custom_nograd)
453!NEC$ SHORTLOOP
454 DO j=1,mesh%GJNUM
455 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
456 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
457 END DO
458 CASE(custom_period)
459!NEC$ SHORTLOOP
460 DO j=1,mesh%GJNUM
461 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
462 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX,m)
463 END DO
464 CASE(custom_reflect)
465!NEC$ SHORTLOOP
466 DO j=1,mesh%GJNUM
467 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
468 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX,m)
469 END DO
470 CASE(custom_reflneg)
471!NEC$ SHORTLOOP
472 DO j=1,mesh%GJNUM
473 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
474 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX,m)
475 END DO
476 CASE(custom_extrapol)
477!NEC$ SHORTLOOP
478 DO j=1,mesh%GJNUM
479 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
480 = (j+1)*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
481 - j*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-1,mesh%KMIN:mesh%KMAX,m)
482 END DO
483 CASE(custom_fixed)
484!NEC$ SHORTLOOP
485 DO j=1,mesh%GJNUM
486 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
487 = this%data(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,m)
488 END DO
489 CASE(custom_logexpol)
490!NEC$ SHORTLOOP
491 DO j=1,mesh%GJNUM
492 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
493 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
494 * abs(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
495 / pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-1,mesh%KMIN:mesh%KMAX,m))**j
496 END DO
497 CASE(custom_outflow)
498! REMARK: this would work for any GJNUM, but with lower performance
499! !NEC$ SHORTLOOP
500! DO j=1,Mesh%GJNUM
501! WHERE (pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m).LE.0.0 )
502! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMAX+j,Mesh%KMIN:Mesh%KMAX,m) &
503! = -pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMAX-j+1,Mesh%KMIN:Mesh%KMAX,m)
504! ELSEWHERE
505! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMAX+j,Mesh%KMIN:Mesh%KMAX,m) &
506! = pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m)
507! END WHERE
508! END DO
509 WHERE (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m).LE.0.0 )
510 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMAX,mesh%KMIN:mesh%KMAX,m) &
511 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-1,mesh%KMIN:mesh%KMAX,m)
512 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMAX-1,mesh%KMIN:mesh%KMAX,m) &
513 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
514 ELSEWHERE
515 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMAX,mesh%KMIN:mesh%KMAX,m) &
516 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
517 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMAX-1,mesh%KMIN:mesh%KMAX,m) &
518 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
519 END WHERE
520 CASE(custom_kepler)
521!NEC$ SHORTLOOP
522 DO j=1,mesh%GJNUM
523 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
524 = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
525 + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX)*mesh%OMEGA)&
526 * this%invRscale(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX) &
527 - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX)*mesh%OMEGA
528 END DO
529 CASE(custom_angkepler)
530!NEC$ SHORTLOOP
531 DO j=1,mesh%GJNUM
532 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
533 = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
534 + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA)&
535 * this%Rscale(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX) &
536 - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA
537 END DO
538 CASE DEFAULT
539 this%err = ior(custom,int(z'0400'))
540 END SELECT
541 END DO
542 CASE(bottom)
543 DO m=1,physics%VNUM
544 SELECT CASE(this%cbtype(m))
545 CASE(custom_nograd)
546!NEC$ SHORTLOOP
547 DO k=1,mesh%GKNUM
548 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
549 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m)
550 END DO
551 CASE(custom_period)
552!NEC$ SHORTLOOP
553 DO k=1,mesh%GKNUM
554 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
555 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-k+1,m)
556 END DO
557 CASE(custom_reflect)
558!NEC$ SHORTLOOP
559 DO k=1,mesh%GKNUM
560 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
561 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+k-1,m)
562 END DO
563 CASE(custom_reflneg)
564!NEC$ SHORTLOOP
565 DO k=1,mesh%GKNUM
566 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
567 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+k-1,m)
568 END DO
569 CASE(custom_extrapol)
570!NEC$ SHORTLOOP
571 DO k=1,mesh%GKNUM
572 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
573 = (k+1)*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m) &
574 - k*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+1,m)
575 END DO
576 CASE(custom_fixed)
577!NEC$ SHORTLOOP
578 DO k=1,mesh%GKNUM
579 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
580 = this%data(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,m)
581 END DO
582 CASE(custom_logexpol)
583!NEC$ SHORTLOOP
584 DO k=1,mesh%GKNUM
585 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
586 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m) &
587 * abs(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m) &
588 / pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+1,m))**k
589 END DO
590 CASE(custom_outflow)
591! REMARK: this would work for any GKNUM, but with lower performance
592! !NEC$ SHORTLOOP
593! DO k=1,Mesh%GKNUM
594! WHERE (pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN,m).GE.0.0 )
595! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN-k,m) &
596! = -pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN+k-1,m)
597! ELSEWHERE
598! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN-k,m) &
599! = pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN,m)
600! END WHERE
601! END DO
602!NEC$ SHORTLOOP
603 DO k=1,mesh%GKNUM
604 WHERE (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m).GE.0.0 )
605 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMIN,m) &
606 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+1,m)
607 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMIN+1,m) &
608 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m)
609 ELSEWHERE
610 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMIN,m) &
611 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m)
612 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMIN+1,m) &
613 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m)
614 END WHERE
615 END DO
616 CASE(custom_kepler)
617!NEC$ SHORTLOOP
618 DO k=1,mesh%GKNUM
619 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
620 = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m) &
621 + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN)*mesh%OMEGA)&
622 * this%invRscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k) &
623 - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k)*mesh%OMEGA
624 END DO
625 CASE(custom_angkepler)
626!NEC$ SHORTLOOP
627 DO k=1,mesh%GKNUM
628 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
629 = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m) &
630 + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN)**2*mesh%OMEGA)&
631 * this%invRscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k) &
632 - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k)**2*mesh%OMEGA
633 END DO
634 CASE DEFAULT
635 this%err = ior(custom,int(z'0500'))
636 END SELECT
637 END DO
638 CASE(top)
639 DO m=1,physics%VNUM
640 SELECT CASE(this%cbtype(m))
641 CASE(custom_nograd)
642!NEC$ SHORTLOOP
643 DO k=1,mesh%GKNUM
644 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
645 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m)
646 END DO
647 CASE(custom_period)
648!NEC$ SHORTLOOP
649 DO k=1,mesh%GKNUM
650 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
651 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+k-1,m)
652 END DO
653 CASE(custom_reflect)
654!NEC$ SHORTLOOP
655 DO k=1,mesh%GKNUM
656 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
657 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-k+1,m)
658 END DO
659 CASE(custom_reflneg)
660!NEC$ SHORTLOOP
661 DO k=1,mesh%GKNUM
662 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
663 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-k+1,m)
664 END DO
665 CASE(custom_extrapol)
666!NEC$ SHORTLOOP
667 DO k=1,mesh%GKNUM
668 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
669 = (k+1)*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m) &
670 - k*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-1,m)
671 END DO
672 CASE(custom_fixed)
673!NEC$ SHORTLOOP
674 DO k=1,mesh%GKNUM
675 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
676 = this%data(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,m)
677 END DO
678 CASE(custom_logexpol)
679!NEC$ SHORTLOOP
680 DO k=1,mesh%GKNUM
681 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
682 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m) &
683 * abs(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m) &
684 / pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-1,m))**k
685 END DO
686 CASE(custom_outflow)
687! REMARK: this would work for any GKNUM, but with lower performance
688! !NEC$ SHORTLOOP
689! DO k=1,Mesh%GKNUM
690! WHERE (pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMAX,m).LE.0.0)
691! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMAX+k,m) &
692! = -pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMAX-k+1,m)
693! ELSEWHERE
694! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMAX+k,m) &
695! = pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMAX,m)
696! END WHERE
697!NEC$ SHORTLOOP
698 DO k=1,mesh%GKNUM
699 WHERE (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m).LE.0.0)
700 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMAX,m) &
701 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-1,m)
702 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMAX-1,m) &
703 = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m)
704 ELSEWHERE
705 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMAX,m) &
706 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m)
707 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMAX-1,m) &
708 = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m)
709 END WHERE
710 END DO
711 CASE(custom_kepler)
712!NEC$ SHORTLOOP
713 DO k=1,mesh%GKNUM
714 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
715 = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m) &
716 + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX)*mesh%OMEGA)&
717 * this%invRscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k) &
718 - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k)*mesh%OMEGA
719 END DO
720 CASE(custom_angkepler)
721!NEC$ SHORTLOOP
722 DO k=1,mesh%GKNUM
723 pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
724 = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m) &
725 + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX)**2*mesh%OMEGA)&
726 * this%invRscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k) &
727 - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k)**2*mesh%OMEGA
728 END DO
729 CASE DEFAULT
730 this%err = ior(custom,int(z'0600'))
731 END SELECT
732 END DO
733 END SELECT
734 END SUBROUTINE setboundarydata
735
736 SUBROUTINE setcustomboundaries(this,Mesh,Physics,cbtype,kepler_radius)
737 IMPLICIT NONE
738 !------------------------------------------------------------------------!
739 CLASS(boundary_custom), INTENT(INOUT) :: this
740 CLASS(mesh_base), INTENT(IN) :: Mesh
741 CLASS(physics_base), INTENT(IN) :: Physics
742 INTEGER, DIMENSION(1:Physics%VNUM) :: cbtype
743 REAL, DIMENSION(:,:,:), POINTER, OPTIONAL :: kepler_radius
744 !------------------------------------------------------------------------!
745 INTEGER :: i,j,k,err = 0
746 !------------------------------------------------------------------------!
747 this%cbtype(:) = cbtype(:)
748 IF (any(this%cbtype(:).EQ.custom_fixed)) THEN
749 SELECT CASE(this%GetDirection())
750 CASE(west,east)
751 ALLOCATE(this%data(mesh%GINUM,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,physics%VNUM), &
752 stat=err)
753 CASE(south,north)
754 ALLOCATE(this%data(mesh%IMIN:mesh%IMAX,mesh%GJNUM,mesh%KMIN:mesh%KMAX,physics%VNUM), &
755 stat=err)
756 CASE(bottom,top)
757 ALLOCATE(this%data(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%GKNUM,physics%VNUM), &
758 stat=err)
759 END SELECT
760 IF (err.NE.0) CALL this%Error("SetCustomBoundaries","Unable to allocate memory.")
761 this%data(:,:,:,:) = 0.0
762 END IF
763
764 IF (any(this%cbtype(:).EQ.custom_outflow)) THEN
765 SELECT CASE(this%GetDirection())
766 CASE(west,east)
767 IF (mesh%GINUM.NE.2) CALL this%Error("SetCustomBoundaries","GINUM must be 2 for outflow boundaries")
768 CASE(south,north)
769 IF (mesh%GJNUM.NE.2) CALL this%Error("SetCustomBoundaries","GJNUM must be 2 for outflow boundaries")
770 CASE(bottom,top)
771 IF (mesh%GKNUM.NE.2) CALL this%Error("SetCustomBoundaries","GKNUM must be 2 for outflow boundaries")
772 END SELECT
773 END IF
774
775 IF (any(this%cbtype(:).EQ.custom_kepler)) THEN
776 SELECT CASE(this%GetDirection())
777 CASE(west,east)
778 ALLOCATE(this%invRscale(mesh%GINUM,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
779 stat=err)
780 CASE(south,north)
781 ALLOCATE(this%invRscale(mesh%IMIN:mesh%IMAX,mesh%GJNUM,mesh%KMIN:mesh%KMAX), &
782 stat=err)
783 CASE(bottom,top)
784 ALLOCATE(this%invRscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%GKNUM), &
785 stat=err)
786 END SELECT
787 IF (err.NE.0) CALL this%Error("SetCustomBoundaries","Unable to allocate memory.")
788
789 IF (PRESENT(kepler_radius)) THEN
790 this%radius => kepler_radius
791 IF (mesh%OMEGA.GT.tiny(mesh%OMEGA)) &
792 CALL this%Warning("SetCustomBoundaries", &
793 "user supplied radius and rotating frame may yield unexpected/wrong results")
794 ELSE
795 this%radius => mesh%radius%bcenter
796 END IF
797
798 SELECT CASE(this%GetDirection())
799 CASE(west)
800 FORALL (i=1:mesh%GINUM,j=mesh%JMIN:mesh%JMAX,k=mesh%KMIN:mesh%KMAX) &
801 this%invRscale(i,j,k) = this%radius(mesh%IMIN,j,k) / this%radius(mesh%IMIN-i,j,k)
802 CASE(east)
803 FORALL (i=1:mesh%GINUM,j=mesh%JMIN:mesh%JMAX,k=mesh%KMIN:mesh%KMAX) &
804 this%invRscale(i,j,k) = this%radius(mesh%IMAX,j,k) / this%radius(mesh%IMAX+i,j,k)
805 CASE(south)
806 FORALL (i=mesh%IMIN:mesh%IMAX,j=1:mesh%GjNUM,k=mesh%KMIN:mesh%KMAX) &
807 this%invRscale(i,j,k) = this%radius(i,mesh%JMIN,k) / this%radius(i,mesh%JMIN-j,k)
808 CASE(north)
809 FORALL (i=mesh%IMIN:mesh%IMAX,j=1:mesh%GjNUM,k=mesh%KMIN:mesh%KMAX) &
810 this%invRscale(i,j,k) = this%radius(i,mesh%JMAX,k) / this%radius(i,mesh%JMAX+j,k)
811 CASE(bottom)
812 FORALL (i=mesh%IMIN:mesh%IMAX,j=mesh%JMIN:mesh%JMAX,k=1:mesh%GKNUM) &
813 this%invRscale(i,j,k) = this%radius(i,j,mesh%KMIN) / this%radius(i,j,mesh%KMIN-k)
814 CASE(top)
815 FORALL (i=mesh%IMIN:mesh%IMAX,j=mesh%JMIN:mesh%JMAX,k=1:mesh%GKNUM) &
816 this%invRscale(i,j,k) = this%radius(i,j,mesh%KMAX) / this%radius(i,j,mesh%KMAX+k)
817 END SELECT
818 this%invRscale(:,:,:) = sqrt(this%invRscale(:,:,:))
819 END IF
820
821 IF (any(this%cbtype(:).EQ.custom_angkepler)) THEN
822 SELECT CASE(this%GetDirection())
823 CASE(west,east)
824 ALLOCATE(this%Rscale(mesh%GINUM,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
825 stat=err)
826 CASE(south,north)
827 ALLOCATE(this%Rscale(mesh%IMIN:mesh%IMAX,mesh%GJNUM,mesh%KMIN:mesh%KMAX), &
828 stat=err)
829 CASE(bottom,top)
830 ALLOCATE(this%Rscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%GKNUM), &
831 stat=err)
832 END SELECT
833
834 IF (PRESENT(kepler_radius)) THEN
835 this%radius => kepler_radius
836 IF (mesh%OMEGA.GT.tiny(mesh%OMEGA)) &
837 CALL this%Warning("SetCustomBoundaries", &
838 "user supplied radius and rotating frame may yield unexpected/wrong results")
839 ELSE
840 this%radius => mesh%radius%bcenter
841 END IF
842
843 IF (err.NE.0) CALL this%Error("SetCustomBoundaries","Unable to allocate memory.")
844 SELECT CASE(this%GetDirection())
845 CASE(west)
846 FORALL (i=1:mesh%GINUM,j=mesh%JMIN:mesh%JMAX,k=mesh%KMIN:mesh%KMAX) &
847 this%Rscale(i,j,k) = this%radius(mesh%IMIN-i,j,k) / this%radius(mesh%IMIN,j,k)
848 CASE(east)
849 FORALL (i=1:mesh%GINUM,j=mesh%JMIN:mesh%JMAX,k=mesh%KMIN:mesh%KMAX) &
850 this%Rscale(i,j,k) = this%radius(mesh%IMAX+i,j,k) / this%radius(mesh%IMAX,j,k)
851 CASE(south)
852 FORALL (i=mesh%IMIN:mesh%IMAX,j=1:mesh%GjNUM,k=mesh%KMIN:mesh%KMAX) &
853 this%Rscale(i,j,k) = this%radius(i,mesh%JMIN-j,k) / this%radius(i,mesh%JMIN,k)
854 CASE(north)
855 FORALL (i=mesh%IMIN:mesh%IMAX,j=1:mesh%GjNUM,k=mesh%KMIN:mesh%KMAX) &
856 this%Rscale(i,j,k) = this%radius(i,mesh%JMAX+j,k) / this%radius(i,mesh%JMAX,k)
857 CASE(bottom)
858 FORALL (i=mesh%IMIN:mesh%IMAX,j=mesh%JMIN:mesh%JMAX,k=1:mesh%GKNUM) &
859 this%Rscale(i,j,k) = this%radius(i,j,mesh%KMIN-k) / this%radius(i,j,mesh%KMIN)
860 CASE(top)
861 FORALL (i=mesh%IMIN:mesh%IMAX,j=mesh%JMIN:mesh%JMAX,k=1:mesh%GKNUM) &
862 this%Rscale(i,j,k) = this%radius(i,j,mesh%KMAX+k) / this%radius(i,j,mesh%KMAX)
863 END SELECT
864 this%Rscale(:,:,:) = sqrt(this%Rscale(:,:,:))
865 END IF
866 END SUBROUTINE
867
868 SUBROUTINE finalize(this)
869 IMPLICIT NONE
870 !------------------------------------------------------------------------!
871 CLASS(boundary_custom), INTENT(INOUT) :: this
872 !------------------------------------------------------------------------!
873 IF (ALLOCATED(this%data)) DEALLOCATE(this%data)
874 IF (ALLOCATED(this%cbtype)) DEALLOCATE(this%cbtype)
875 IF (ALLOCATED(this%Rscale)) DEALLOCATE(this%Rscale)
876 IF (ALLOCATED(this%invRscale)) DEALLOCATE(this%invRscale)
877
878 CALL this%Finalize_base()
879 END SUBROUTINE finalize
880END MODULE boundary_custom_mod
boundary_base_mod
Definition: boundary_base.f90:40
boundary_base_mod::custom
integer, parameter custom
user defined
Definition: boundary_base.f90:110
boundary_custom_mod
Boundary module for custom conditions.
Definition: boundary_custom.f90:36
boundary_custom_mod::custom_reflneg
@, public custom_reflneg
Definition: boundary_custom.f90:61
boundary_custom_mod::custom_extrapol
@, public custom_extrapol
Definition: boundary_custom.f90:61
boundary_custom_mod::setboundarydata
subroutine setboundarydata(this, Mesh, Physics, time, pvar)
Applies the custom boundary conditions.
Definition: boundary_custom.f90:115
boundary_custom_mod::custom_outflow
@, public custom_outflow
Definition: boundary_custom.f90:61
boundary_custom_mod::custom_fixed
@, public custom_fixed
Definition: boundary_custom.f90:61
boundary_custom_mod::custom_undefined
@ custom_undefined
Definition: boundary_custom.f90:61
boundary_custom_mod::finalize
subroutine finalize(this)
Definition: boundary_custom.f90:869
boundary_custom_mod::custom_kepler
@, public custom_kepler
Definition: boundary_custom.f90:61
boundary_custom_mod::boundcond_name
character(len=32), parameter boundcond_name
Definition: boundary_custom.f90:57
boundary_custom_mod::setcustomboundaries
subroutine setcustomboundaries(this, Mesh, Physics, cbtype, kepler_radius)
Definition: boundary_custom.f90:737
boundary_custom_mod::custom_period
@, public custom_period
Definition: boundary_custom.f90:61
boundary_custom_mod::custom_angkepler
@, public custom_angkepler
Definition: boundary_custom.f90:61
boundary_custom_mod::initboundary_custom
subroutine initboundary_custom(this, Mesh, Physics, dir, config)
Constructor for custom boundary conditions.
Definition: boundary_custom.f90:86
boundary_custom_mod::custom_logexpol
@, public custom_logexpol
Definition: boundary_custom.f90:61
boundary_custom_mod::custom_reflect
@, public custom_reflect
Definition: boundary_custom.f90:61
boundary_custom_mod::custom_nograd
@, public custom_nograd
Definition: boundary_custom.f90:61
boundary_fixed_mod
Boundary module for fixed in/outflow conditions.
Definition: boundary_fixed.f90:37
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
mesh_base_mod
basic mesh module
Definition: mesh_base.f90:72
mesh_base_mod::east
integer, parameter east
named constant for eastern boundary
Definition: mesh_base.f90:101
mesh_base_mod::bottom
integer, parameter bottom
named constant for bottom boundary
Definition: mesh_base.f90:101
mesh_base_mod::south
integer, parameter south
named constant for southern boundary
Definition: mesh_base.f90:101
mesh_base_mod::top
integer, parameter top
named constant for top boundary
Definition: mesh_base.f90:101
mesh_base_mod::north
integer, parameter north
named constant for northern boundary
Definition: mesh_base.f90:101
mesh_base_mod::west
integer, parameter west
named constant for western boundary
Definition: mesh_base.f90:101
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
boundary_custom_mod::boundary_custom
Definition: boundary_custom.f90:46
boundary_fixed_mod::boundary_fixed
Definition: boundary_fixed.f90:46
common_dict::dict_typ
Definition: common_dict.f90:123
marray_compound_mod::marray_compound
data types and methods
Definition: marray_compound.f90:46
mesh_base_mod::mesh_base
mesh data structure
Definition: mesh_base.f90:122
physics_base_mod::physics_base
Definition: physics_base.f90:69
Generated by   doxygen 1.9.4