DISPATCH
boundaries_mod.f90
1 !===============================================================================
2 !> Boundary conditions for centered variables, which have the physical boundary
3 !> in-between mesh points lo and li, or/and in-between mesh points ui and uo.
4 !===============================================================================
6  USE io_mod
7  USE trace_mod
8  USE kinds_mod
9  USE mpi_mod
10  USE bits_mod
11  USE mesh_mod
12  USE index_mod
13  USE boundary_wavekill_mod
14  implicit none
15  private
16  type, public:: boundaries_t
17  integer:: id = 0
18  integer:: bits = 0
19  logical:: unset = .true.
20  class(mesh_t), pointer:: m(:)
21  real(8):: position(3)=0d0, radius=0d0
22  real:: vmax
23  logical(kind=2), pointer:: filled(:,:,:,:)
24  logical:: check_filled=.false.
25  contains
26  procedure:: init
27  procedure:: set
28  procedure:: is_set
29  procedure:: is_clear
30  procedure:: mask
31  procedure:: condition
32  procedure:: symmetric
33  procedure:: antisymmetric
34  procedure:: extrapolate
35  procedure:: log_extrapolate
36  procedure:: constant
37  procedure:: spherical0
38  procedure:: spherical3
39  procedure:: spherical_p_ns
40  procedure:: spherical_ns
41  procedure:: spherical_jagged_d
42  procedure:: spherical_jagged_s0
43  procedure:: spherical_jagged_s
44  procedure:: spherical_jagged_p
45  generic:: spherical => spherical0, spherical3
46  procedure:: spherical_extrapolate
47  procedure:: spherical_constant_entropy
48  procedure, nopass:: sponge => wave_killing
49  end type
50 CONTAINS
51 
52 !===============================================================================
53 !> Initialize the boundary data type, while testing also if the patch contains
54 !> a spherical boundary, or possibly is completely immersed in such a boundary.
55 !===============================================================================
56 SUBROUTINE init (self, mesh, position, radius)
57  class(boundaries_t) :: self
58  class(mesh_t), pointer, optional :: mesh(:)
59  real(8), optional :: position(3), radius
60  !.............................................................................
61  real(8) :: size(3), p(3)
62  integer :: i, j, k
63  logical :: one_inside, all_inside
64  !-----------------------------------------------------------------------------
65  call trace%begin ('boundaries_t%init')
66  if (present(radius)) self%radius = radius
67  if (present(position)) self%position = position
68  if (present(mesh)) then
69  self%m => mesh
70  if (present(position).and.present(radius)) then
71  one_inside = .false.
72  all_inside = .true.
73  size(:) = mesh(:)%s
74  do k=-1,1; do j=-1,1; do i=-1,1
75  p(:) = position(:) + 0.5*size(:)*[i,j,k]
76  if (sum(p**2) < radius**2) then
77  one_inside = .true.
78  else
79  all_inside = .false.
80  end if
81  end do; end do; end do
82  if (one_inside) call self%set (bits%spherical)
83  end if
84  end if
85  self%unset = .false.
86  call trace%end()
87 END SUBROUTINE init
88 
89 !===============================================================================
90 !> Set bits in the status word
91 !===============================================================================
92 SUBROUTINE set (self, bits)
93  class(boundaries_t) :: self
94  integer :: bits
95  !-----------------------------------------------------------------------------
96  self%bits = ior(self%bits, bits)
97 END SUBROUTINE set
98 
99 !===============================================================================
100 !> Check bits in the status word
101 !===============================================================================
102 FUNCTION is_set (self, bits) RESULT (out)
103  class(boundaries_t) :: self
104  integer :: bits
105  logical :: out
106  !-----------------------------------------------------------------------------
107  out = (iand(self%bits, bits) /=0)
108 END FUNCTION is_set
109 
110 !===============================================================================
111 !> Check bits in the status word
112 !===============================================================================
113 FUNCTION is_clear (self, bits) RESULT (out)
114  class(boundaries_t) :: self
115  integer :: bits
116  logical :: out
117  !-----------------------------------------------------------------------------
118  out = (iand(self%bits, bits) ==0)
119 END FUNCTION is_clear
120 
121 !===============================================================================
122 !> Check bits in the status word
123 !===============================================================================
124 FUNCTION mask (self, m) RESULT (out)
125  class(boundaries_t) :: self
126  class(mesh_t), pointer :: m(:)
127  logical :: out(m(1)%gn,m(2)%gn,m(3)%gn)
128  !............................................................................
129  integer :: ix, iy, iz
130  real(8) :: x, y, z, r
131  !-----------------------------------------------------------------------------
132  call trace_begin('boundaries_t%mask')
133  out(:,:,:) = .false.
134  out(m(1)%li:m(1)%ui,m(2)%li:m(2)%ui,m(3)%li:m(3)%ui) = .true.
135  if (self%is_set (bits%spherical)) then
136  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
137  do iz=m(3)%li,m(3)%ui
138  z = m(3)%p + r3(iz) - self%position(3)
139  do iy=m(2)%li,m(2)%ui
140  y = m(2)%p + r2(iy) - self%position(2)
141  do ix=m(1)%li,m(1)%ui
142  x = m(1)%p + r1(ix) - self%position(1)
143  r = sqrt(x**2+y**2+z**2)
144  out(ix,iy,iz) = merge(out(ix,iy,iz), .false., r>self%radius)
145  end do
146  end do
147  end do
148  end associate
149  end if
150  call trace_end
151 END FUNCTION mask
152 
153 !===============================================================================
154 !> Apply a boundary condition, using the procedure pointer "operator"
155 !===============================================================================
156 SUBROUTINE condition (self, f, iv, operator, value, select)
157  class(boundaries_t) :: self
158  real(kind=KindScalarVar) :: f(:,:,:)
159  integer :: iv
160  real, optional :: value
161  procedure(make_symmetric), pointer :: operator
162  integer, optional :: select
163  !.............................................................................
164  class(mesh_t), pointer:: m(:), m1, m2
165  integer:: i, saved_status
166  !-----------------------------------------------------------------------------
167  call trace%begin ('boundaries_t%condition')
168  if (self%unset) then
169  call mpi%abort ('boundaries_t%init has not been called')
170  end if
171  m => self%m
172  if (present(select)) then
173  saved_status = self%bits
174  self%bits = iand(self%bits, select)
175  if (self%is_set(bits%xl+bits%xu)) call x_condition
176  if (self%is_set(bits%yl+bits%yu)) call y_condition
177  if (self%is_set(bits%zl+bits%zu)) call z_condition
178  self%bits = saved_status
179  else
180  if (self%is_set(bits%xl+bits%xu)) call x_condition
181  if (self%is_set(bits%yl+bits%yu)) call y_condition
182  if (self%is_set(bits%zl+bits%zu)) call z_condition
183  end if
184  call trace%end()
185  return
186 contains
187  subroutine x_condition
188  if (self%check_filled) call trace_begin ('x_condition')
189  if (self%id==io%id_debug) &
190  print *, 'MK BOUNDARY CONDITION: x-direction', self%bits, present(value)
191  m1 => self%m(2)
192  m2 => self%m(1)
193  do i=m(3)%lb,m(3)%ub
194  call operator (self, f(:,:,i),m1,m2,iv,lower=self%is_set(bits%xl), &
195  upper=self%is_set(bits%xu), transpose=.true., value=value)
196  end do
197  if (self%check_filled) then
198  if (self%is_set(bits%xl)) then
199  self%filled(self%m(1)%lb:self%m(1)%lo,:,:,iv) = .true.
200  end if
201  if (self%is_set(bits%xu)) then
202  self%filled(self%m(1)%uo:self%m(1)%ub,:,:,iv) = .true.
203  end if
204  end if
205  if (self%check_filled) call trace_end
206  end subroutine x_condition
207  subroutine y_condition
208  if (self%check_filled) call trace_begin ('y_condition')
209  if (self%id==io%id_debug) &
210  print *, 'MK BOUNDARY CONDITION: y-direction', self%bits, present(value)
211  m1 => self%m(1)
212  m2 => self%m(2)
213  do i=m(3)%lb,m(3)%ub
214  call operator (self, f(:,:,i),m1,m2,iv,lower=self%is_set(bits%yl), &
215  upper=self%is_set(bits%yu), transpose=.false., value=value)
216  end do
217  if (self%check_filled) then
218  if (self%is_set(bits%yl)) then
219  self%filled(:,self%m(2)%lb:self%m(2)%lo,:,iv) = .true.
220  end if
221  if (self%is_set(bits%yu)) then
222  self%filled(:,self%m(2)%uo:self%m(2)%ub,:,iv) = .true.
223  end if
224  end if
225  if (self%check_filled) call trace_end
226  end subroutine y_condition
227  subroutine z_condition
228  if (self%check_filled) call trace_begin ('z_condition')
229  if (self%id==io%id_debug) &
230  print *, 'MK BOUNDARY CONDITION: z-direction', self%bits, present(value)
231  m1 => self%m(1)
232  m2 => self%m(3)
233  do i=m(2)%lb,m(2)%ub
234  call operator (self, f(:,i,:),m1,m2,iv,lower=self%is_set(bits%zl), &
235  upper=self%is_set(bits%zu), transpose=.false., value=value)
236  end do
237  if (self%check_filled) then
238  if (self%is_set(bits%zl)) then
239  self%filled(:,:,self%m(3)%lb:self%m(3)%lo,iv) = .true.
240  end if
241  if (self%is_set(bits%zu)) then
242  self%filled(:,:,self%m(3)%uo:self%m(3)%ub,iv) = .true.
243  end if
244  end if
245  if (self%check_filled) call trace_end
246  end subroutine z_condition
247 END SUBROUTINE condition
248 
249 !===============================================================================
250 !> Apply symmetric boundary conditions in the directions indicated by status bits
251 !===============================================================================
252 SUBROUTINE symmetric (self, mem, iv, select)
253  class(boundaries_t) :: self
254  real(kind=KindScalarVar), pointer :: mem(:,:,:)
255  integer :: iv
256  integer, optional :: select
257  procedure(make_symmetric), pointer :: operator
258  !-----------------------------------------------------------------------------
259  operator => make_symmetric
260  call trace%begin ('boundaries_t%symmetric')
261  if (self%id==io%id_debug) &
262  print *, self%id, 'SYMM', iv, self%is_set(bits%yl), self%is_set(bits%yu)
263  call self%condition (mem, iv, operator)
264  call trace%end()
265 END SUBROUTINE symmetric
266 
267 !===============================================================================
268 !> Apply antisymmetric boundary conditions
269 !===============================================================================
270 SUBROUTINE antisymmetric (self, mem, iv, select)
271  class(boundaries_t) :: self
272  real(kind=KindScalarVar), pointer :: mem(:,:,:)
273  integer :: iv
274  integer, optional :: select
275  procedure(make_symmetric), pointer :: operator
276  !-----------------------------------------------------------------------------
277  operator => make_antisymmetric
278  call trace%begin ('boundaries_t%antisymmetric')
279  print *, self%id, 'ANTI', iv, self%is_set(bits%yl), self%is_set(bits%yu)
280  call self%condition (mem, iv, operator)
281  call trace%end()
282 END SUBROUTINE antisymmetric
283 
284 !===============================================================================
285 !> Apply extrapolating boundary conditions
286 !===============================================================================
287 SUBROUTINE extrapolate (self, mem, iv, select)
288  class(boundaries_t) :: self
289  real(kind=KindScalarVar), pointer :: mem(:,:,:)
290  integer :: iv
291  integer, optional :: select
292  procedure(make_symmetric), pointer :: operator
293  !-----------------------------------------------------------------------------
294  operator => make_extrapolated
295  call trace%begin ('boundaries_t%extrapolate')
296  if (self%id==io%id_debug) &
297  print *, self%id, 'EXTRAP', iv, self%is_set(bits%yl), self%is_set(bits%yu)
298  call self%condition (mem, iv, operator)
299  call trace%end()
300 END SUBROUTINE extrapolate
301 
302 !===============================================================================
303 !> Apply extrapolating boundary conditions
304 !===============================================================================
305 SUBROUTINE log_extrapolate (self, mem, iv, select)
306  class(boundaries_t) :: self
307  real(kind=KindScalarVar), pointer :: mem(:,:,:)
308  integer :: iv
309  integer, optional :: select
310  procedure(make_symmetric), pointer :: operator
311  !-----------------------------------------------------------------------------
312  operator => make_log_extrapolated
313  call trace%begin ('boundaries_t%extrapolate')
314  call self%condition (mem, iv, operator)
315  call trace%end()
316 END SUBROUTINE log_extrapolate
317 
318 !===============================================================================
319 !> Apply constant boundary conditions
320 !===============================================================================
321 SUBROUTINE constant (self, mem, iv, value, select)
322  class(boundaries_t) :: self
323  real(kind=KindScalarVar), pointer :: mem(:,:,:)
324  integer :: iv
325  real :: value
326  integer, optional :: select
327  procedure(make_symmetric), pointer :: operator
328  !-----------------------------------------------------------------------------
329  operator => make_constant
330  call trace%begin ('boundaries_t%constant')
331  call self%condition (mem, iv, operator, value=value)
332  call trace%end()
333 END SUBROUTINE constant
334 
335 !===============================================================================
336 !> Set indices that define the location of the physical boundary, by determining
337 !> how symmtry conditions are applied across boundaries. lo and uo are the points
338 !> in the boundary zones closest to the interior, while li and ui are the neartest
339 !> points in the interior, which are being used to construct the boundary zone
340 !> values.
341 !===============================================================================
342 SUBROUTINE bndry_indices (m, iv, lo, li, ui, uo)
343  class(mesh_t), pointer :: m
344  integer :: iv, lo, li, ui, uo
345  !-----------------------------------------------------------------------------
346  ! no-mans-land mesh, with staggered points at patch boundaries
347  ! | |
348  ! centered: x---x-|-o---o-- ... --o---o-|-x---x
349  ! lo li ui uo
350  ! staggered: x---x---o---o-- ... --o---o---x---x
351  ! lo li ui uo
352  !-----------------------------------------------------------------------------
353  if (m%no_mans_land) then
354  ! --- staggered variables ---
355  if (m%h(iv) < 0.0) then
356  lo = m%lo
357  li = lo+2
358  uo = m%uo+1
359  ui = uo-2
360  ! --- centered variables ---
361  else
362  lo = m%lo
363  li = lo+1
364  uo = m%uo
365  ui = uo-1
366  end if
367  !-----------------------------------------------------------------------------
368  ! no-no-mans-land mesh, with centered points at patch boundaries
369  ! | |
370  ! centered: x---x---o---o-- ... --o---o---o---x---x
371  ! lo li ui uo
372  ! staggered: x---x---o-|-o-- ... --o---o---o-|-x---x
373  ! lo li ui uo
374  !-----------------------------------------------------------------------------
375  else
376  ! --- staggered variables ---
377  if (m%h(iv) < 0.0) then
378  lo = m%li
379  li = lo+1
380  uo = m%uo
381  ui = uo-1
382  ! --- centered variables ---
383  else
384  lo = m%lo
385  li = lo+2
386  uo = m%uo
387  ui = uo-2
388  end if
389  end if
390 END SUBROUTINE bndry_indices
391 
392 !===============================================================================
393 !> Apply symmtric boundary condition to the 2nd index, looping over the 1st
394 !===============================================================================
395 SUBROUTINE make_symmetric (self, f, m1, m2, iv, lower, upper, transpose, value)
396  class(boundaries_t) :: self
397  real(kind=KindScalarVar), dimension(:,:) :: f
398  class(mesh_t), pointer :: m1, m2
399  integer :: iv
400  logical :: lower, upper, transpose
401  real, optional :: value
402  integer :: i, j, lo, li, ui, uo
403  !-----------------------------------------------------------------------------
404  call bndry_indices (m2, iv, lo, li, ui, uo)
405  if (lower) then
406  if (transpose) then
407  do j=1,lo
408  do i=m1%lb,m1%ub
409  f(j,i) = f(li+(lo-j),i)
410  end do
411  end do
412  else
413  do j=1,lo
414  do i=m1%lb,m1%ub
415  f(i,j) = f(i,li+(lo-j))
416  end do
417  end do
418  end if
419  end if
420  if (upper) then
421  if (transpose) then
422  do j=uo,m2%ub
423  do i=m1%lb,m1%ub
424  f(j,i) = f(ui-(j-uo),i)
425  end do
426  end do
427  else
428  do j=uo,m2%ub
429  do i=m1%lb,m1%ub
430  f(i,j) = f(i,ui-(j-uo))
431  end do
432  end do
433  end if
434  end if
435 END SUBROUTINE make_symmetric
436 
437 !===============================================================================
438 !> Make the solution antisymmetric across the boundary
439 !===============================================================================
440 SUBROUTINE make_antisymmetric (self, f, m1, m2, iv, lower, upper, transpose, value)
441  class(boundaries_t) :: self
442  real(kind=KindScalarVar), dimension(:,:) :: f
443  class(mesh_t), pointer :: m1, m2
444  integer :: iv
445  logical :: lower, upper, transpose
446  real, optional :: value
447  integer :: i, j, lo, li, ui, uo
448  !-----------------------------------------------------------------------------
449  call bndry_indices (m2, iv, lo, li, ui, uo)
450  if (lower) then
451  if (transpose) then
452  do j=1,lo
453  do i=m1%lb,m1%ub
454  f(j,i) = -f(li+(lo-j),i)
455  end do
456  end do
457  else
458  do j=1,lo
459  do i=m1%lb,m1%ub
460  f(i,j) = -f(i,li+(lo-j))
461  end do
462  end do
463  end if
464  end if
465  if (upper) then
466  if (transpose) then
467  do j=uo,m2%ub
468  do i=m1%lb,m1%ub
469  f(j,i) = -f(ui-(j-uo),i)
470  end do
471  end do
472  else
473  do j=uo,m2%ub
474  do i=m1%lb,m1%ub
475  f(i,j) = -f(i,ui-(j-uo))
476  end do
477  end do
478  end if
479  end if
480 END SUBROUTINE make_antisymmetric
481 
482 !===============================================================================
483 !> Make an extrapolated boundary condition, possibly transposing for vector
484 !> speedup
485 !===============================================================================
486 SUBROUTINE make_extrapolated (self, f, m1, m2, iv, lower, upper, transpose, value)
487  class(boundaries_t) :: self
488  real(kind=KindScalarVar), dimension(:,:) :: f
489  class(mesh_t), pointer :: m1, m2
490  integer :: iv
491  logical :: lower, upper, transpose
492  real, optional :: value
493  integer :: i, j, lo, li, ui, uo
494  !-----------------------------------------------------------------------------
495  call bndry_indices (m2, iv, lo, li, ui, uo)
496  if (lower) then
497  if (transpose) then
498  do j=1,lo
499  do i=m1%lb,m1%ub
500  f(j,i) = 2.*f(lo+1,i) - f(lo+2+(lo-j),i)
501  end do
502  end do
503  else
504  do j=1,lo
505  do i=m1%lb,m1%ub
506  f(i,j) = 2.*f(i,lo+1) - f(i,lo+2+(lo-j))
507  end do
508  end do
509  end if
510  end if
511  if (upper) then
512  if (transpose) then
513  do j=uo,m2%ub
514  do i=m1%lb,m1%ub
515  f(j,i) = 2.*f(uo-1,i) - f(uo-2-(j-uo),i)
516  end do
517  end do
518  else
519  do j=uo,m2%ub
520  do i=m1%lb,m1%ub
521  f(i,j) = 2.*f(i,uo-1) - f(i,uo-2-(j-uo))
522  end do
523  end do
524  end if
525  end if
526 END SUBROUTINE make_extrapolated
527 
528 !===============================================================================
529 !> Make an extrapolated boundary condition, possibly transposing for vector
530 !> speedup
531 !===============================================================================
532 SUBROUTINE make_log_extrapolated (self, f, m1, m2, iv, lower, upper, transpose, value)
533  class(boundaries_t) :: self
534  real(kind=KindScalarVar), dimension(:,:) :: f
535  class(mesh_t), pointer :: m1, m2
536  integer :: iv
537  logical :: lower, upper, transpose
538  real, optional :: value
539  integer :: i, j, lo, li, ui, uo
540  !-----------------------------------------------------------------------------
541  call bndry_indices (m2, iv, lo, li, ui, uo)
542  if (lower) then
543  if (transpose) then
544  do j=1,lo
545  do i=m1%lb,m1%ub
546  f(j,i) = exp(2.*log(f(lo+1,i)) - log(f(lo+2+(lo-j),i)))
547  end do
548  end do
549  else
550  do j=1,lo
551  do i=m1%lb,m1%ub
552  f(i,j) = exp(2.*log(f(i,lo+1)) - log(f(i,lo+2+(lo-j))))
553  end do
554  end do
555  end if
556  end if
557  if (upper) then
558  if (transpose) then
559  do j=uo,m2%ub
560  do i=m1%lb,m1%ub
561  f(j,i) = exp(2.*log(f(uo-1,i)) - log(f(uo-2-(j-uo),i)))
562  end do
563  end do
564  else
565  do j=uo,m2%ub
566  do i=m1%lb,m1%ub
567  f(i,j) = exp(2.*log(f(i,uo-1)) - log(f(i,uo-2-(j-uo))))
568  end do
569  end do
570  end if
571  end if
572 END SUBROUTINE make_log_extrapolated
573 
574 !===============================================================================
575 !> Make an constant boundary condition, possibly transposing for vector
576 !> speedup
577 !===============================================================================
578 SUBROUTINE make_constant (self, f, m1, m2, iv, lower, upper, transpose, value)
579  class(boundaries_t) :: self
580  real(kind=KindScalarVar), dimension(:,:) :: f
581  class(mesh_t), pointer :: m1, m2
582  integer :: iv
583  logical :: lower, upper, transpose
584  real, optional :: value
585  integer :: i, j, lo, li, ui, uo
586  !-----------------------------------------------------------------------------
587  call bndry_indices (m2, iv, lo, li, ui, uo)
588  if (lower) then
589  if (transpose) then
590  do j=1,lo
591  do i=m1%lb,m1%ub
592  f(j,i) = value
593  end do
594  end do
595  else
596  do j=1,lo
597  do i=m1%lb,m1%ub
598  f(i,j) = value
599  end do
600  end do
601  end if
602  end if
603  if (upper) then
604  if (transpose) then
605  do j=uo,m2%ub
606  do i=m1%lb,m1%ub
607  f(j,i) = value
608  end do
609  end do
610  else
611  do j=uo,m2%ub
612  do i=m1%lb,m1%ub
613  f(i,j) = value
614  end do
615  end do
616  end if
617  end if
618 END SUBROUTINE make_constant
619 
620 !===============================================================================
621 !> Compute stagger offsets, taking care with buggy compilers, which may need
622 !> help from specific mesh information.
623 !===============================================================================
624 SUBROUTINE set_stagger (m, h, iv)
625  class(mesh_t), pointer :: m(:)
626  real(8) :: h(3)
627  integer :: iv, i
628  !-----------------------------------------------------------------------------
629  select type (m)
630  type is (cartesian_mesh)
631  do i=1,3
632  h(i) = m(i)%h(iv)
633  end do
634  type is (cylindrical_mesh)
635  do i=1,3
636  h(i) = m(i)%h(iv)
637  end do
638  type is (spherical_mesh)
639  do i=1,3
640  h(i) = m(i)%h(iv)
641  end do
642  class default
643  call mpi%abort ('boundaries_mod::set_stagger: staggering h undefined')
644  end select
645 END SUBROUTINE set_stagger
646 
647 !===============================================================================
648 !> Impose fixed boundary condition; scalar values
649 !===============================================================================
650 SUBROUTINE spherical0 (self, mem, v, m, iv, maxfill)
651  class(boundaries_t) :: self
652  real(kind=KindScalarVar), dimension(:,:,:) :: mem
653  real :: v
654  class(mesh_t), pointer :: m(:)
655  integer :: iv
656  logical, optional :: maxfill
657  !............................................................................
658  integer :: ix, iy, iz
659  real(8) :: x, y, z, r, h(3)
660  !-----------------------------------------------------------------------------
661  call trace_begin('boundaries_t%spherical0')
662  call set_stagger (m, h, iv)
663  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
664  do iz=m(3)%lb,m(3)%ub
665  z = m(3)%p + r3(iz) + h(3)*m(3)%d - self%position(3)
666  do iy=m(2)%lb,m(2)%ub
667  y = m(2)%p + r2(iy) + h(2)*m(2)%d - self%position(2)
668  do ix=m(1)%lb,m(1)%ub
669  x = m(1)%p + r1(ix) + h(1)*m(1)%d - self%position(1)
670  r = sqrt(x**2+y**2+z**2)
671  if (r < self%radius+m(1)%d*0.5) then
672  mem(ix,iy,iz) = v
673  end if
674  end do
675  end do
676  end do
677  end associate
678  call trace_end
679 END SUBROUTINE spherical0
680 
681 !===============================================================================
682 !> Impose vanishing mass flux on any face that is inside the radius; this
683 !> automatically makes density and entropy well behaved there also
684 !===============================================================================
685 SUBROUTINE spherical_p_ns (self, p, m, position, radius)
686  class(boundaries_t) :: self
687  real(kind=KindScalarVar), pointer :: p(:,:,:,:)
688  class(mesh_t), pointer :: m(:)
689  real(8) :: position(3), radius
690  !............................................................................
691  integer :: ix, iy, iz, iv, i
692  real :: h(3), r(3)
693  class(mesh_t), pointer :: mm
694  !-----------------------------------------------------------------------------
695  call trace_begin('boundaries_t%spherical_p_ns')
696  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
697  do iv=1,3
698  do i=1,3
699  mm => m(i)
700  h(i) = mm%h(idx%px+iv-1)*mm%d
701  end do
702  do iz=m(3)%lb,m(3)%ub
703  do iy=m(2)%lb,m(2)%ub
704  do ix=m(1)%lb,m(1)%ub
705  r(1) = m(1)%p - position(1) + r1(ix) + h(1)
706  r(2) = m(2)%p - position(2) + r2(iy) + h(2)
707  r(3) = m(3)%p - position(3) + r3(iz) + h(3)
708  if (sum(r**2) < radius**2) p(ix,iy,iz,iv) = 0.0
709  end do
710  end do
711  end do
712  end do
713  end associate
714  call trace_end
715 END SUBROUTINE spherical_p_ns
716 
717 !===============================================================================
718 !> Impose no-slip boundary condition; scalar values
719 !===============================================================================
720 SUBROUTINE spherical_ns (self, mem, v, m, iv, maxfill)
721  class(boundaries_t) :: self
722  real(kind=KindScalarVar), dimension(:,:,:) :: mem
723  real :: v
724  class(mesh_t), pointer :: m(:)
725  integer :: iv
726  logical, optional :: maxfill
727  !............................................................................
728  integer :: ix, iy, iz
729  real(8) :: x, y, z, r, h(3)
730  real(8) :: d1, dw, x1, rp, p1, p2
731  !-----------------------------------------------------------------------------
732  call trace_begin('boundaries_t%spherical_ns')
733  call set_stagger (m, h, iv)
734  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
735  do iz=m(3)%lb,m(3)%ub
736  z = m(3)%p + r3(iz) + h(3)*m(3)%d - self%position(3)
737  do iy=m(2)%lb,m(2)%ub
738  y = m(2)%p + r2(iy) + h(2)*m(2)%d - self%position(2)
739  do ix=m(1)%lb,m(1)%ub
740  x = m(1)%p + r1(ix) + h(1)*m(1)%d - self%position(1)
741  r = sqrt(x**2+y**2+z**2)
742  if (iv == idx%px) then ! mass flux in X-direction
743  if (r <= self%radius) then
744  p1 = sqrt(y**2+z**2+(m(1)%p + r1(ix-1) + h(1)*m(1)%d - self%position(1))**2)
745  p2 = sqrt(y**2+z**2+(m(1)%p + r1(ix+1) + h(1)*m(1)%d - self%position(1))**2)
746  if ((p1 >= self%radius+m(1)%d*0.5).and.(p2 <= self%radius+m(1)%d*0.5)) then
747  !entering the boundary
748  rp = sqrt(abs(self%radius**2 - y**2 - z**2)) ! boundary location at X axis
749  if (x<0) rp = -rp
750  x1 = m(1)%p + r1(ix-1) + h(1)*m(1)%d - self%position(1) ! last point OUTSIDE the bnd
751  d1 = abs(rp-x1) ! distance between bnd and last point before it
752  dw = d1/(d1+m(1)%d) ! interpolation weight for the last point before the bnd
753  mem(ix-1,iy,iz) = mem(ix-2,iy,iz) * dw ! linear interpolation for the last point before the bnd
754  mem(ix ,iy,iz) = &
755  extrapolate_antisym(x,x1,rp,mem(ix-1,iy,iz)) ! antisymmetric extrapolation from the last point before the bnd
756  else if ((p1 <= self%radius+m(1)%d*0.5).and.(p2 >= self%radius+m(1)%d*0.5)) then
757  ! leaving the boundary
758  rp = sqrt(abs(self%radius**2 - y**2 - z**2)) ! boundary location at X axis
759  if (x<0) rp = -rp
760  x1 = m(1)%p + r1(ix+1) + h(1)*m(1)%d - self%position(1) ! first point OUTSIDE the bnd
761  d1 = abs(rp-x1) ! distance between bnd and the point outside it
762  dw = d1/(d1+m(1)%d) ! interpolation weight for the first point after the bnd
763  mem(ix+1,iy,iz) = mem(ix+2,iy,iz)*dw ! linear interpolation for the first point after the bnd
764  mem(ix ,iy,iz) = &
765  extrapolate_antisym(x,x1,rp,mem(ix+1,iy,iz)) ! antisymmetric extrapolation from the first point after the bnd
766  else if ((p1 <= self%radius+m(1)%d*0.5).and.(p2 <= self%radius+m(1)%d*0.5)) then
767  mem(ix,iy,iz) = 0.0 ! cell is well inside the boundary
768  else
769  rp = sqrt(abs(self%radius**2 - y**2 - z**2)) ! boundary location at X axis
770  if (x<0) rp = -rp
771  x1 = m(1)%p + r1(ix-1) + h(1)*m(1)%d - self%position(1) ! last point INSIDE the bnd
772  d1 = abs(rp-x1) ! distance between bnd and last point before it
773  dw = d1/(d1+m(1)%d) ! interpolation weight for the last point before the bnd
774  mem(ix-1,iy,iz) = mem(ix-2,iy,iz) * dw ! linear interpolation for the last point before the bnd
775  x1 = m(1)%p + r1(ix+1) + h(1)*m(1)%d - self%position(1) ! first point OUTSIDE the bnd
776  d1 = abs(rp-x1) ! distance between bnd and last point before it
777  dw = d1/(d1+m(1)%d) ! interpolation weight for the last point before the bnd
778  mem(ix+1,iy,iz) = mem(ix+2,iy,iz) * dw ! linear interpolation for the last point before the bnd
779  mem(ix,iy,iz) = 0.0 ! zero for the bnd cell
780  end if
781  end if
782  else if (iv == idx%py) then ! mass flux in Y-direction
783  if (r <= self%radius) then
784  p1 = sqrt(x**2+(m(2)%p + r2(iy-1) + h(2)*m(2)%d - self%position(2))**2+z**2)
785  p2 = sqrt(x**2+(m(2)%p + r2(iy+1) + h(2)*m(2)%d - self%position(2))**2+z**2)
786  if ((p1 > self%radius+m(2)%d*0.5).and.(p2 < self%radius+m(2)%d*0.5)) then
787  !entering the boundary
788  rp = sqrt(abs(self%radius**2 - x**2 - z**2)) ! boundary location at Y axis
789  if (y<0) rp = -rp
790  x1 = m(2)%p + r2(iy-1) + h(2)*m(2)%d - self%position(2) ! last point OUTSIDE the bnd
791  d1 = abs(rp-x1) ! distance between bnd and last point before it
792  dw = d1/(d1+m(2)%d) ! interpolation weight for the last point before the bnd
793  mem(ix,iy-1,iz) = mem(ix,iy-2,iz) * dw ! linear interpolation for the last point before the bnd
794  mem(ix,iy ,iz) = &
795  extrapolate_antisym(y,x1,rp,mem(ix,iy-1,iz)) ! antisymmetric extrapolation from the last point before the bnd
796  else if ((p1 < self%radius+m(2)%d*0.5).and.(p2 > self%radius+m(2)%d*0.5)) then
797  ! leaving the boundary
798  rp = sqrt(abs(self%radius**2 - x**2 - z**2)) ! boundary location at Y axis
799  if (y<0) rp = -rp
800  x1 = m(2)%p + r2(iy+1) + h(2)*m(2)%d - self%position(2) ! first point OUTSIDE the bnd
801  d1 = abs(rp-x1) ! distance between bnd and the point outside it
802  dw = d1/(d1+m(2)%d) ! interpolation weight for the first point after the bnd
803  mem(ix,iy+1,iz) = mem(ix,iy+2,iz)*dw ! linear interpolation for the first point after the bnd
804  mem(ix,iy,iz) = &
805  extrapolate_antisym(y,x1,rp,mem(ix,iy+1,iz)) ! antisymmetric extrapolation from the first point after the bnd
806  else if ((p1 < self%radius+m(2)%d*0.5).and.(p2 < self%radius+m(2)%d*0.5)) then
807  mem(ix,iy,iz) = 0.0 ! cell is well inside the boundary
808  else ! special case, where a single cell gets a bc
809  rp = sqrt(abs(self%radius**2 - x**2 - z**2)) ! boundary location at Y axis
810  if (y<0) rp = -rp
811  x1 = m(2)%p + r2(iy-1) + h(2)*m(2)%d - self%position(2) ! last point OUTSIDE the bnd
812  d1 = abs(rp-x1) ! distance between bnd and last point before it
813  dw = d1/(d1+m(2)%d) ! interpolation weight for the last point before the bnd
814  mem(ix,iy-1,iz) = mem(ix,iy-2,iz) * dw ! linear interpolation for the last point before the bnd
815  x1 = m(2)%p + r2(iy+1) + h(2)*m(2)%d - self%position(2) ! first point OUTSIDE the bnd
816  d1 = abs(rp-x1) ! distance between bnd and last point before it
817  dw = d1/(d1+m(2)%d) ! interpolation weight for the last point before the bnd
818  mem(ix,iy+1,iz) = mem(ix,iy+2,iz) * dw ! linear interpolation for the last point before the bnd
819  mem(ix,iy,iz) = 0.0 ! zero for the bnd cell
820  end if
821  end if
822  else if (iv == idx%pz) then ! mass flux in Z-direction
823  if (r <= self%radius) then
824  p1 = sqrt(x**2+y**2+(m(3)%p + r3(iz-1) + h(3)*m(3)%d - self%position(3))**2)
825  p2 = sqrt(x**2+y**2+(m(3)%p + r3(iz+1) + h(3)*m(3)%d - self%position(3))**2)
826  if ((p1 > self%radius+m(3)%d*0.5).and.(p2 < self%radius+m(3)%d*0.5)) then
827  !entering the boundary
828  rp = sqrt(abs(self%radius**2 - x**2 - y**2)) ! boundary location at Z axis
829  if (z<0) rp = -rp
830  x1 = m(3)%p + r3(iz-1) + h(3)*m(3)%d - self%position(3) ! last point OUTSIDE the bnd
831  d1 = abs(rp-x1) ! distance between bnd and last point before it
832  dw = d1/(d1+m(3)%d) ! interpolation weight for the last point before the bnd
833  mem(ix,iy,iz-1) = mem(ix,iy,iz-2) * dw ! linear interpolation for the last point before the bnd
834  mem(ix,iy,iz ) = &
835  extrapolate_antisym(z,x1,rp,mem(ix,iy,iz-1)) ! antisymmetric extrapolation from the last point before the bnd
836  else if ((p1 < self%radius+m(3)%d*0.5).and.(p2 > self%radius+m(3)%d*0.5)) then
837  ! leaving the boundary
838  rp = sqrt(abs(self%radius**2 - x**2 - y**2)) ! boundary location at Z axis
839  if (y<0) rp = -rp
840  x1 = m(3)%p + r3(iz+1) + h(3)*m(3)%d - self%position(3) ! last point INSIDE the bnd
841  d1 = abs(rp-x1) ! distance between bnd and the point outside it
842  dw = d1/(d1+m(3)%d) ! interpolation weight for the first point after the bnd
843  mem(ix,iy,iz+1) = mem(ix,iy,iz+2)*dw ! linear interpolation for the first point after the bnd
844  mem(ix,iy,iz) = &
845  extrapolate_antisym(z,x1,rp,mem(ix,iy,iz+1)) ! antisymmetric extrapolation from the first point after the bnd
846  else if ((p1 < self%radius+m(3)%d*0.5).and.(p2 < self%radius+m(3)%d*0.5)) then
847  mem(ix,iy,iz) = 0.0 ! cell is well inside the boundary
848  else
849  rp = sqrt(abs(self%radius**2 - x**2 - y**2)) ! boundary location at Z axis
850  if (z<0) rp = -rp
851  x1 = m(3)%p + r3(iz-1) + h(3)*m(3)%d - self%position(3) ! last point OUTSIDE the bnd
852  d1 = abs(rp-x1) ! distance between bnd and last point before it
853  dw = d1/(d1+m(3)%d) ! interpolation weight for the last point before the bnd
854  mem(ix,iy,iz-1) = mem(ix,iy,iz-2) * dw ! linear interpolation for the last point before the bnd
855  x1 = m(3)%p + r3(iz+1) + h(3)*m(3)%d - self%position(3) ! first point OUTSIDE the bnd
856  d1 = abs(rp-x1) ! distance between bnd and last point before it
857  dw = d1/(d1+m(3)%d) ! interpolation weight for the last point before the bnd
858  mem(ix,iy,iz+1) = mem(ix,iy,iz+2) * dw ! linear interpolation for the last point before the bnd
859  mem(ix,iy,iz) = 0.0 ! zero for the bnd cell
860  end if
861  end if
862  end if
863  end do
864  end do
865  end do
866  ! error stop "end"
867  end associate
868  call trace_end
869 END SUBROUTINE spherical_ns
870 
871 !===============================================================================
872 !> Impose spherical density boundary condition, assuming that the sphere is jagged.
873 !===============================================================================
874 SUBROUTINE spherical_jagged_d (self, mem, m, iv, it, new)
875  class(boundaries_t) :: self
876  real(kind=KindScalarVar), dimension(:,:,:,:) :: mem
877  class(mesh_t), pointer :: m(:), m1, m2, m3
878  integer :: iv, it, new
879  !............................................................................
880  integer :: ix, iy, iz
881  real(8) :: x1, x2, y1, y2, z1, z2
882  real(8) :: r(3)
883  logical :: inside(8)
884  !-----------------------------------------------------------------------------
885  call trace_begin('boundaries_t%spherical_jagged_d')
886  m1 => m(1)
887  m2 => m(2)
888  m3 => m(3)
889  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
890  do iz=m(3)%lb,m(3)%ub
891  do iy=m(2)%lb,m(2)%ub
892  do ix=m(1)%lb,m(1)%ub
893  r(1) = m(1)%p - self%position(1) + r1(ix)
894  r(2) = m(2)%p - self%position(2) + r2(iy)
895  r(3) = m(3)%p - self%position(3) + r3(iz)
896  if (sum(r**2) < self%radius**2) mem(ix,iy,iz,new) = &
897  mem(ix,iy,iz,it)
898  end do
899  end do
900  end do
901  end associate
902  nullify (m1, m2, m3)
903  call trace_end
904 END SUBROUTINE spherical_jagged_d
905 
906 !===============================================================================
907 !> Impose spherical density boundary condition, assuming that the sphere is jagged.
908 !===============================================================================
909 SUBROUTINE spherical_jagged_s (self, new, old, m)
910  class(boundaries_t) :: self
911  real(kind=KindScalarVar), dimension(:,:,:) :: old, new
912  class(mesh_t), pointer :: m(:), m1, m2, m3
913  !............................................................................
914  type(index_t):: idx
915  integer :: ix, iy, iz
916  real(8) :: r(3)
917  logical :: inside(8)
918  !-----------------------------------------------------------------------------
919  call trace_begin('boundaries_t%spherical_jagged_s')
920  m1 => m(1)
921  m2 => m(2)
922  m3 => m(3)
923  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
924  do iz=m(3)%lb,m(3)%ub
925  do iy=m(2)%lb,m(2)%ub
926  do ix=m(1)%lb,m(1)%ub
927  r(1) = m(1)%p - self%position(1) + r1(ix)
928  r(2) = m(2)%p - self%position(2) + r2(iy)
929  r(3) = m(3)%p - self%position(3) + r3(iz)
930  if (sum(r**2) < self%radius**2) new(ix,iy,iz) = &
931  old(ix,iy,iz)
932  end do
933  end do
934  end do
935  end associate
936  nullify (m1, m2, m3)
937  call trace_end
938 END SUBROUTINE spherical_jagged_s
939 
940 !> =============================================================================
941 !> Impose spherical density boundary condition, assuming that the sphere is jagged.
942 !> =============================================================================
943 SUBROUTINE spherical_jagged_s0 (self, mem, m, iv, v)
944  class(boundaries_t) :: self
945  real(kind=KindScalarVar), dimension(:,:,:) :: mem
946  class(mesh_t), pointer :: m(:), m1, m2, m3
947  integer :: iv, it
948  real :: v
949  !............................................................................
950  integer :: ix, iy, iz
951  real(8) :: x1, x2, y1, y2, z1, z2
952  real(8) :: r2
953  logical :: inside(8)
954  !-----------------------------------------------------------------------------
955  call trace_begin('boundaries_t%spherical_jagged_s0')
956  m1 => m(1)
957  m2 => m(2)
958  m3 => m(3)
959  r2 = self%radius**2
960  if (iv == idx%s) then
961  do iz = m3%lb, m3%ub
962  do iy = m2%lb, m2%ub
963  do ix = m1%lb, m1%ub
964  x1 = m1%p - self%position(1) + m1%r(ix) + 0.5*m1%d
965  x2 = m1%p - self%position(1) + m1%r(ix) - 0.5*m1%d
966  y1 = m2%p - self%position(2) + m2%r(iy) + 0.5*m2%d
967  y2 = m2%p - self%position(2) + m2%r(iy) - 0.5*m2%d
968  z1 = m3%p - self%position(3) + m3%r(iz) + 0.5*m3%d
969  z2 = m3%p - self%position(3) + m3%r(iz) - 0.5*m3%d
970  inside(1) = x1**2+y1**2+z1**2 < r2
971  inside(2) = x2**2+y1**2+z1**2 < r2
972  inside(3) = x1**2+y2**2+z1**2 < r2
973  inside(4) = x2**2+y2**2+z1**2 < r2
974  inside(5) = x1**2+y1**2+z2**2 < r2
975  inside(6) = x2**2+y1**2+z2**2 < r2
976  inside(7) = x1**2+y2**2+z2**2 < r2
977  inside(8) = x2**2+y2**2+z2**2 < r2
978  if (all(inside)) mem(ix,iy,iz) = v
979  end do
980  end do
981  end do
982  else
983  print *, 'iv passed',iv, 'iv required', idx%s
984  call mpi%abort ("boundaries_mod::spherical_jagged_d:: invalid memory index iv")
985  end if
986  nullify (m1, m2, m3)
987  call trace_end
988 END SUBROUTINE spherical_jagged_s0
989 
990 !===============================================================================
991 !> Same as above, but for momentum
992 !===============================================================================
993 SUBROUTINE spherical_jagged_p (self, mem, m, iv)
994  class(boundaries_t) :: self
995  real(kind=KindScalarVar), dimension(:,:,:) :: mem
996  class(mesh_t), pointer :: m(:), m1, m2, m3
997  integer :: iv
998  !............................................................................
999  integer :: ix, iy, iz
1000  real(8) :: x1, x2, y1, y2, z1, z2
1001  real(8) :: r2, h(3)
1002  logical :: inside(4)
1003  !-----------------------------------------------------------------------------
1004  call trace_begin('boundaries_t%spherical_jagged_p')
1005  m1 => m(1)
1006  m2 => m(2)
1007  m3 => m(3)
1008  call set_stagger (m, h, iv)
1009  r2 = self%radius**2
1010  if (iv == idx%px) then
1011  do iz = m3%lb, m3%ub
1012  do iy = m2%lb, m2%ub
1013  do ix = m1%lb, m1%ub
1014  y1 = m2%p - self%position(2) + m2%r(iy) + 0.5*m2%d
1015  y2 = m2%p - self%position(2) + m2%r(iy) - 0.5*m2%d
1016  z1 = m3%p - self%position(3) + m3%r(iz) + 0.5*m3%d
1017  z2 = m3%p - self%position(3) + m3%r(iz) - 0.5*m3%d
1018  x1 = m1%p - self%position(1) + m1%r(ix) + h(1)*m1%d
1019  inside(1) = y1**2+z1**2 + x1**2 < r2
1020  inside(2) = y2**2+z1**2 + x1**2 < r2
1021  inside(3) = y1**2+z2**2 + x1**2 < r2
1022  inside(4) = y2**2+z2**2 + x1**2 < r2
1023  if (any(inside)) mem(ix,iy,iz) = 0.0
1024  end do
1025  end do
1026  end do
1027  else if (iv == idx%py) then
1028  do iz = m3%lb, m3%ub
1029  do iy = m2%lb, m2%ub
1030  do ix = m1%lb, m1%ub
1031  x1 = m1%p - self%position(1) + m1%r(ix) + 0.5*m1%d
1032  x2 = m1%p - self%position(1) + m1%r(ix) - 0.5*m1%d
1033  z1 = m3%p - self%position(3) + m3%r(iz) + 0.5*m3%d
1034  z2 = m3%p - self%position(3) + m3%r(iz) - 0.5*m3%d
1035  y1 = m2%p - self%position(2) + m2%r(iy) + h(2)*m2%d
1036  inside(1) = x1**2+z1**2 + y1**2 < r2
1037  inside(2) = x2**2+z1**2 + y1**2 < r2
1038  inside(3) = x1**2+z2**2 + y1**2 < r2
1039  inside(4) = x2**2+z2**2 + y1**2 < r2
1040  if (any(inside)) mem(ix,iy,iz) = 0.0
1041  end do
1042  end do
1043  end do
1044  else if (iv == idx%pz) then
1045  do iz = m3%lb, m3%ub
1046  do iy = m2%lb, m2%ub
1047  do ix = m1%lb, m1%ub
1048  x1 = m1%p - self%position(1) + m1%r(ix) + 0.5*m1%d
1049  x2 = m1%p - self%position(1) + m1%r(ix) - 0.5*m1%d
1050  y1 = m2%p - self%position(2) + m2%r(iy) + 0.5*m2%d
1051  y2 = m2%p - self%position(2) + m2%r(iy) - 0.5*m2%d
1052  z1 = m3%p - self%position(3) + m3%r(iz) + h(3)*m3%d
1053  inside(1) = x1**2+y1**2 + z1**2 < r2
1054  inside(2) = x2**2+y1**2 + z1**2 < r2
1055  inside(3) = x1**2+y2**2 + z1**2 < r2
1056  inside(4) = x2**2+y2**2 + z1**2 < r2
1057  if (any(inside)) mem(ix,iy,iz) = 0.0
1058  end do
1059  end do
1060  end do
1061  else
1062  print *, 'iv passed',iv, 'iv required', idx%px, idx%py, idx%pz
1063  call mpi%abort ("boundaries_mod::spherical_jagged_p:: invalid memory index iv")
1064  end if
1065  nullify (m1, m2, m3)
1066  call trace_end
1067 END SUBROUTINE spherical_jagged_p
1068 
1069 !=======================================================================
1070 !> Linear antisymmetric extrapolation
1071 !=======================================================================
1072 FUNCTION extrapolate_antisym(ip, op, bp, yop) result(out)
1073  real(8) :: ip, op, bp
1074  real(kind=KindScalarVar):: yop, out
1075  !.....................................................................
1076  out = yop+(ip-op)/(bp-op)*(-yop)
1077 END FUNCTION extrapolate_antisym
1078 
1079 !===============================================================================
1080 !> Impose fixed boundary condition; array values
1081 !===============================================================================
1082 SUBROUTINE spherical3 (self, mem, v, m, iv, maxfill)
1083  class(boundaries_t) :: self
1084  real(kind=KindScalarVar), dimension(:,:,:) :: mem, v
1085  class(mesh_t), pointer :: m(:)
1086  integer :: iv
1087  logical, optional :: maxfill
1088  !............................................................................
1089  integer :: ix, iy, iz, n1, n2
1090  real(8) :: x, y, z, r, h(3)
1091  !-----------------------------------------------------------------------------
1092  call trace_begin('boundaries_t%spherical3')
1093  call set_stagger (m, h, iv)
1094  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
1095  self%vmax = 0.0
1096  do iz=m(3)%lb,m(3)%ub
1097  z = m(3)%p + r3(iz) + h(3)*m(3)%d - self%position(3)
1098  do iy=m(2)%lb,m(2)%ub
1099  y = m(2)%p + r2(iy) + h(2)*m(2)%d - self%position(2)
1100  do ix=m(1)%lb,m(1)%ub
1101  x = m(1)%p + r1(ix) + h(1)*m(1)%d - self%position(1)
1102  r = sqrt(x**2+y**2+z**2)
1103  if (r < self%radius+m(1)%d*0.5) then
1104  mem(ix,iy,iz) = v(ix,iy,iz)
1105  else if (present(maxfill)) then
1106  self%vmax = max(self%vmax,mem(ix,iy,iz))
1107  end if
1108  end do
1109  end do
1110  end do
1111  if (present(maxfill)) then
1112  do iz=m(3)%lb,m(3)%ub
1113  z = m(3)%p + r3(iz) + h(3)*m(3)%d - self%position(3)
1114  do iy=m(2)%lb,m(2)%ub
1115  y = m(2)%p + r2(iy) + h(2)*m(2)%d - self%position(2)
1116  do ix=m(1)%lb,m(1)%ub
1117  x = m(1)%p + r1(ix) + h(1)*m(1)%d - self%position(1)
1118  r = sqrt(x**2+y**2+z**2)
1119  if (r < self%radius-m(1)%d*0.1) then
1120  mem(ix,iy,iz) = max(mem(ix,iy,iz),self%vmax)
1121  end if
1122  end do
1123  end do
1124  end do
1125  end if
1126  end associate
1127  call trace_end
1128 END SUBROUTINE spherical3
1129 
1130 !===============================================================================
1131 !> Extrapolate values across a spherical boundary, sampling values 4 mesh points
1132 !> from the boundary, and at the boundary, to establish a stable slope
1133 !===============================================================================
1134 SUBROUTINE spherical_extrapolate (self, mem, m, iv)
1135  class(boundaries_t) :: self
1136  real(kind=KindScalarVar), dimension(:,:,:) :: mem
1137  class(mesh_t), pointer :: m(:)
1138  integer :: iv
1139  !............................................................................
1140  integer :: ix, iy, iz, n1, n2
1141  real(8) :: x, y, z, r, h(3), a1, a2, ds, p
1142  !-----------------------------------------------------------------------------
1143  call trace_begin('boundaries_t%spherical_extrapolate')
1144  call set_stagger (m, h, iv)
1145  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
1146  ds = m(1)%d
1147  a1 = 0d0; n1 = 0
1148  a2 = 0d0; n2 = 0
1149  !-----------------------------------------------------------------------------
1150  ! a1 and a2 sample the values in annuli around R-3*ds and R, respectively
1151  !-----------------------------------------------------------------------------
1152  do iz=m(3)%lb,m(3)%ub
1153  z = m(3)%p + r3(iz) + h(3)*m(3)%d - self%position(3)
1154  do iy=m(2)%lb,m(2)%ub
1155  y = m(2)%p + r2(iy) + h(2)*m(2)%d - self%position(2)
1156  do ix=m(1)%lb,m(1)%ub
1157  x = m(1)%p + r1(ix) + h(1)*m(1)%d - self%position(1)
1158  r = sqrt(x**2+y**2+z**2)
1159  if (abs(r-(self%radius+2.*ds)) <= ds) then
1160  n1 = n1+1
1161  a1 = a1 + mem(ix,iy,iz)
1162  else if (abs(r-self%radius) <= ds) then
1163  n2 = n2+1
1164  a2 = a2 + mem(ix,iy,iz)
1165  end if
1166  end do
1167  end do
1168  end do
1169  if (n1>0) a1 = log(a1/n1)
1170  if (n2>0) a2 = log(a2/n2)
1171  a2 = max(a2,a1)
1172  !-----------------------------------------------------------------------------
1173  ! For radii inside R-0.5*ds, set the value to be log-interpolated from a1, a2
1174  ! The value at a2 would be chosen for r=R
1175  !-----------------------------------------------------------------------------
1176  do iz=m(3)%lb,m(3)%ub
1177  z = m(3)%p + r3(iz) + h(3)*m(3)%d - self%position(3)
1178  do iy=m(2)%lb,m(2)%ub
1179  y = m(2)%p + r2(iy) + h(2)*m(2)%d - self%position(2)
1180  do ix=m(1)%lb,m(1)%ub
1181  x = m(1)%p + r1(ix) + h(1)*m(1)%d - self%position(1)
1182  r = sqrt(x**2+y**2+z**2)
1183  r = max(r,self%radius-4.0*ds)
1184  if (r < self%radius) then
1185  if (a1 /= 0d0 .and. a2 /= 0d0) then
1186  p = (self%radius-r)/(2.0*ds) + 1.0
1187  mem(ix,iy,iz) = exp(a1*(1d0-p) + a2*p)
1188  end if
1189  end if
1190  end do
1191  end do
1192  end do
1193  end associate
1194  call trace_end
1195 END SUBROUTINE spherical_extrapolate
1196 
1197 !===============================================================================
1198 !> Set constant entropy inside the boundary radius (included)
1199 !===============================================================================
1200 SUBROUTINE spherical_constant_entropy (self, d, e, m, s, gamma)
1201  class(boundaries_t) :: self
1202  real(kind=KindScalarVar), dimension(:,:,:) :: d, e
1203  class(mesh_t), pointer :: m(:)
1204  real :: s, gamma
1205  !............................................................................
1206  integer :: ix, iy, iz
1207  real(8) :: x, y, z, r, ds, c
1208  !-----------------------------------------------------------------------------
1209  call trace_begin('boundaries_t%sperical_entropy')
1210  !-----------------------------------------------------------------------------
1211  ! a1 and a2 sample the values in annuli around R-3*ds and R, respectively
1212  !-----------------------------------------------------------------------------
1213  c = exp(s*(gamma-1.0))/(gamma-1.0)
1214  ds = m(1)%d
1215  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
1216  do iz=m(3)%lb,m(3)%ub
1217  z = m(3)%p + r3(iz) - self%position(3)
1218  do iy=m(2)%lb,m(2)%ub
1219  y = m(2)%p + r2(iy) - self%position(2)
1220  do ix=m(1)%lb,m(1)%ub
1221  x = m(1)%p + r1(ix)- self%position(1)
1222  r = sqrt(x**2+y**2+z**2)
1223  if (r < self%radius) then
1224  e(ix,iy,iz) = c*d(ix,iy,iz)**gamma
1225  if (iy==18.and.iz==18) print *,'BC: s', m(1)%id, ix, s
1226  end if
1227  end do
1228  end do
1229  end do
1230  end associate
1231  call trace_end
1232 END SUBROUTINE spherical_constant_entropy
1233 
1234 END MODULE boundaries_mod
Boundary conditions for centered variables, which have the physical boundary in-between mesh points l...
Template module for mesh.
Definition: mesh_mod.f90:4
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
Definition: index_mod.f90:7
Definition: io_mod.f90:4