13 USE boundary_wavekill_mod
19 logical:: unset = .true.
20 class(
mesh_t),
pointer:: m(:)
21 real(8):: position(3)=0d0, radius=0d0
23 logical(kind=2),
pointer:: filled(:,:,:,:)
24 logical:: check_filled=.false.
33 procedure:: antisymmetric
34 procedure:: extrapolate
35 procedure:: log_extrapolate
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
56 SUBROUTINE init (self, mesh, position, radius)
58 class(
mesh_t),
pointer,
optional :: mesh(:)
59 real(8),
optional :: position(3), radius
61 real(8) :: size(3), p(3)
63 logical :: one_inside, all_inside
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 70 if (
present(position).and.
present(radius))
then 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 81 end do; end do; end do
82 if (one_inside)
call self%set (bits%spherical)
92 SUBROUTINE set (self, bits)
96 self%bits = ior(self%bits, bits)
102 FUNCTION is_set (self, bits)
RESULT (out)
107 out = (iand(self%bits, bits) /=0)
113 FUNCTION is_clear (self, bits)
RESULT (out)
118 out = (iand(self%bits, bits) ==0)
119 END FUNCTION is_clear
124 FUNCTION mask (self, m)
RESULT (out)
126 class(mesh_t),
pointer :: m(:)
127 logical :: out(m(1)%gn,m(2)%gn,m(3)%gn)
129 integer :: ix, iy, iz
130 real(8) :: x, y, z, r
132 call trace_begin(
'boundaries_t%mask')
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)
156 SUBROUTINE condition (self, f, iv, operator, value, select)
158 real(kind=KindScalarVar) :: f(:,:,:)
160 real,
optional :: value
161 procedure(make_symmetric),
pointer :: operator
162 integer,
optional :: select
164 class(mesh_t),
pointer:: m(:), m1, m2
165 integer:: i, saved_status
167 call trace%begin (
'boundaries_t%condition')
169 call mpi%abort (
'boundaries_t%init has not been called')
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
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
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)
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)
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.
201 if (self%is_set(bits%xu))
then 202 self%filled(self%m(1)%uo:self%m(1)%ub,:,:,iv) = .true.
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)
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)
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.
221 if (self%is_set(bits%yu))
then 222 self%filled(:,self%m(2)%uo:self%m(2)%ub,:,iv) = .true.
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)
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)
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.
241 if (self%is_set(bits%zu))
then 242 self%filled(:,:,self%m(3)%uo:self%m(3)%ub,iv) = .true.
245 if (self%check_filled)
call trace_end
246 end subroutine z_condition
247 END SUBROUTINE condition
252 SUBROUTINE symmetric (self, mem, iv, select)
254 real(kind=KindScalarVar),
pointer :: mem(:,:,:)
256 integer,
optional :: select
257 procedure(make_symmetric),
pointer :: operator
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)
265 END SUBROUTINE symmetric
270 SUBROUTINE antisymmetric (self, mem, iv, select)
272 real(kind=KindScalarVar),
pointer :: mem(:,:,:)
274 integer,
optional :: select
275 procedure(make_symmetric),
pointer :: operator
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)
282 END SUBROUTINE antisymmetric
287 SUBROUTINE extrapolate (self, mem, iv, select)
289 real(kind=KindScalarVar),
pointer :: mem(:,:,:)
291 integer,
optional :: select
292 procedure(make_symmetric),
pointer :: operator
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)
300 END SUBROUTINE extrapolate
305 SUBROUTINE log_extrapolate (self, mem, iv, select)
307 real(kind=KindScalarVar),
pointer :: mem(:,:,:)
309 integer,
optional :: select
310 procedure(make_symmetric),
pointer :: operator
312 operator => make_log_extrapolated
313 call trace%begin (
'boundaries_t%extrapolate')
314 call self%condition (mem, iv, operator)
316 END SUBROUTINE log_extrapolate
321 SUBROUTINE constant (self, mem, iv, value, select)
323 real(kind=KindScalarVar),
pointer :: mem(:,:,:)
326 integer,
optional :: select
327 procedure(make_symmetric),
pointer :: operator
329 operator => make_constant
330 call trace%begin (
'boundaries_t%constant')
331 call self%condition (mem, iv, operator,
value=
value)
333 END SUBROUTINE constant
342 SUBROUTINE bndry_indices (m, iv, lo, li, ui, uo)
343 class(mesh_t),
pointer :: m
344 integer :: iv, lo, li, ui, uo
353 if (m%no_mans_land)
then 355 if (m%h(iv) < 0.0)
then 377 if (m%h(iv) < 0.0)
then 390 END SUBROUTINE bndry_indices
395 SUBROUTINE make_symmetric (self, f, m1, m2, iv, lower, upper, transpose, value)
397 real(kind=KindScalarVar),
dimension(:,:) :: f
398 class(mesh_t),
pointer :: m1, m2
400 logical :: lower, upper, transpose
401 real,
optional :: value
402 integer :: i, j, lo, li, ui, uo
404 call bndry_indices (m2, iv, lo, li, ui, uo)
409 f(j,i) = f(li+(lo-j),i)
415 f(i,j) = f(i,li+(lo-j))
424 f(j,i) = f(ui-(j-uo),i)
430 f(i,j) = f(i,ui-(j-uo))
435 END SUBROUTINE make_symmetric
440 SUBROUTINE make_antisymmetric (self, f, m1, m2, iv, lower, upper, transpose, value)
442 real(kind=KindScalarVar),
dimension(:,:) :: f
443 class(mesh_t),
pointer :: m1, m2
445 logical :: lower, upper, transpose
446 real,
optional :: value
447 integer :: i, j, lo, li, ui, uo
449 call bndry_indices (m2, iv, lo, li, ui, uo)
454 f(j,i) = -f(li+(lo-j),i)
460 f(i,j) = -f(i,li+(lo-j))
469 f(j,i) = -f(ui-(j-uo),i)
475 f(i,j) = -f(i,ui-(j-uo))
480 END SUBROUTINE make_antisymmetric
486 SUBROUTINE make_extrapolated (self, f, m1, m2, iv, lower, upper, transpose, value)
488 real(kind=KindScalarVar),
dimension(:,:) :: f
489 class(mesh_t),
pointer :: m1, m2
491 logical :: lower, upper, transpose
492 real,
optional :: value
493 integer :: i, j, lo, li, ui, uo
495 call bndry_indices (m2, iv, lo, li, ui, uo)
500 f(j,i) = 2.*f(lo+1,i) - f(lo+2+(lo-j),i)
506 f(i,j) = 2.*f(i,lo+1) - f(i,lo+2+(lo-j))
515 f(j,i) = 2.*f(uo-1,i) - f(uo-2-(j-uo),i)
521 f(i,j) = 2.*f(i,uo-1) - f(i,uo-2-(j-uo))
526 END SUBROUTINE make_extrapolated
532 SUBROUTINE make_log_extrapolated (self, f, m1, m2, iv, lower, upper, transpose, value)
534 real(kind=KindScalarVar),
dimension(:,:) :: f
535 class(mesh_t),
pointer :: m1, m2
537 logical :: lower, upper, transpose
538 real,
optional :: value
539 integer :: i, j, lo, li, ui, uo
541 call bndry_indices (m2, iv, lo, li, ui, uo)
546 f(j,i) = exp(2.*log(f(lo+1,i)) - log(f(lo+2+(lo-j),i)))
552 f(i,j) = exp(2.*log(f(i,lo+1)) - log(f(i,lo+2+(lo-j))))
561 f(j,i) = exp(2.*log(f(uo-1,i)) - log(f(uo-2-(j-uo),i)))
567 f(i,j) = exp(2.*log(f(i,uo-1)) - log(f(i,uo-2-(j-uo))))
572 END SUBROUTINE make_log_extrapolated
578 SUBROUTINE make_constant (self, f, m1, m2, iv, lower, upper, transpose, value)
580 real(kind=KindScalarVar),
dimension(:,:) :: f
581 class(mesh_t),
pointer :: m1, m2
583 logical :: lower, upper, transpose
584 real,
optional :: value
585 integer :: i, j, lo, li, ui, uo
587 call bndry_indices (m2, iv, lo, li, ui, uo)
618 END SUBROUTINE make_constant
624 SUBROUTINE set_stagger (m, h, iv)
625 class(mesh_t),
pointer :: m(:)
630 type is (cartesian_mesh)
634 type is (cylindrical_mesh)
638 type is (spherical_mesh)
643 call mpi%abort (
'boundaries_mod::set_stagger: staggering h undefined')
645 END SUBROUTINE set_stagger
650 SUBROUTINE spherical0 (self, mem, v, m, iv, maxfill)
652 real(kind=KindScalarVar),
dimension(:,:,:) :: mem
654 class(mesh_t),
pointer :: m(:)
656 logical,
optional :: maxfill
658 integer :: ix, iy, iz
659 real(8) :: x, y, z, r, h(3)
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 679 END SUBROUTINE spherical0
685 SUBROUTINE spherical_p_ns (self, p, m, position, radius)
687 real(kind=KindScalarVar),
pointer :: p(:,:,:,:)
688 class(mesh_t),
pointer :: m(:)
689 real(8) :: position(3), radius
691 integer :: ix, iy, iz, iv, i
693 class(mesh_t),
pointer :: mm
695 call trace_begin(
'boundaries_t%spherical_p_ns')
696 associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
700 h(i) = mm%h(idx%px+iv-1)*mm%d
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
715 END SUBROUTINE spherical_p_ns
720 SUBROUTINE spherical_ns (self, mem, v, m, iv, maxfill)
722 real(kind=KindScalarVar),
dimension(:,:,:) :: mem
724 class(mesh_t),
pointer :: m(:)
726 logical,
optional :: maxfill
728 integer :: ix, iy, iz
729 real(8) :: x, y, z, r, h(3)
730 real(8) :: d1, dw, x1, rp, p1, p2
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 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 748 rp = sqrt(abs(self%radius**2 - y**2 - z**2))
750 x1 = m(1)%p + r1(ix-1) + h(1)*m(1)%d - self%position(1)
753 mem(ix-1,iy,iz) = mem(ix-2,iy,iz) * dw
755 extrapolate_antisym(x,x1,rp,mem(ix-1,iy,iz))
756 else if ((p1 <= self%radius+m(1)%d*0.5).and.(p2 >= self%radius+m(1)%d*0.5))
then 758 rp = sqrt(abs(self%radius**2 - y**2 - z**2))
760 x1 = m(1)%p + r1(ix+1) + h(1)*m(1)%d - self%position(1)
763 mem(ix+1,iy,iz) = mem(ix+2,iy,iz)*dw
765 extrapolate_antisym(x,x1,rp,mem(ix+1,iy,iz))
766 else if ((p1 <= self%radius+m(1)%d*0.5).and.(p2 <= self%radius+m(1)%d*0.5))
then 769 rp = sqrt(abs(self%radius**2 - y**2 - z**2))
771 x1 = m(1)%p + r1(ix-1) + h(1)*m(1)%d - self%position(1)
774 mem(ix-1,iy,iz) = mem(ix-2,iy,iz) * dw
775 x1 = m(1)%p + r1(ix+1) + h(1)*m(1)%d - self%position(1)
778 mem(ix+1,iy,iz) = mem(ix+2,iy,iz) * dw
782 else if (iv == idx%py)
then 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 788 rp = sqrt(abs(self%radius**2 - x**2 - z**2))
790 x1 = m(2)%p + r2(iy-1) + h(2)*m(2)%d - self%position(2)
793 mem(ix,iy-1,iz) = mem(ix,iy-2,iz) * dw
795 extrapolate_antisym(y,x1,rp,mem(ix,iy-1,iz))
796 else if ((p1 < self%radius+m(2)%d*0.5).and.(p2 > self%radius+m(2)%d*0.5))
then 798 rp = sqrt(abs(self%radius**2 - x**2 - z**2))
800 x1 = m(2)%p + r2(iy+1) + h(2)*m(2)%d - self%position(2)
803 mem(ix,iy+1,iz) = mem(ix,iy+2,iz)*dw
805 extrapolate_antisym(y,x1,rp,mem(ix,iy+1,iz))
806 else if ((p1 < self%radius+m(2)%d*0.5).and.(p2 < self%radius+m(2)%d*0.5))
then 809 rp = sqrt(abs(self%radius**2 - x**2 - z**2))
811 x1 = m(2)%p + r2(iy-1) + h(2)*m(2)%d - self%position(2)
814 mem(ix,iy-1,iz) = mem(ix,iy-2,iz) * dw
815 x1 = m(2)%p + r2(iy+1) + h(2)*m(2)%d - self%position(2)
818 mem(ix,iy+1,iz) = mem(ix,iy+2,iz) * dw
822 else if (iv == idx%pz)
then 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 828 rp = sqrt(abs(self%radius**2 - x**2 - y**2))
830 x1 = m(3)%p + r3(iz-1) + h(3)*m(3)%d - self%position(3)
833 mem(ix,iy,iz-1) = mem(ix,iy,iz-2) * dw
835 extrapolate_antisym(z,x1,rp,mem(ix,iy,iz-1))
836 else if ((p1 < self%radius+m(3)%d*0.5).and.(p2 > self%radius+m(3)%d*0.5))
then 838 rp = sqrt(abs(self%radius**2 - x**2 - y**2))
840 x1 = m(3)%p + r3(iz+1) + h(3)*m(3)%d - self%position(3)
843 mem(ix,iy,iz+1) = mem(ix,iy,iz+2)*dw
845 extrapolate_antisym(z,x1,rp,mem(ix,iy,iz+1))
846 else if ((p1 < self%radius+m(3)%d*0.5).and.(p2 < self%radius+m(3)%d*0.5))
then 849 rp = sqrt(abs(self%radius**2 - x**2 - y**2))
851 x1 = m(3)%p + r3(iz-1) + h(3)*m(3)%d - self%position(3)
854 mem(ix,iy,iz-1) = mem(ix,iy,iz-2) * dw
855 x1 = m(3)%p + r3(iz+1) + h(3)*m(3)%d - self%position(3)
858 mem(ix,iy,iz+1) = mem(ix,iy,iz+2) * dw
869 END SUBROUTINE spherical_ns
874 SUBROUTINE spherical_jagged_d (self, mem, m, iv, it, new)
876 real(kind=KindScalarVar),
dimension(:,:,:,:) :: mem
877 class(mesh_t),
pointer :: m(:), m1, m2, m3
878 integer :: iv, it, new
880 integer :: ix, iy, iz
881 real(8) :: x1, x2, y1, y2, z1, z2
885 call trace_begin(
'boundaries_t%spherical_jagged_d')
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) = &
904 END SUBROUTINE spherical_jagged_d
909 SUBROUTINE spherical_jagged_s (self, new, old, m)
911 real(kind=KindScalarVar),
dimension(:,:,:) :: old, new
912 class(mesh_t),
pointer :: m(:), m1, m2, m3
915 integer :: ix, iy, iz
919 call trace_begin(
'boundaries_t%spherical_jagged_s')
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) = &
938 END SUBROUTINE spherical_jagged_s
943 SUBROUTINE spherical_jagged_s0 (self, mem, m, iv, v)
945 real(kind=KindScalarVar),
dimension(:,:,:) :: mem
946 class(mesh_t),
pointer :: m(:), m1, m2, m3
950 integer :: ix, iy, iz
951 real(8) :: x1, x2, y1, y2, z1, z2
955 call trace_begin(
'boundaries_t%spherical_jagged_s0')
960 if (iv == idx%s)
then 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
983 print *,
'iv passed',iv,
'iv required', idx%s
984 call mpi%abort (
"boundaries_mod::spherical_jagged_d:: invalid memory index iv")
988 END SUBROUTINE spherical_jagged_s0
993 SUBROUTINE spherical_jagged_p (self, mem, m, iv)
995 real(kind=KindScalarVar),
dimension(:,:,:) :: mem
996 class(mesh_t),
pointer :: m(:), m1, m2, m3
999 integer :: ix, iy, iz
1000 real(8) :: x1, x2, y1, y2, z1, z2
1002 logical :: inside(4)
1004 call trace_begin(
'boundaries_t%spherical_jagged_p')
1008 call set_stagger (m, h, iv)
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
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
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
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")
1065 nullify (m1, m2, m3)
1067 END SUBROUTINE spherical_jagged_p
1072 FUNCTION extrapolate_antisym(ip, op, bp, yop)
result(out)
1073 real(8) :: ip, op, bp
1074 real(kind=KindScalarVar):: yop, out
1076 out = yop+(ip-op)/(bp-op)*(-yop)
1077 END FUNCTION extrapolate_antisym
1082 SUBROUTINE spherical3 (self, mem, v, m, iv, maxfill)
1084 real(kind=KindScalarVar),
dimension(:,:,:) :: mem, v
1085 class(mesh_t),
pointer :: m(:)
1087 logical,
optional :: maxfill
1089 integer :: ix, iy, iz, n1, n2
1090 real(8) :: x, y, z, r, h(3)
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)
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))
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)
1128 END SUBROUTINE spherical3
1134 SUBROUTINE spherical_extrapolate (self, mem, m, iv)
1136 real(kind=KindScalarVar),
dimension(:,:,:) :: mem
1137 class(mesh_t),
pointer :: m(:)
1140 integer :: ix, iy, iz, n1, n2
1141 real(8) :: x, y, z, r, h(3), a1, a2, ds, p
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)
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 1161 a1 = a1 + mem(ix,iy,iz)
1162 else if (abs(r-self%radius) <= ds)
then 1164 a2 = a2 + mem(ix,iy,iz)
1169 if (n1>0) a1 = log(a1/n1)
1170 if (n2>0) a2 = log(a2/n2)
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)
1195 END SUBROUTINE spherical_extrapolate
1200 SUBROUTINE spherical_constant_entropy (self, d, e, m, s, gamma)
1202 real(kind=KindScalarVar),
dimension(:,:,:) :: d, e
1203 class(mesh_t),
pointer :: m(:)
1206 integer :: ix, iy, iz
1207 real(8) :: x, y, z, r, ds, c
1209 call trace_begin(
'boundaries_t%sperical_entropy')
1213 c = exp(s*(gamma-1.0))/(gamma-1.0)
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
1232 END SUBROUTINE spherical_constant_entropy
Boundary conditions for centered variables, which have the physical boundary in-between mesh points l...
Template module for mesh.
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...