23 integer,
intent(in) :: ixi^
l, ix^
l, idims
24 double precision,
intent(in) :: q(ixi^s),qct(ixi^s)
26 double precision,
intent(inout) :: qrc(ixg^t),qlc(ixg^t)
28 double precision,
dimension(ixG^T) :: dqc,d2qc,ldq
29 double precision,
dimension(ixG^T) :: qmin,qmax,tmp
31 integer :: lxc^
l,lxr^
l
32 integer :: ixll^
l,ixl^
l,ixo^
l,ixr^
l,ixrr^
l
33 integer :: hxl^
l,hxc^
l,hxr^
l
34 integer :: kxll^
l,kxl^
l,kxc^
l,kxr^
l,kxrr^
l
36 ixomin^
d=ixmin^
d-
kr(idims,^
d);ixomax^
d=ixmax^
d+
kr(idims,^
d);
37 ixl^
l=ixo^
l-
kr(idims,^
d);
38 ixll^
l=ixl^
l-
kr(idims,^
d);
39 ixr^
l=ixo^
l+
kr(idims,^
d);
40 ixrr^
l=ixr^
l+
kr(idims,^
d);
42 hxcmin^
d=ixomin^
d;hxcmax^
d=ixmax^
d;
43 hxl^
l=hxc^
l-
kr(idims,^
d);
44 hxr^
l=hxc^
l+
kr(idims,^
d);
46 kxcmin^
d=ixllmin^
d; kxcmax^
d=ixrmax^
d;
47 kxl^
l=kxc^
l-
kr(idims,^
d);
48 kxr^
l=kxc^
l+
kr(idims,^
d);
50 lxcmin^
d=ixllmin^
d-
kr(idims,^
d);lxcmax^
d=ixrrmax^
d;
51 lxr^
l=lxc^
l+
kr(idims,^
d);
54 dqc(lxc^s)=q(lxr^s)-q(lxc^s)
56 d2qc(kxc^s)=half*(q(kxr^s)-q(kxl^s))
57 where(dqc(kxc^s)*dqc(kxl^s)>zero)
59 qmin(kxc^s)= sign(one,d2qc(kxc^s))
61 ldq(kxc^s)= qmin(kxc^s)*min(dabs(d2qc(kxc^s)),&
62 2.0d0*dabs(dqc(kxl^s)),&
63 2.0d0*dabs(dqc(kxc^s)))
69 qlc(ixo^s)=qlc(ixo^s)+half*dqc(ixo^s)&
70 +(ldq(ixo^s)-ldq(ixr^s))/6.0d0
71 qrc(ixl^s)=qrc(ixl^s) -(half*dqc(ixl^s)&
72 -(ldq(ixl^s)-ldq(ixo^s))/6.0d0)
78 qrc(ixl^s)=max(qmin(ixo^s),min(qmax(ixo^s),qrc(ixl^s)))
79 qlc(ixo^s)=max(qmin(ixo^s),min(qmax(ixo^s),qlc(ixo^s)))
82 where((qrc(ixl^s)-qct(ixo^s))*(qct(ixo^s)-qlc(ixo^s))<=zero)
87 qmax(ixo^s)=(qlc(ixo^s)-qrc(ixl^s))&
88 *(qct(ixo^s)-(qlc(ixo^s)+qrc(ixl^s))/2.0d0)
89 qmin(ixo^s)=(qlc(ixo^s)-qrc(ixl^s))**2.0d0/6.0d0
93 where(qmax(hxr^s)>qmin(hxr^s))
94 qrc(hxc^s)= 3.0d0*qct(hxr^s)-2.0d0*qlc(hxr^s)
97 where(qmax(hxc^s)<-qmin(hxc^s))
98 qlc(hxc^s)= 3.0d0*qct(hxc^s)-2.0d0*tmp(hxl^s)
114 use mod_physics,
only: phys_ppm_flatcd, phys_ppm_flatsh
116 integer,
intent(in) :: ixi^
l, ix^
l, idims
117 double precision,
intent(in) :: w(ixi^s,1:nw),wct(ixi^s,1:nw)
119 double precision,
intent(inout) :: wrc(ixg^t,1:nw),wlc(ixg^t,1:nw)
121 double precision,
dimension(ixG^T,1:nwflux) :: dwc,d2wc,ldw
122 double precision,
dimension(ixG^T,1:nwflux) :: wmin,wmax,tmp
123 double precision,
dimension(ixG^T) :: aa, ab, ac, dv
124 double precision,
dimension(ixG^T,1:ndim) :: exi
126 integer :: lxc^
l,lxr^
l
127 integer :: ixll^
l,ixl^
l,ixo^
l,ixr^
l,ixrr^
l
128 integer :: hxl^
l,hxc^
l,hxr^
l
129 integer :: kxll^
l,kxl^
l,kxc^
l,kxr^
l,kxrr^
l
130 integer :: iw, idimss
132 double precision,
parameter :: betamin=0.75d0, betamax=0.85d0,&
133 zmin=0.25d0, zmax=0.75d0,&
134 eta1=20.0d0,eta2=0.05d0,eps=0.01d0,kappa=0.1d0
136 ixomin^
d=ixmin^
d-
kr(idims,^
d);ixomax^
d=ixmax^
d+
kr(idims,^
d);
137 ixl^
l=ixo^
l-
kr(idims,^
d);
138 ixll^
l=ixl^
l-
kr(idims,^
d);
139 ixr^
l=ixo^
l+
kr(idims,^
d);
140 ixrr^
l=ixr^
l+
kr(idims,^
d);
142 hxcmin^
d=ixomin^
d;hxcmax^
d=ixmax^
d;
143 hxl^
l=hxc^
l-
kr(idims,^
d);
144 hxr^
l=hxc^
l+
kr(idims,^
d);
146 kxcmin^
d=ixllmin^
d; kxcmax^
d=ixrmax^
d;
147 kxl^
l=kxc^
l-
kr(idims,^
d);
148 kxr^
l=kxc^
l+
kr(idims,^
d);
150 lxcmin^
d=ixllmin^
d-
kr(idims,^
d);lxcmax^
d=ixrrmax^
d;
151 lxr^
l=lxc^
l+
kr(idims,^
d);
153 dwc(lxc^s,1:nwflux)=w(lxr^s,1:nwflux)-w(lxc^s,1:nwflux)
155 d2wc(kxc^s,1:nwflux)=half*(w(kxr^s,1:nwflux)-w(kxl^s,1:nwflux))
156 where(dwc(kxc^s,1:nwflux)*dwc(kxl^s,1:nwflux)>zero)
158 wmin(kxc^s,1:nwflux)= sign(one,d2wc(kxc^s,1:nwflux))
160 ldw(kxc^s,1:nwflux)= wmin(kxc^s,1:nwflux)*min(dabs(d2wc(kxc^s,1:nwflux)),&
161 2.0d0*dabs(dwc(kxl^s,1:nwflux)),&
162 2.0d0*dabs(dwc(kxc^s,1:nwflux)))
164 ldw(kxc^s,1:nwflux)=zero
168 wlc(ixo^s,1:nwflux)=wlc(ixo^s,1:nwflux)+half*dwc(ixo^s,1:nwflux)&
169 +(ldw(ixo^s,1:nwflux)-ldw(ixr^s,1:nwflux))/6.0d0
171 wrc(ixl^s,1:nwflux)=wrc(ixl^s,1:nwflux)-(half*dwc(ixl^s,1:nwflux)&
172 -(ldw(ixl^s,1:nwflux)-ldw(ixo^s,1:nwflux))/6.0d0)
178 wrc(ixl^s,1:nwflux)=max(wmin(ixo^s,1:nwflux)&
179 ,min(wmax(ixo^s,1:nwflux),wrc(ixl^s,1:nwflux)))
180 wlc(ixo^s,1:nwflux)=max(wmin(ixo^s,1:nwflux)&
181 ,min(wmax(ixo^s,1:nwflux),wlc(ixo^s,1:nwflux)))
185 where((wrc(ixl^s,1:nwflux)-wct(ixo^s,1:nwflux))&
186 *(wct(ixo^s,1:nwflux)-wlc(ixo^s,1:nwflux))<=zero)
187 wrc(ixl^s,1:nwflux)=wct(ixo^s,1:nwflux)
188 wlc(ixo^s,1:nwflux)=wct(ixo^s,1:nwflux)
191 wmax(ixo^s,1:nwflux)=(wlc(ixo^s,1:nwflux)-wrc(ixl^s,1:nwflux))*&
192 (wct(ixo^s,1:nwflux)-&
193 (wlc(ixo^s,1:nwflux)+wrc(ixl^s,1:nwflux))/2.0d0)
194 wmin(ixo^s,1:nwflux)=(wlc(ixo^s,1:nwflux)-wrc(ixl^s,1:nwflux))**2.0d0/6.0d0
195 tmp(ixl^s,1:nwflux)=wrc(ixl^s,1:nwflux)
197 where(wmax(hxr^s,1:nwflux)>wmin(hxr^s,1:nwflux))
198 wrc(hxc^s,1:nwflux)= 3.0d0*wct(hxr^s,1:nwflux)&
199 -2.0d0*wlc(hxr^s,1:nwflux)
202 where(wmax(hxc^s,1:nwflux)<-wmin(hxc^s,1:nwflux))
203 wlc(hxc^s,1:nwflux)= 3.0d0*wct(hxc^s,1:nwflux)&
204 -2.0d0*tmp(hxl^s,1:nwflux)
209 call phys_ppm_flatcd(ixi^
l,kxc^
l,kxl^
l,kxr^
l,wct,d2wc,aa,ab)
210 if(any(kappa*aa(kxc^s)>=ab(kxc^s)))
then
212 where(kappa*aa(kxc^s)>=ab(kxc^s).and. dabs(dwc(kxc^s,iw))>smalldouble)
213 wmax(kxc^s,iw) = wct(kxr^s,iw)-2.0d0*wct(kxc^s,iw)+wct(kxl^s,iw)
216 where(wmax(ixr^s,iw)*wmax(ixl^s,iw)<zero .and.&
217 dabs(wct(ixr^s,iw)-wct(ixl^s,iw))&
218 -eps*min(dabs(wct(ixr^s,iw)),dabs(wct(ixl^s,iw)))>zero &
219 .and. kappa*aa(ixo^s)>=ab(ixo^s)&
220 .and. dabs(dwc(ixo^s,iw))>smalldouble)
222 ac(ixo^s)=(wct(ixll^s,iw)-wct(ixrr^s,iw)+4.0d0*dwc(ixo^s,iw))&
223 /(12.0d0*dwc(ixo^s,iw))
224 wmin(ixo^s,iw)=max(zero,min(eta1*(ac(ixo^s)-eta2),one))
229 where(wmin(hxc^s,iw)>zero)
230 wlc(hxc^s,iw) = wlc(hxc^s,iw)*(one-wmin(hxc^s,iw))&
231 +(wct(hxc^s,iw)+half*ldw(hxc^s,iw))*wmin(hxc^s,iw)
233 where(wmin(hxr^s,iw)>zero)
234 wrc(hxc^s,iw) = wrc(hxc^s,iw)*(one-wmin(hxr^s,iw))&
235 +(wct(hxr^s,iw)-half*ldw(hxr^s,iw))*wmin(hxr^s,iw)
244 kxcmin^
d=ixmin^
d-2; kxcmax^
d=ixmax^
d+2;
246 kxl^
l=kxc^
l-
kr(idimss,^
d);
247 kxr^
l=kxc^
l+
kr(idimss,^
d);
248 kxll^
l=kxl^
l-
kr(idimss,^
d);
249 kxrr^
l=kxr^
l+
kr(idimss,^
d);
251 call phys_ppm_flatsh(ixi^
l,kxc^
l,kxll^
l,kxl^
l,kxr^
l,kxrr^
l,idimss,wct,aa,ab,dv)
254 ac(kxc^s) = max(zero,min(one,(betamax-aa(kxc^s))/(betamax-betamin)))
257 where (dabs(dv(kxc^s))<smalldouble)
258 aa(kxc^s) = max(ac(kxc^s), min(one,(zmax-ab(kxc^s))/(zmax-zmin)))
264 {^nooned
call extremaa(ixi^
l,ixo^
l,aa,1,exi(ixg^t,idimss))}
266 {^nooned ab(ixo^s)=min(^
d&exi(ixo^s,^
d))}
269 where(dabs(ab(ixo^s)-one)>smalldouble)
270 wmax(ixo^s,iw) = (one-ab(ixo^s))*wct(ixo^s,iw)
273 where(dabs(ab(hxc^s)-one)>smalldouble)
274 wlc(hxc^s,iw) = ab(hxc^s)*wlc(hxc^s,iw)+wmax(hxc^s,iw)
277 where(dabs(ab(hxr^s)-one)>smalldouble)
278 wrc(hxc^s,iw) = ab(hxr^s)*wrc(hxc^s,iw)+wmax(hxr^s,iw)
285 subroutine extremaq(ixI^L,ixO^L,q,nshift,qMax,qMin)
289 integer,
intent(in) :: ixI^L,ixO^L
290 double precision,
intent(in) :: q(ixI^S)
291 integer,
intent(in) :: nshift
293 double precision,
intent(out) :: qMax(ixI^S),qMin(ixI^S)
295 integer :: ixs^L,ixsR^L,ixsL^L,idims,jdims,kdims,ishift,i,j
299 ixsr^l=ixo^l+ishift*
kr(idims,^
d);
300 ixsl^l=ixo^l-ishift*
kr(idims,^
d);
302 qmax(ixo^s)=max(q(ixo^s),q(ixsr^s),q(ixsl^s))
303 qmin(ixo^s)=min(q(ixo^s),q(ixsr^s),q(ixsl^s))
305 qmax(ixo^s)=max(qmax(ixo^s),q(ixsr^s),q(ixsl^s))
306 qmin(ixo^s)=min(qmin(ixo^s),q(ixsr^s),q(ixsl^s))
312 ixs^l=ixo^l+i*ishift*
kr(idims,^
d);
313 ixsr^l=ixs^l+ishift*
kr(jdims,^
d);
314 ixsl^l=ixs^l-ishift*
kr(jdims,^
d);
315 qmax(ixo^s)=max(qmax(ixo^s),q(ixsr^s),q(ixsl^s))
316 qmin(ixo^s)=min(qmin(ixo^s),q(ixsr^s),q(ixsl^s))
324 ixs^l=ixo^l+i*ishift*
kr(idims,^
d);
326 ixs^l=ixo^l+j*ishift*
kr(jdims,^
d);
327 ixsr^l=ixs^l+ishift*
kr(kdims,^
d);
328 ixsl^l=ixs^l-ishift*
kr(kdims,^
d);
329 qmax(ixo^s)=max(qmax(ixo^s),q(ixsr^s),q(ixsl^s))
330 qmin(ixo^s)=min(qmin(ixo^s),q(ixsr^s),q(ixsl^s))
341 integer,
intent(in) :: ixI^L,ixO^L
342 double precision,
intent(in) :: a(ixI^S)
343 integer,
intent(in) :: nshift
345 double precision,
intent(out) :: aMin(ixI^S)
347 integer :: ixs^L,ixsR^L,ixsL^L,idims,jdims,kdims,ishift,i,j
351 ixsr^l=ixo^l+ishift*
kr(idims,^
d);
352 ixsl^l=ixo^l-ishift*
kr(idims,^
d);
353 amin(ixo^s)=min(a(ixsr^s),a(ixo^s),a(ixsl^s))
358 ixs^l=ixo^l+i*ishift*
kr(idims,^
d);
359 ixsr^l=ixs^l+ishift*
kr(jdims,^
d);
360 ixsl^l=ixs^l-ishift*
kr(jdims,^
d);
361 amin(ixo^s)=min(amin(ixo^s),a(ixsr^s),a(ixsl^s))
369 ixs^l=ixo^l+i*ishift*
kr(idims,^
d);
371 ixs^l=ixo^l+j*ishift*
kr(jdims,^
d);
372 ixsr^l=ixs^l+ishift*
kr(kdims,^
d);
373 ixsl^l=ixs^l-ishift*
kr(kdims,^
d);
374 amin(ixo^s)=min(amin(ixo^s),a(ixsr^s),a(ixsl^s))
382 subroutine extremaw(ixI^L,ixO^L,w,nshift,wMax,wMin)
385 integer,
intent(in) :: ixI^L,ixO^L
386 double precision,
intent(in) :: w(ixI^S,1:nw)
387 integer,
intent(in) :: nshift
389 double precision,
intent(out) :: wMax(ixI^S,1:nwflux),wMin(ixI^S,1:nwflux)
391 integer :: ixs^L,ixsR^L,ixsL^L,idims,jdims,kdims,ishift,i,j
395 ixsr^l=ixo^l+ishift*
kr(idims,^
d);
396 ixsl^l=ixo^l-ishift*
kr(idims,^
d);
398 wmax(ixo^s,1:nwflux)= &
399 max(w(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
400 wmin(ixo^s,1:nwflux)= &
401 min(w(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
403 wmax(ixo^s,1:nwflux)= &
404 max(wmax(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
405 wmin(ixo^s,1:nwflux)= &
406 min(wmin(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
412 ixs^l=ixo^l+i*ishift*
kr(idims,^
d);
413 ixsr^l=ixs^l+ishift*
kr(jdims,^
d);
414 ixsl^l=ixs^l-ishift*
kr(jdims,^
d);
415 wmax(ixo^s,1:nwflux)= &
416 max(wmax(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
417 wmin(ixo^s,1:nwflux)= &
418 min(wmin(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
426 ixs^l=ixo^l+i*ishift*
kr(idims,^
d);
428 ixs^l=ixo^l+j*ishift*
kr(jdims,^
d);
429 ixsr^l=ixs^l+ishift*
kr(kdims,^
d);
430 ixsl^l=ixs^l-ishift*
kr(kdims,^
d);
431 wmax(ixo^s,1:nwflux)= &
432 max(wmax(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
433 wmin(ixo^s,1:nwflux)= &
434 min(wmin(ixo^s,1:nwflux),w(ixsr^s,1:nwflux),w(ixsl^s,1:nwflux))
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer, dimension(:), allocatable, parameter d
This module defines the procedures of a physics module. It contains function pointers for the various...
subroutine, public ppmlimitervar(ixIL, ixL, idims, q, qCT, qLC, qRC)
subroutine extremaq(ixIL, ixOL, q, nshift, qMax, qMin)
subroutine extremaw(ixIL, ixOL, w, nshift, wMax, wMin)
subroutine, public ppmlimiter(ixIL, ixL, idims, w, wCT, wLC, wRC)
subroutine extremaa(ixIL, ixOL, a, nshift, aMin)