37 character(len=10):: riemann=
'hlld' 38 real,
pointer,
dimension(:,:,:,:):: q
39 real,
pointer,
dimension(:,:,:,:):: bf, em
40 real,
pointer,
dimension(:,:,:,:,:):: qp, qm
41 real,
pointer,
dimension(:,:,:,:,:):: dq, dq2, dbf
42 real,
pointer,
dimension(:,:,:,:,:):: qlb,qrb,qlt,qrt
43 real,
pointer,
dimension(:,:,:,:):: uin
44 real,
pointer,
dimension(:,:,:,:):: unew
45 real,
pointer,
dimension(:,:,:):: ex,ey,ez
46 real,
pointer,
dimension(:,:,:):: c
50 procedure:: gas_pressure
52 logical,
save:: first_time=.true.
53 logical,
save:: detailed_timer=.false.
54 logical,
save:: unsigned=.true.
55 integer,
save:: verbose=0
56 integer,
save:: divb_cleaner=2
57 integer,
save:: b_slope_type=0
58 integer,
save:: u_slope_type=0
59 character(len=16):: eqn=
'mhd' 63 SUBROUTINE init (self)
65 character(len=10):: riemann=
'hlld' 66 real,
save:: max_dlogd=20.
67 integer,
save:: n_dump=0
70 namelist /ramses_params/ gamma, riemann, slope_type, b_slope_type, &
71 u_slope_type, courant_factor, smallr, smallc, isothermal, csound, &
72 verbose, detailed_timer, unsigned, courant_type, max_dlogd, n_dump, &
74 character(len=120):: id = &
75 '$Id: 4b4ff75d6255b6cfcb2ea17e1d438c0ccadfb8ba $ solvers/ramses/mhd/mhd_mod.f90' 77 call trace%begin(
'mhd_mod%init')
78 call trace%print_id (id)
82 self%kind =
'ramses_mhd_patch' 83 call self%idx%init (5, self%mhd)
84 call self%patch_t%init
92 read (io_unit%input, ramses_params)
93 isothermal = isothermal .or. gamma==1d0
97 if (isothermal) gamma=1.0d0
98 if (mpi%master)
write (*, ramses_params)
101 rieman%max_dlogd = max_dlogd
103 self%riemann = riemann
104 self%courant = courant_factor
107 self%staggered = .false.
110 self%mesh(i)%h(self%idx%bx+i-1) = -0.5
113 self%unsigned(self%idx%d) = unsigned
114 self%unsigned(self%idx%e) = unsigned
115 self%pervolume(self%idx%px) = unsigned
116 self%pervolume(self%idx%py) = unsigned
117 self%pervolume(self%idx%pz) = unsigned
118 call self%gpatch_t%init
119 call self%extras_t%init
125 SUBROUTINE update (self)
127 real,
pointer,
dimension(:,:,:,:):: p
128 real,
allocatable,
dimension(:,:,:,:):: em
129 real,
allocatable,
dimension(:,:,:,:,:):: fl
130 real,
allocatable,
dimension(:,:,:,:)::fluxambdiff,emfambdiff
131 real,
pointer,
dimension(:,:,:):: d, e, ux=>null(), uy=>null(), uz=>null()
132 real(dp):: dtdx, dtdy, dtdz
133 integer:: n(3), nv3, i
134 integer,
parameter:: ndim=3
135 integer,
save:: itimer=0
137 call trace%begin (
'mhd_t%update', itimer=itimer)
138 n = self%mesh%ui-self%mesh%li+1
139 if (io%smallr > 0) smallr = io%smallr
140 if (io%courant > 0) self%courant = io%courant
151 self%uin => self%mem(:,:,:,:,self%it ,1)
152 self%unew => self%mem(:,:,:,:,self%new,1)
155 allocate (self%c (n(1),n(2),n(3)))
156 allocate (self%Ex (n(1),n(2),n(3)))
157 allocate (self%Ey (n(1),n(2),n(3)))
158 allocate (self%Ez (n(1),n(2),n(3)))
159 allocate (self%q (n(1),n(2),n(3),self%nv))
160 allocate (self%qp (n(1),n(2),n(3),self%nv,ndim))
161 allocate (self%qm (n(1),n(2),n(3),self%nv,ndim))
162 allocate (self%dq (n(1),n(2),n(3),self%nv,ndim))
163 allocate (self%dq2(n(1),n(2),n(3),self%nv,ndim))
164 allocate (self%bf (n(1),n(2),n(3),3))
165 allocate (self%dbf(n(1),n(2),n(3),3,2))
166 allocate (self%qLB(n(1),n(2),n(3),self%nv,ndim))
167 allocate (self%qLT(n(1),n(2),n(3),self%nv,ndim))
168 allocate (self%qRB(n(1),n(2),n(3),self%nv,ndim))
169 allocate (self%qRT(n(1),n(2),n(3),self%nv,ndim))
170 d => self%mem(:,:,:, self%idx%d ,self%it,1)
171 e => self%mem(:,:,:, self%idx%e ,self%it,1)
172 p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%it,1)
176 call check_small (self)
179 call uslope(self,self%bf,self%q,self%dq,self%dbf,slope_type)
180 if (verbose>0) print *,
'uslope', self%id, minval(self%dbf), maxval(self%dbf)
185 if (slope_type <= 3.0)
then 188 call uslope(self,self%bf,self%q,self%dq2,self%dbf,3.0)
191 if (verbose>0) print *,
'uslope', self%id, minval(self%dq), maxval(self%dq)
192 if (verbose>0) print *,
'uslope', self%id, minval(self%dq2), maxval(self%dq2)
194 if (non_ideal%is_used)
then 195 allocate (fluxambdiff(self%gn(1),self%gn(2),self%gn(3),3))
196 allocate (emfambdiff(self%gn(1),self%gn(2),self%gn(3),3))
197 call non_ideal%non_ideal_comp(self%mesh,self%idx,self%bf,self%uin,self%q,&
198 fluxambdiff,emfambdiff,self%u_max,self%ds,self%dtime,self%gamma)
200 call self%courant_condition
204 allocate (fl(self%gn(1),self%gn(2),self%gn(3),self%nv,3))
205 allocate (em(self%gn(1),self%gn(2),self%gn(3),3))
208 call cmpflxm(self,self%gn(1),2,3,4,6,7,8,fl(:,:,:,:,1))
209 call cmpflxm(self,self%gn(2),3,2,4,7,6,8,fl(:,:,:,:,2))
210 call cmpflxm(self,self%gn(3),4,2,3,8,6,7,fl(:,:,:,:,3))
212 call cmp_mag_flx(self, self%qRT, self%qRB, self%qLT, self%qLB,2,3,4,6,7,8,em(:,:,:,3),self%gn(3)-1,self%nv)
213 call cmp_mag_flx(self, self%qRT, self%qLT, self%qRB, self%qLB,4,2,3,8,6,7,em(:,:,:,2),self%gn(2)-1,self%nv)
214 call cmp_mag_flx(self, self%qRT, self%qRB, self%qLT, self%qLB,3,4,2,7,8,6,em(:,:,:,1),self%gn(1)-1,self%nv)
216 dtdx = self%dtime/self%ds(1)
217 dtdy = self%dtime/self%ds(2)
218 dtdz = self%dtime/self%ds(3)
219 em(:,:,:,1) = em(:,:,:,1)*dtdy
220 em(:,:,:,2) = em(:,:,:,2)*dtdz
221 em(:,:,:,3) = em(:,:,:,3)*dtdx
223 if (non_ideal%is_used)
then 224 call non_ideal%flux_upd(fl,fluxambdiff,self%ds(1),self%dtime,1)
225 call non_ideal%flux_upd(fl,fluxambdiff,self%ds(2),self%dtime,2)
226 call non_ideal%flux_upd(fl,fluxambdiff,self%ds(3),self%dtime,3)
227 call non_ideal%non_ideal_emf_up(em(:,:,:,3),emfambdiff,self%ds(3),self%dtime,3)
228 call non_ideal%non_ideal_emf_up(em(:,:,:,2),emfambdiff,self%ds(2),self%dtime,2)
229 call non_ideal%non_ideal_emf_up(em(:,:,:,1),emfambdiff,self%ds(1),self%dtime,1)
232 fl(:,:,:,:,1) = fl(:,:,:,:,1) * dtdx
233 fl(:,:,:,:,2) = fl(:,:,:,:,2) * dtdy
234 fl(:,:,:,:,3) = fl(:,:,:,:,3) * dtdz
239 associate(i0=>self%mesh(1)%li, j0=>self%mesh(2)%li, k0=>self%mesh(3)%li, &
240 i1=>self%mesh(1)%ui, j1=>self%mesh(2)%ui, k1=>self%mesh(3)%ui, &
241 uin=>self%uin,unew=>self%unew)
243 call self%lock%set (
'update')
246 unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,1:nv3 ) = &
247 unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,1:nv3 ) + &
248 fl(i0 :i1 ,j0 :j1 ,k0 :k1 ,1:nv3,1) - &
249 fl(i0+1:i1+1,j0 :j1 ,k0 :k1 ,1:nv3,1) + &
250 fl(i0 :i1 ,j0 :j1 ,k0 :k1 ,1:nv3,2) - &
251 fl(i0 :i1 ,j0+1:j1+1,k0 :k1 ,1:nv3,2) + &
252 fl(i0 :i1 ,j0 :j1 ,k0 :k1 ,1:nv3,3) - &
253 fl(i0 :i1 ,j0 :j1 ,k0+1:k1+1,1:nv3,3)
266 unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,6) = &
267 unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,6) + &
268 ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,2) - &
269 em(i0 :i1 ,j0 :j1 ,k0+1:k1+1,2) ) - &
270 ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,3) - &
271 em(i0 :i1 ,j0+1:j1+1,k0 :k1 ,3) )
276 unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,7) = &
277 unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,7) + &
278 ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,3) - &
279 em(i0+1:i1+1,j0 :j1 ,k0 :k1 ,3) ) - &
280 ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,1) - &
281 em(i0 :i1 ,j0 :j1 ,k0+1:k1+1,1) )
286 unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,8) = &
287 unew(i0 :i1 ,j0 :j1 ,k0 :k1 ,8) + &
288 ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,1) - &
289 em(i0 :i1 ,j0+1:j1+1,k0 :k1 ,1) ) - &
290 ( em(i0 :i1 ,j0 :j1 ,k0 :k1 ,2) - &
291 em(i0+1:i1+1,j0 :j1 ,k0 :k1 ,2) )
293 if (verbose>0) print *,
'unew', self%id, minval(unew), maxval(unew)
299 self%mem(:,:,:,self%idx%e,self%new,1) = &
300 self%mem(:,:,:,self%idx%d,self%new,1)
304 if (
allocated(self%force_per_unit_mass))
then 305 p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%new,1)
307 p(:,:,:,i) = p(:,:,:,i) + self%dtime*self%force_per_unit_mass(:,:,:,i)*d
313 call self%counter_update
314 if (non_ideal%is_used) &
315 deallocate (fluxambdiff,emfambdiff)
316 deallocate (fl,em,self%q,self%qp,self%qm,self%dq,self%dq2,self%bf,self%dbf, &
317 self%qLT,self%qLB,self%qRT,self%qRB,self%c, &
318 self%Ex,self%Ey,self%Ez)
322 select case (divb_cleaner)
324 call divb_clean1 (self)
326 call divb_clean2 (self)
329 call self%lock%unset (
'update')
330 call trace%end (itimer)
331 if (rieman%ndiff > 0)
then 332 write(io_unit%mpi,*)
'id,ndiff', self%id, rieman%ndiff, rieman%nsolu
338 END SUBROUTINE update
351 SUBROUTINE divb_clean1 (self)
354 real,
dimension(:,:,:),
pointer:: bx, by, bz
355 real,
dimension(:,:,:),
allocatable:: divb
356 real,
parameter:: eps=0.11
359 call trace%begin(
'mhd_t%divb_clean1', itimer=itimer)
360 bx => self%mem(:,:,:,self%idx%bx,self%it,1)
361 by => self%mem(:,:,:,self%idx%by,self%it,1)
362 bz => self%mem(:,:,:,self%idx%bz,self%it,1)
363 allocate (divb(self%mesh(1)%gn,self%mesh(2)%gn,self%mesh(3)%gn))
364 divb = ddup(bx,1) + ddup(by,2) + ddup(bz,3)
365 bx = bx + eps*dddn(divb,1)
366 by = by + eps*dddn(divb,2)
367 bz = bz + eps*dddn(divb,3)
369 call trace%end (itimer)
373 function ddup(f,i)
result (g)
374 real:: f(:,:,:), g(size(f,1),size(f,2),size(f,3))
377 g = cshift(f,1,i) - f
380 function dddn(f,i)
result (g)
381 real:: f(:,:,:), g(size(f,1),size(f,2),size(f,3))
384 g = f - cshift(f,-1,i)
386 END SUBROUTINE divb_clean1
392 SUBROUTINE divb_clean2 (self)
395 real,
dimension(:,:,:),
pointer:: bx, by, bz
396 real,
dimension(:,:,:),
allocatable:: divb
397 real,
parameter:: eps=0.11
398 integer:: l(3), u(3), ix, iy, iz
401 call trace%begin(
'mhd_t%divb_clean2', itimer=itimer)
402 bx => self%mem(:,:,:,self%idx%bx,self%it,1)
403 by => self%mem(:,:,:,self%idx%by,self%it,1)
404 bz => self%mem(:,:,:,self%idx%bz,self%it,1)
408 do ix=self%mesh(1)%lo,self%mesh(1)%lb,-1
409 bx(ix ,l(2):u(2),l(3):u(3)) = &
410 bx(ix+1,l(2):u(2),l(3):u(3)) + ( &
411 by(ix,l(2)+1:u(2)+1,l(3) :u(3) ) - by(ix,l(2):u(2),l(3):u(3)) + &
412 bz(ix,l(2) :u(2) ,l(3)+1:u(3)+1) - bz(ix,l(2):u(2),l(3):u(3)))
414 do ix=self%mesh(1)%ui,self%mesh(1)%ub-1
415 bx(ix+1,l(2):u(2),l(3):u(3)) = &
416 bx(ix ,l(2):u(2),l(3):u(3)) - ( &
417 by(ix,l(2)+1:u(2)+1,l(3) :u(3) ) - by(ix,l(2):u(2),l(3):u(3)) + &
418 bz(ix,l(2) :u(2) ,l(3)+1:u(3)+1) - bz(ix,l(2):u(2),l(3):u(3)))
421 do iy=self%mesh(2)%lo,self%mesh(2)%lb,-1
422 by(l(1):u(1),iy ,l(3):u(3)) = &
423 by(l(1):u(1),iy+1,l(3):u(3)) + ( &
424 bx(l(1)+1:u(1)+1,iy,l(3) :u(3) ) - bx(l(1):u(1),iy,l(3):u(3)) + &
425 bz(l(1) :u(1) ,iy,l(3)+1:u(3)+1) - bz(l(1):u(1),iy,l(3):u(3)))
427 do iy=self%mesh(2)%ui,self%mesh(2)%ub-1
428 by(l(1):u(1),iy+1,l(3):u(3)) = &
429 by(l(1):u(1),iy ,l(3):u(3)) - ( &
430 bx(l(1)+1:u(1)+1,iy,l(3) :u(3) ) - bx(l(1):u(1),iy,l(3):u(3)) + &
431 bz(l(1) :u(1) ,iy,l(3)+1:u(3)+1) - bz(l(1):u(1),iy,l(3):u(3)))
434 do iz=self%mesh(3)%lo,self%mesh(3)%lb,-1
435 bz(l(1):u(1),l(2):u(2),iz ) = &
436 bz(l(1):u(1),l(2):u(2),iz+1) + ( &
437 bx(l(1)+1:u(1)+1,l(2) :u(2) ,iz) - bx(l(1):u(1),l(2):u(2),iz) + &
438 by(l(1) :u(1) ,l(2)+1:u(2)+1,iz) - by(l(1):u(1),l(2):u(2),iz))
440 do iz=self%mesh(3)%ui,self%mesh(3)%ub-1
441 bz(l(1):u(1),l(2):u(2),iz+1) = &
442 bz(l(1):u(1),l(2):u(2),iz ) - ( &
443 bx(l(1)+1:u(1)+1,l(2) :u(2) ,iz) - bx(l(1):u(1),l(2):u(2),iz) + &
444 by(l(1) :u(1) ,l(2)+1:u(2)+1,iz) - by(l(1):u(1),l(2):u(2),iz))
446 call trace%end (itimer)
447 END SUBROUTINE divb_clean2
451 subroutine check_small (self)
453 integer:: i,j,k,lo(3),uo(3),n_smalld,n_smallp
454 real(dp):: c,d,u,v,w,p,ekin,b1,b2,b3,emag
455 integer,
save:: nprint=10
456 integer,
save:: itimer=0
459 associate(uu=>self%uin)
461 call trace%begin(
'mhd_t%check_small', itimer=itimer, detailed_timer=detailed_timer)
470 if (uu(i,j,k,1) < smallr)
then 474 write (io_unit%log,
'(i7,3i5,1p,2e12.3)') self%id,i,j,k,uu(i,j,k,1),smallr
476 n_smalld = n_smalld+1
484 ekin = 0.5*d*(u**2+v**2+w**2)
485 b1 = 0.5*(uu(i,j,k,6)+uu(i+1,j,k,6))
486 b2 = 0.5*(uu(i,j,k,7)+uu(i,j+1,k,7))
487 b3 = 0.5*(uu(i,j,k,8)+uu(i,j,k+1,8))
488 emag = 0.5*(b1**2+b2**2+b3**2)
491 uu(i,j,k,5) = p+ekin+emag
493 p = (uu(i,j,k,5)-ekin)*(gamma-one)
494 if (p < smallr*smallc**2)
then 496 uu(i,j,k,5) = p/(gamma-one)+ekin+emag
497 n_smallp = n_smallp+1
503 if (n_smalld+n_smallp > 0)
then 504 write (io_unit%log,*)
'small: id,nr,np', &
505 wallclock(),self%id,self%time,n_smalld,n_smallp
508 if (panic .and. self%n_dump > 0)
then 510 self%n_dump = self%n_dump-1
511 write (io_unit%dump) self%id, self%gn, self%nv, self%nt, self%nw, self%istep
512 write (io_unit%dump) self%time, self%dtime
516 write (io_unit%dump) self%mem(:,:,:,:,jt,iw)
519 write (io_unit%mpi,*)
'density < SMALLR, dumped', self%id, self%time
522 call trace%end (itimer, detailed_timer=detailed_timer)
523 end subroutine check_small
529 subroutine ctoprim(self)
532 real(dp):: smallp, c, entho
533 real(dp):: eken, emag, etot, eint, b1, b2, b3
534 integer,
save:: itimer=0
535 integer:: i, j, k, loc(3), l(3), u(3)
536 associate(uin=>self%uin,q=>self%q,bf=>self%bf,dt=>self%dtime,c=>self%c)
538 call trace%begin(
'mhd_t%ctoprim', itimer=itimer, detailed_timer=detailed_timer)
542 entho = one/max(gamma-one,1e-6)
544 smallp = smallr*smallc**2/self%gamma
548 uin(:,:,:,1) = max(uin(:,:,:,1),smallr)
549 q(:,:,:,1) = uin(:,:,:,1)
550 q(:,:,:,2) = uin(:,:,:,2)/q(:,:,:,1)
551 q(:,:,:,3) = uin(:,:,:,3)/q(:,:,:,1)
552 q(:,:,:,4) = uin(:,:,:,4)/q(:,:,:,1)
554 q(:,:,:,6) = uin(:,:,:,6)
555 q(:,:,:,7) = uin(:,:,:,7)
556 q(:,:,:,8) = uin(:,:,:,8)
557 bf(:,:,:,1) = uin(:,:,:,6)
558 bf(:,:,:,2) = uin(:,:,:,7)
559 bf(:,:,:,3) = uin(:,:,:,8)
560 if (verbose>0) print *,
'ctoprim:rh', self%id, minval(q(:,:,:,1)), maxval(q(:,:,:,1))
561 if (verbose>0) print *,
'ctoprim:bf', self%id, minval(bf), maxval(bf)
566 eken = half*(q(i,j,k,2)**2 + q(i,j,k,3)**2 + q(i,j,k,4)**2)*q(i,j,k,1)
567 b1 = half*(q(i,j,k,6)+q(i+1,j,k,6))
568 b2 = half*(q(i,j,k,7)+q(i,j+1,k,7))
569 b3 = half*(q(i,j,k,8)+q(i,j,k+1,8))
570 emag = half*(b1**2 + b2**2 + b3**2)
571 q(i,j,k,5) = (uin(i,j,k,5) - emag - eken)*(gamma-1d0)
572 q(i,j,k,5) = merge(q(i,j,k,1)*self%csound**2, max(q(i,j,k,5),smallp), isothermal)
573 c(i,j,k) = sqrt((q(i,j,k,5)*gamma + 2.0*emag)/q(i,j,k,1))
580 if (courant_type==1)
then 581 c(:,:,:) = c(:,:,:) + max(abs(q(:,:,:,2)), abs(q(:,:,:,3)), abs(q(:,:,:,4)))
583 c(:,:,:) = c(:,:,:) + (abs(q(:,:,:,2)) + abs(q(:,:,:,3)) + abs(q(:,:,:,4)))/3.0
589 self%u_max = maxval(c(l(1):u(1),l(2):u(2),l(3):u(3)))
590 if (self%get_umax_location)
then 591 loc = maxloc(c(l(1):u(1),l(2):u(2),l(3):u(3))) + l-1
592 write (io_unit%log,
'(a,i6,f12.6,2(1x,3i4),3f9.4,2(1x,a,1p,4g10.2))')
'u_max:', &
593 self%id, self%time, self%ipos, loc, &
594 (self%mesh(1)%p + self%mesh(1)%r(loc(1))-self%mesh(2)%llc_cart)/self%mesh(1)%b, &
595 (self%mesh(2)%p + self%mesh(2)%r(loc(2))-self%mesh(2)%llc_cart)/self%mesh(2)%b, &
596 (self%mesh(3)%p + self%mesh(3)%r(loc(3))-self%mesh(2)%llc_cart)/self%mesh(3)%b, &
597 'P,B:', q(loc(1),loc(2),loc(3),5:8), &
598 'D,U:', q(loc(1),loc(2),loc(3),1:4)
602 call trace%end (itimer, detailed_timer=detailed_timer)
603 end subroutine ctoprim
608 subroutine uslope(self,bf,q,dq,dbf,rslope_type)
610 real,
dimension(:,:,:,:):: q
611 real,
dimension(:,:,:,:):: bf
612 real,
dimension(:,:,:,:,:):: dq
613 real,
dimension(:,:,:,:,:):: dbf
615 integer::nv, islope_type
616 integer::i, j, k, n, l(3), u(3)
617 real(dp)::dsgn, dlim, dcen, dlft, drgt, slop, xslope_type
618 real(dp)::dfll,dflm,dflr,dfml,dfmm,dfmr,dfrl,dfrm,dfrr
619 real(dp)::dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl
620 real(dp)::dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm
621 real(dp)::dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr
622 real(dp)::vmin,vmax,dfx,dfy,dfz,dff,idff
623 integer,
save:: itimer=0
625 call trace%begin (
'mhd_t%uslope', itimer=itimer, detailed_timer=detailed_timer)
627 islope_type = int(rslope_type)
631 if(slope_type==-1)
then 636 dq(i,j,k,n,1) = half*(q(i+1,j,k,n) - q(i-1,j,k,n))
637 dq(i,j,k,n,2) = half*(q(i,j+1,k,n) - q(i,j-1,k,n))
638 dq(i,j,k,n,3) = half*(q(i,j,k+1,n) - q(i,j,k-1,n))
643 else if(islope_type==0)
then 646 else if (islope_type==1)
then 652 dlft = q(i ,j,k,n) - q(i-1,j,k,n)
653 drgt = q(i+1,j,k,n) - q(i ,j,k,n)
654 dq(i,j,k,n,1) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
656 dlft = q(i,j ,k,n) - q(i,j-1,k,n)
657 drgt = q(i,j+1,k,n) - q(i,j ,k,n)
658 dq(i,j,k,n,2) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
660 dlft = q(i,j,k ,n) - q(i,j,k-1,n)
661 drgt = q(i,j,k+1,n) - q(i,j,k ,n)
662 dq(i,j,k,n,3) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
667 else if (islope_type==2)
then 674 dlft = (q(i ,j,k,n) - q(i-1,j,k,n))
675 drgt = (q(i+1,j,k,n) - q(i ,j,k,n))
676 dcen = half*(dlft+drgt)
677 dsgn = sign(one, dcen)
678 slop = xslope_type*min(abs(dlft),abs(drgt))
679 slop = merge(slop,zero,(dlft*drgt)>zero)
680 dq(i,j,k,n,1) = dsgn*min(slop,abs(dcen))
681 dlft = (q(i,j ,k,n) - q(i,j-1,k,n))
682 drgt = (q(i,j+1,k,n) - q(i,j ,k,n))
683 dcen = half*(dlft+drgt)
684 dsgn = sign(one, dcen)
685 slop = xslope_type*min(abs(dlft),abs(drgt))
686 slop = merge(slop,zero,(dlft*drgt)>zero)
687 dq(i,j,k,n,2) = dsgn*min(slop,abs(dcen))
689 dlft = (q(i,j,k ,n) - q(i,j,k-1,n))
690 drgt = (q(i,j,k+1,n) - q(i,j,k ,n))
691 dcen = half*(dlft+drgt)
692 dsgn = sign(one, dcen)
693 slop = xslope_type*min(abs(dlft),abs(drgt))
694 slop = merge(slop,zero,(dlft*drgt)>zero)
695 dq(i,j,k,n,3) = dsgn*min(slop,abs(dcen))
705 dlft = (bf(i,j ,k,1) - bf(i,j-1,k,1))
706 drgt = (bf(i,j+1,k,1) - bf(i,j ,k,1))
707 dcen = half*(dlft+drgt)
708 dsgn = sign(one, dcen)
709 slop = xslope_type*min(abs(dlft),abs(drgt))
710 slop = merge(slop,zero,(dlft*drgt)>zero)
711 dbf(i,j,k,1,1) = dsgn*min(slop,abs(dcen))
713 dlft = (bf(i,j,k ,1) - bf(i,j,k-1,1))
714 drgt = (bf(i,j,k+1,1) - bf(i,j,k ,1))
715 dcen = half*(dlft+drgt)
716 dsgn = sign(one, dcen)
717 slop = xslope_type*min(abs(dlft),abs(drgt))
718 slop = merge(slop,zero,(dlft*drgt)>zero)
719 dbf(i,j,k,1,2) = dsgn*min(slop,abs(dcen))
722 dlft = (bf(i ,j,k,2) - bf(i-1,j,k,2))
723 drgt = (bf(i+1,j,k,2) - bf(i ,j,k,2))
724 dcen = half*(dlft+drgt)
725 dsgn = sign(one, dcen)
726 slop = xslope_type*min(abs(dlft),abs(drgt))
727 slop = merge(slop,zero,(dlft*drgt)>zero)
728 dbf(i,j,k,2,1) = dsgn*min(slop,abs(dcen))
730 dlft = (bf(i,j,k ,2) - bf(i,j,k-1,2))
731 drgt = (bf(i,j,k+1,2) - bf(i,j,k ,2))
732 dcen = half*(dlft+drgt)
733 dsgn = sign(one, dcen)
734 slop = xslope_type*min(abs(dlft),abs(drgt))
735 slop = merge(slop,zero,(dlft*drgt)>zero)
736 dbf(i,j,k,2,2) = dsgn*min(slop,abs(dcen))
739 dlft = (bf(i ,j,k,3) - bf(i-1,j,k,3))
740 drgt = (bf(i+1,j,k,3) - bf(i ,j,k,3))
741 dcen = half*(dlft+drgt)
742 dsgn = sign(one, dcen)
743 slop = xslope_type*min(abs(dlft),abs(drgt))
744 slop = merge(slop,zero,(dlft*drgt)>zero)
745 dbf(i,j,k,3,1) = dsgn*min(slop,abs(dcen))
747 dlft = (bf(i,j ,k,3) - bf(i,j-1,k,3))
748 drgt = (bf(i,j+1,k,3) - bf(i,j ,k,3))
749 dcen = half*(dlft+drgt)
750 dsgn = sign(one, dcen)
751 slop = xslope_type*min(abs(dlft),abs(drgt))
752 slop = merge(slop,zero,(dlft*drgt)>zero)
753 dbf(i,j,k,3,2) = dsgn*min(slop,abs(dcen))
758 else if (islope_type==3)
then 759 xslope_type =
real(slope_type - 2_dp,kind=dp)
764 dflll = q(i-1,j-1,k-1,n)-q(i,j,k,n)
765 dflml = q(i-1,j ,k-1,n)-q(i,j,k,n)
766 dflrl = q(i-1,j+1,k-1,n)-q(i,j,k,n)
767 dfmll = q(i ,j-1,k-1,n)-q(i,j,k,n)
768 dfmml = q(i ,j ,k-1,n)-q(i,j,k,n)
769 dfmrl = q(i ,j+1,k-1,n)-q(i,j,k,n)
770 dfrll = q(i+1,j-1,k-1,n)-q(i,j,k,n)
771 dfrml = q(i+1,j ,k-1,n)-q(i,j,k,n)
772 dfrrl = q(i+1,j+1,k-1,n)-q(i,j,k,n)
774 dfllm = q(i-1,j-1,k ,n)-q(i,j,k,n)
775 dflmm = q(i-1,j ,k ,n)-q(i,j,k,n)
776 dflrm = q(i-1,j+1,k ,n)-q(i,j,k,n)
777 dfmlm = q(i ,j-1,k ,n)-q(i,j,k,n)
778 dfmmm = q(i ,j ,k ,n)-q(i,j,k,n)
779 dfmrm = q(i ,j+1,k ,n)-q(i,j,k,n)
780 dfrlm = q(i+1,j-1,k ,n)-q(i,j,k,n)
781 dfrmm = q(i+1,j ,k ,n)-q(i,j,k,n)
782 dfrrm = q(i+1,j+1,k ,n)-q(i,j,k,n)
784 dfllr = q(i-1,j-1,k+1,n)-q(i,j,k,n)
785 dflmr = q(i-1,j ,k+1,n)-q(i,j,k,n)
786 dflrr = q(i-1,j+1,k+1,n)-q(i,j,k,n)
787 dfmlr = q(i ,j-1,k+1,n)-q(i,j,k,n)
788 dfmmr = q(i ,j ,k+1,n)-q(i,j,k,n)
789 dfmrr = q(i ,j+1,k+1,n)-q(i,j,k,n)
790 dfrlr = q(i+1,j-1,k+1,n)-q(i,j,k,n)
791 dfrmr = q(i+1,j ,k+1,n)-q(i,j,k,n)
792 dfrrr = q(i+1,j+1,k+1,n)-q(i,j,k,n)
794 vmin = min(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
795 & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
796 & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
797 vmax = max(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
798 & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
799 & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
801 dfx = half*(q(i+1,j,k,n)-q(i-1,j,k,n))
802 dfy = half*(q(i,j+1,k,n)-q(i,j-1,k,n))
803 dfz = half*(q(i,j,k+1,n)-q(i,j,k-1,n))
804 idff = xslope_type / sqrt(dfx**2+dfy**2+dfz**2+tiny(1.0_dp))
815 dlim = merge(min(one,min(abs(vmin),abs(vmax))*idff), one, idff<0.001_dp * huge(1.0_dp))
817 dq(i,j,k,n,1) = dlim*dfx
818 dq(i,j,k,n,2) = dlim*dfy
819 dq(i,j,k,n,3) = dlim*dfz
828 dlft = (bf(i,j ,k,1) - bf(i,j-1,k,1))
829 drgt = (bf(i,j+1,k,1) - bf(i,j ,k,1))
830 dcen = half*(dlft+drgt)
831 dsgn = sign(one, dcen)
832 slop = xslope_type*min(abs(dlft),abs(drgt))
833 slop = merge(slop,zero,(dlft*drgt)>zero)
834 dbf(i,j,k,1,1) = dsgn*min(slop,abs(dcen))
836 dlft = (bf(i,j,k ,1) - bf(i,j,k-1,1))
837 drgt = (bf(i,j,k+1,1) - bf(i,j,k ,1))
838 dcen = half*(dlft+drgt)
839 dsgn = sign(one, dcen)
840 slop = xslope_type*min(abs(dlft),abs(drgt))
841 slop = merge(slop,zero,(dlft*drgt)>zero)
842 dbf(i,j,k,1,2) = dsgn*min(slop,abs(dcen))
845 dlft = (bf(i ,j,k,2) - bf(i-1,j,k,2))
846 drgt = (bf(i+1,j,k,2) - bf(i ,j,k,2))
847 dcen = half*(dlft+drgt)
848 dsgn = sign(one, dcen)
849 slop = xslope_type*min(abs(dlft),abs(drgt))
850 slop = merge(slop,zero,(dlft*drgt)>zero)
851 dbf(i,j,k,2,1) = dsgn*min(slop,abs(dcen))
853 dlft = (bf(i,j,k ,2) - bf(i,j,k-1,2))
854 drgt = (bf(i,j,k+1,2) - bf(i,j,k ,2))
855 dcen = half*(dlft+drgt)
856 dsgn = sign(one, dcen)
857 slop = xslope_type*min(abs(dlft),abs(drgt))
858 slop = merge(slop,zero,(dlft*drgt)>zero)
859 dbf(i,j,k,2,2) = dsgn*min(slop,abs(dcen))
862 dlft = (bf(i ,j,k,3) - bf(i-1,j,k,3))
863 drgt = (bf(i+1,j,k,3) - bf(i ,j,k,3))
864 dcen = half*(dlft+drgt)
865 dsgn = sign(one, dcen)
866 slop = xslope_type*min(abs(dlft),abs(drgt))
867 slop = merge(slop,zero,(dlft*drgt)>zero)
868 dbf(i,j,k,3,1) = dsgn*min(slop,abs(dcen))
870 dlft = (bf(i,j ,k,3) - bf(i,j-1,k,3))
871 drgt = (bf(i,j+1,k,3) - bf(i,j ,k,3))
872 dcen = half*(dlft+drgt)
873 dsgn = sign(one, dcen)
874 slop = xslope_type*min(abs(dlft),abs(drgt))
875 slop = merge(slop,zero,(dlft*drgt)>zero)
876 dbf(i,j,k,3,2) = dsgn*min(slop,abs(dcen))
882 write(*,*)
'Unknown slope type' 885 call trace%end (itimer, detailed_timer=detailed_timer)
886 end subroutine uslope
891 SUBROUTINE trace3d(self)
894 real(dp):: dx, dy, dz, dt
895 integer :: i, j, k, n
896 integer,
parameter ::ir=1, iu=2, iv=3, iw=4, ip=5, ia=6, ib=7, ic=8
897 real(dp)::dtdx, dtdy, dtdz, smallp
899 real(dp) :: r, u, v, w, p, a, b, c
900 real(dp) :: ell, elr, erl, err
901 real(dp) :: fll, flr, frl, frr
902 real(dp) :: gll, glr, grl, grr
903 real(dp) :: drx, dux, dvx, dwx, dpx, dax, dbx, dcx
904 real(dp) :: dry, duy, dvy, dwy, dpy, day, dby, dcy
905 real(dp) :: drz, duz, dvz, dwz, dpz, daz, dbz, dcz
906 real(dp) :: drx2,dpx2,dry2,dpy2,drz2,dpz2
907 real(dp) :: sr0, su0, sv0, sw0, sp0, sa0, sb0, sc0
908 real(dp) :: al, ar, bl, br, cl, cr
909 real(dp) :: daly, dary, dalz, darz
910 real(dp) :: dblx, dbrx, dblz, dbrz
911 real(dp) :: dclx, dcrx, dcly, dcry
912 real(dp) :: sal0, sar0, sbl0, sbr0, scl0, scr0
913 real(dp),
dimension(1:self%gn(1),1:self%gn(2),1:self%gn(3)):: rho_t, p_t
914 real(dp) :: eff_gamma, off
915 integer,
save:: itimer(4)=0
916 real:: ff(self%gn(1),3)
918 associate(q=>self%q,bf=>self%bf,dq=>self%dq,dq2=>self%dq2,qm=>self%qm,qp=>self%qp, &
919 dbf=>self%dbf,qrt=>self%qRT,qrb=>self%qRB,qlt=>self%qLT,qlb=>self%qLB, &
920 bx=>self%bf(:,:,:,1),by=>self%bf(:,:,:,2),bz=>self%bf(:,:,:,3), &
921 ux=>self%q(:,:,:,2),uy=>self%q(:,:,:,3),uz=>self%q(:,:,:,4), &
922 ex=>self%Ex,ey=>self%Ey,ez=>self%Ez)
924 call trace%begin (
'mhd_t%trace3d', itimer=itimer(1), detailed_timer=detailed_timer)
925 su0=0.0_dp; sv0=0.0_dp; sw0=0.0_dp
927 dtdx = self%dtime/self%ds(1)
928 dtdy = self%dtime/self%ds(2)
929 dtdz = self%dtime/self%ds(3)
930 smallp = smallr*smallc**2/gamma
933 ex = ydn(zdn(uy))*ydn(bz) - zdn(ydn(uz))*zdn(by)
934 ey = zdn(xdn(uz))*zdn(bx) - xdn(zdn(ux))*xdn(bz)
935 ez = xdn(ydn(ux))*xdn(by) - ydn(xdn(uy))*ydn(bx)
938 rho_t(:,:,:) = q(:,:,:,ir)
939 p_t(:,:,:) = q(:,:,:,ip)
949 if (
allocated(self%force_per_unit_mass))
then 950 ff(:,1) = self%force_per_unit_mass(:,j,k,1)
951 ff(:,2) = self%force_per_unit_mass(:,j,k,2)
952 ff(:,3) = self%force_per_unit_mass(:,j,k,3)
954 if (
allocated(self%force_per_unit_volume))
then 955 ff(:,1) = self%force_per_unit_mass(:,j,k,1)/q(:,j,k,ir)
956 ff(:,2) = self%force_per_unit_mass(:,j,k,2)/q(:,j,k,ir)
957 ff(:,3) = self%force_per_unit_mass(:,j,k,3)/q(:,j,k,ir)
992 drx = half * dq(i,j,k,ir,1)
993 drx2= half *dq2(i,j,k,ir,1)
994 dux = half * dq(i,j,k,iu,1)
995 dvx = half * dq(i,j,k,iv,1)
996 dwx = half * dq(i,j,k,iw,1)
997 dpx = half * dq(i,j,k,ip,1)
998 dpx2= half *dq2(i,j,k,ip,1)
999 dbx = half * dq(i,j,k,ib,1)
1000 dcx = half * dq(i,j,k,ic,1)
1002 dry = half * dq(i,j,k,ir,2)
1003 dry2= half *dq2(i,j,k,ir,2)
1004 duy = half * dq(i,j,k,iu,2)
1005 dvy = half * dq(i,j,k,iv,2)
1006 dwy = half * dq(i,j,k,iw,2)
1007 dpy = half * dq(i,j,k,ip,2)
1008 dpy2= half *dq2(i,j,k,ip,2)
1009 duy = half * dq(i,j,k,iu,2)
1010 day = half * dq(i,j,k,ia,2)
1011 dcy = half * dq(i,j,k,ic,2)
1013 drz = half * dq(i,j,k,ir,3)
1014 drz2= half *dq2(i,j,k,ir,3)
1015 duz = half * dq(i,j,k,iu,3)
1016 dvz = half * dq(i,j,k,iv,3)
1017 dwz = half * dq(i,j,k,iw,3)
1018 dpz = half * dq(i,j,k,ip,3)
1019 dpz2= half *dq2(i,j,k,ip,3)
1020 daz = half * dq(i,j,k,ia,3)
1021 dbz = half * dq(i,j,k,ib,3)
1024 daly = half * dbf(i ,j ,k ,1,1)
1025 dary = half * dbf(i+1,j ,k ,1,1)
1026 dalz = half * dbf(i ,j ,k ,1,2)
1027 darz = half * dbf(i+1,j ,k ,1,2)
1029 dblx = half * dbf(i ,j ,k ,2,1)
1030 dbrx = half * dbf(i ,j+1,k ,2,1)
1031 dblz = half * dbf(i ,j ,k ,2,2)
1032 dbrz = half * dbf(i ,j+1,k ,2,2)
1034 dclx = half * dbf(i ,j ,k ,3,1)
1035 dcrx = half * dbf(i ,j ,k+1,3,1)
1036 dcly = half * dbf(i ,j ,k ,3,2)
1037 dcry = half * dbf(i ,j ,k+1,3,2)
1041 sal0 = +(glr-gll)*dtdy*half -(flr-fll)*dtdz*half
1042 sar0 = +(grr-grl)*dtdy*half -(frr-frl)*dtdz*half
1043 sbl0 = -(grl-gll)*dtdx*half +(elr-ell)*dtdz*half
1044 sbr0 = -(grr-glr)*dtdx*half +(err-erl)*dtdz*half
1045 scl0 = +(frl-fll)*dtdx*half -(erl-ell)*dtdy*half
1046 scr0 = +(frr-flr)*dtdx*half -(err-elr)*dtdy*half
1057 sr0 = (-u*drx-dux*r)*dtdx + (-v*dry-dvy*r)*dtdy + (-w*drz-dwz*r)*dtdz
1058 su0 = (-u*dux-(dpx+b*dbx+c*dcx)*invr)*dtdx + (-v*duy+b*day*invr)*dtdy + (-w*duz+c*daz*invr)*dtdz
1059 sv0 = (-u*dvx+a*dbx*invr)*dtdx + (-v*dvy-(dpy+a*day+c*dcy)*invr)*dtdy + (-w*dvz+c*dbz*invr)*dtdz
1060 sw0 = (-u*dwx+a*dcx*invr)*dtdx + (-v*dwy+b*dcy*invr)*dtdy + (-w*dwz-(dpz+a*daz+b*dbz)*invr)*dtdz
1061 sp0 = (-u*dpx-dux*eff_gamma*p)*dtdx + (-v*dpy-dvy*eff_gamma*p)*dtdy + (-w*dpz-dwz*eff_gamma*p)*dtdz
1064 su0 = su0 + ff(i,1)*half*self%dtime
1065 sv0 = sv0 + ff(i,2)*half*self%dtime
1066 sw0 = sw0 + ff(i,3)*half*self%dtime
1113 qp(i,j,k,ir,1) = r - drx
1114 qp(i,j,k,iu,1) = u - dux
1115 qp(i,j,k,iv,1) = v - dvx
1116 qp(i,j,k,iw,1) = w - dwx
1117 qp(i,j,k,ip,1) = p - dpx
1119 qp(i,j,k,ib,1) = b - dbx
1120 qp(i,j,k,ic,1) = c - dcx
1123 qm(i,j,k,ir,1) = r + drx
1124 qm(i,j,k,iu,1) = u + dux
1125 qm(i,j,k,iv,1) = v + dvx
1126 qm(i,j,k,iw,1) = w + dwx
1127 qm(i,j,k,ip,1) = p + dpx
1129 qm(i,j,k,ib,1) = b + dbx
1130 qm(i,j,k,ic,1) = c + dcx
1133 qp(i,j,k,ir,2) = r - dry
1134 qp(i,j,k,iu,2) = u - duy
1135 qp(i,j,k,iv,2) = v - dvy
1136 qp(i,j,k,iw,2) = w - dwy
1137 qp(i,j,k,ip,2) = p - dpy
1138 qp(i,j,k,ia,2) = a - day
1140 qp(i,j,k,ic,2) = c - dcy
1143 qm(i,j,k,ir,2) = r + dry
1144 qm(i,j,k,iu,2) = u + duy
1145 qm(i,j,k,iv,2) = v + dvy
1146 qm(i,j,k,iw,2) = w + dwy
1147 qm(i,j,k,ip,2) = p + dpy
1148 qm(i,j,k,ia,2) = a + day
1150 qm(i,j,k,ic,2) = c + dcy
1153 qp(i,j,k,ir,3) = r - drz
1154 qp(i,j,k,iu,3) = u - duz
1155 qp(i,j,k,iv,3) = v - dvz
1156 qp(i,j,k,iw,3) = w - dwz
1157 qp(i,j,k,ip,3) = p - dpz
1158 qp(i,j,k,ia,3) = a - daz
1159 qp(i,j,k,ib,3) = b - dbz
1163 qm(i,j,k,ir,3) = r + drz
1164 qm(i,j,k,iu,3) = u + duz
1165 qm(i,j,k,iv,3) = v + dvz
1166 qm(i,j,k,iw,3) = w + dwz
1167 qm(i,j,k,ip,3) = p + dpz
1168 qm(i,j,k,ia,3) = a + daz
1169 qm(i,j,k,ib,3) = b + dbz
1173 qrt(i,j,k,ir,1) = r + (+dry2+drz2)
1174 qrt(i,j,k,iu,1) = u + (+duy+duz)
1175 qrt(i,j,k,iv,1) = v + (+dvy+dvz)
1176 qrt(i,j,k,iw,1) = w + (+dwy+dwz)
1177 qrt(i,j,k,ip,1) = p + (+dpy2+dpz2)
1178 qrt(i,j,k,ia,1) = a + (+day+daz)
1179 qrt(i,j,k,ib,1) = br+ ( +dbrz)
1180 qrt(i,j,k,ic,1) = cr+ (+dcry )
1183 qrb(i,j,k,ir,1) = r + (+dry2-drz2)
1184 qrb(i,j,k,iu,1) = u + (+duy-duz)
1185 qrb(i,j,k,iv,1) = v + (+dvy-dvz)
1186 qrb(i,j,k,iw,1) = w + (+dwy-dwz)
1187 qrb(i,j,k,ip,1) = p + (+dpy2-dpz2)
1188 qrb(i,j,k,ia,1) = a + (+day-daz)
1189 qrb(i,j,k,ib,1) = br+ ( -dbrz)
1190 qrb(i,j,k,ic,1) = cl+ (+dcly )
1193 qlt(i,j,k,ir,1) = r + (-dry2+drz2)
1194 qlt(i,j,k,iu,1) = u + (-duy+duz)
1195 qlt(i,j,k,iv,1) = v + (-dvy+dvz)
1196 qlt(i,j,k,iw,1) = w + (-dwy+dwz)
1197 qlt(i,j,k,ip,1) = p + (-dpy2+dpz2)
1198 qlt(i,j,k,ia,1) = a + (-day+daz)
1199 qlt(i,j,k,ib,1) = bl+ ( +dblz)
1200 qlt(i,j,k,ic,1) = cr+ (-dcry )
1203 qlb(i,j,k,ir,1) = r + (-dry2-drz2)
1204 qlb(i,j,k,iu,1) = u + (-duy-duz)
1205 qlb(i,j,k,iv,1) = v + (-dvy-dvz)
1206 qlb(i,j,k,iw,1) = w + (-dwy-dwz)
1207 qlb(i,j,k,ip,1) = p + (-dpy2-dpz2)
1208 qlb(i,j,k,ia,1) = a + (-day-daz)
1209 qlb(i,j,k,ib,1) = bl+ ( -dblz)
1210 qlb(i,j,k,ic,1) = cl+ (-dcly )
1213 qrt(i,j,k,ir,2) = r + (+drx2+drz2)
1214 qrt(i,j,k,iu,2) = u + (+dux+duz)
1215 qrt(i,j,k,iv,2) = v + (+dvx+dvz)
1216 qrt(i,j,k,iw,2) = w + (+dwx+dwz)
1217 qrt(i,j,k,ip,2) = p + (+dpx2+dpz2)
1218 qrt(i,j,k,ia,2) = ar+ ( +darz)
1219 qrt(i,j,k,ib,2) = b + (+dbx+dbz)
1220 qrt(i,j,k,ic,2) = cr+ (+dcrx )
1223 qrb(i,j,k,ir,2) = r + (+drx2-drz2)
1224 qrb(i,j,k,iu,2) = u + (+dux-duz)
1225 qrb(i,j,k,iv,2) = v + (+dvx-dvz)
1226 qrb(i,j,k,iw,2) = w + (+dwx-dwz)
1227 qrb(i,j,k,ip,2) = p + (+dpx2-dpz2)
1228 qrb(i,j,k,ia,2) = ar+ ( -darz)
1229 qrb(i,j,k,ib,2) = b + (+dbx-dbz)
1230 qrb(i,j,k,ic,2) = cl+ (+dclx )
1233 qlt(i,j,k,ir,2) = r + (-drx2+drz2)
1234 qlt(i,j,k,iu,2) = u + (-dux+duz)
1235 qlt(i,j,k,iv,2) = v + (-dvx+dvz)
1236 qlt(i,j,k,iw,2) = w + (-dwx+dwz)
1237 qlt(i,j,k,ip,2) = p + (-dpx2+dpz2)
1238 qlt(i,j,k,ia,2) = al+ ( +dalz)
1239 qlt(i,j,k,ib,2) = b + (-dbx+dbz)
1240 qlt(i,j,k,ic,2) = cr+ (-dcrx )
1243 qlb(i,j,k,ir,2) = r + (-drx2-drz2)
1244 qlb(i,j,k,iu,2) = u + (-dux-duz)
1245 qlb(i,j,k,iv,2) = v + (-dvx-dvz)
1246 qlb(i,j,k,iw,2) = w + (-dwx-dwz)
1247 qlb(i,j,k,ip,2) = p + (-dpx2-dpz2)
1248 qlb(i,j,k,ia,2) = al+ ( -dalz)
1249 qlb(i,j,k,ib,2) = b + (-dbx-dbz)
1250 qlb(i,j,k,ic,2) = cl+ (-dclx )
1253 qrt(i,j,k,ir,3) = r + (+drx2+dry2)
1254 qrt(i,j,k,iu,3) = u + (+dux+duy)
1255 qrt(i,j,k,iv,3) = v + (+dvx+dvy)
1256 qrt(i,j,k,iw,3) = w + (+dwx+dwy)
1257 qrt(i,j,k,ip,3) = p + (+dpx2+dpy2)
1258 qrt(i,j,k,ia,3) = ar+ ( +dary)
1259 qrt(i,j,k,ib,3) = br+ (+dbrx )
1260 qrt(i,j,k,ic,3) = c + (+dcx+dcy)
1263 qrb(i,j,k,ir,3) = r + (+drx2-dry2)
1264 qrb(i,j,k,iu,3) = u + (+dux-duy)
1265 qrb(i,j,k,iv,3) = v + (+dvx-dvy)
1266 qrb(i,j,k,iw,3) = w + (+dwx-dwy)
1267 qrb(i,j,k,ip,3) = p + (+dpx2-dpy2)
1268 qrb(i,j,k,ia,3) = ar+ ( -dary)
1269 qrb(i,j,k,ib,3) = bl+ (+dblx )
1270 qrb(i,j,k,ic,3) = c + (+dcx-dcy)
1273 qlt(i,j,k,ir,3) = r + (-drx2+dry2)
1274 qlt(i,j,k,iu,3) = u + (-dux+duy)
1275 qlt(i,j,k,iv,3) = v + (-dvx+dvy)
1276 qlt(i,j,k,iw,3) = w + (-dwx+dwy)
1277 qlt(i,j,k,ip,3) = p + (-dpx2+dpy2)
1278 qlt(i,j,k,ia,3) = al+ ( +daly)
1279 qlt(i,j,k,ib,3) = br+ (-dbrx )
1280 qlt(i,j,k,ic,3) = c + (-dcx+dcy)
1283 qlb(i,j,k,ir,3) = r + (-drx2-dry2)
1284 qlb(i,j,k,iu,3) = u + (-dux-duy)
1285 qlb(i,j,k,iv,3) = v + (-dvx-dvy)
1286 qlb(i,j,k,iw,3) = w + (-dwx-dwy)
1287 qlb(i,j,k,ip,3) = p + (-dpx2-dpy2)
1288 qlb(i,j,k,ia,3) = al+ ( -daly)
1289 qlb(i,j,k,ib,3) = bl+ (-dblx )
1290 qlb(i,j,k,ic,3) = c + (-dcx-dcy)
1305 drx = sqrt(sqrt(r*rho_t(i,j+1,k)*rho_t(i,j,k+1)*rho_t(i,j+1,k+1))) - r
1306 dpx = sqrt(sqrt(p* p_t(i,j+1,k)* p_t(i,j,k+1)* p_t(i,j+1,k+1))) - p
1307 dry = qrt(i,j,k,ir,1) - r
1308 dpy = qrt(i,j,k,ip,1) - p
1309 qrt(i,j,k,ir,1) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1310 qrt(i,j,k,ip,1) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1313 drx = sqrt(sqrt(r*rho_t(i,j+1,k)*rho_t(i,j,k-1)*rho_t(i,j+1,k-1))) - r
1314 dpx = sqrt(sqrt(p* p_t(i,j+1,k)* p_t(i,j,k-1)* p_t(i,j+1,k-1))) - p
1315 dry = qrb(i,j,k,ir,1) - r
1316 dpy = qrb(i,j,k,ip,1) - p
1317 qrb(i,j,k,ir,1) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1318 qrb(i,j,k,ip,1) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1321 drx = sqrt(sqrt(r*rho_t(i,j-1,k)*rho_t(i,j,k+1)*rho_t(i,j-1,k+1))) - r
1322 dpx = sqrt(sqrt(p* p_t(i,j-1,k)* p_t(i,j,k+1)* p_t(i,j-1,k+1))) - p
1323 dry = qlt(i,j,k,ir,1) - r
1324 dpy = qlt(i,j,k,ip,1) - p
1325 qlt(i,j,k,ir,1) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1326 qlt(i,j,k,ip,1) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1329 drx = sqrt(sqrt(r*rho_t(i,j-1,k)*rho_t(i,j,k-1)*rho_t(i,j-1,k-1))) - r
1330 dpx = sqrt(sqrt(p* p_t(i,j-1,k)* p_t(i,j,k-1)* p_t(i,j-1,k-1))) - p
1331 dry = qlb(i,j,k,ir,1) - r
1332 dpy = qlb(i,j,k,ip,1) - p
1333 qlb(i,j,k,ir,1) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1334 qlb(i,j,k,ip,1) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1337 drx = sqrt(sqrt(r*rho_t(i+1,j,k)*rho_t(i,j,k+1)*rho_t(i+1,j,k+1))) - r
1338 dpx = sqrt(sqrt(p* p_t(i+1,j,k)* p_t(i,j,k+1)* p_t(i+1,j,k+1))) - p
1339 dry = qrt(i,j,k,ir,2) - r
1340 dpy = qrt(i,j,k,ip,2) - p
1341 qrt(i,j,k,ir,2) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1342 qrt(i,j,k,ip,2) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1345 drx = sqrt(sqrt(r*rho_t(i+1,j,k)*rho_t(i,j,k-1)*rho_t(i+1,j,k-1))) - r
1346 dpx = sqrt(sqrt(p* p_t(i+1,j,k)* p_t(i,j,k-1)* p_t(i+1,j,k-1))) - p
1347 dry = qrb(i,j,k,ir,2) - r
1348 dpy = qrb(i,j,k,ip,2) - p
1349 qrb(i,j,k,ir,2) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1350 qrb(i,j,k,ip,2) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1354 drx = sqrt(sqrt(r*rho_t(i-1,j,k)*rho_t(i,j,k+1)*rho_t(i-1,j,k+1))) - r
1355 dpx = sqrt(sqrt(p* p_t(i-1,j,k)* p_t(i,j,k+1)* p_t(i-1,j,k+1))) - p
1356 dry = qlt(i,j,k,ir,2) - r
1357 dpy = qlt(i,j,k,ip,2) - p
1358 qlt(i,j,k,ir,2) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1359 qlt(i,j,k,ip,2) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1362 drx = sqrt(sqrt(r*rho_t(i-1,j,k)*rho_t(i,j,k-1)*rho_t(i-1,j,k-1))) - r
1363 dpx = sqrt(sqrt(p* p_t(i-1,j,k)* p_t(i,j,k-1)* p_t(i-1,j,k-1))) - p
1364 dry = qlb(i,j,k,ir,2) - r
1365 dpy = qlb(i,j,k,ip,2) - p
1366 qlb(i,j,k,ir,2) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1367 qlb(i,j,k,ip,2) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1370 drx = sqrt(sqrt(r*rho_t(i+1,j,k)*rho_t(i,j+1,k)*rho_t(i+1,j+1,k))) - r
1371 dpx = sqrt(sqrt(p* p_t(i+1,j,k)* p_t(i,j+1,k)* p_t(i+1,j+1,k))) - p
1372 dry = qrt(i,j,k,ir,3) - r
1373 dpy = qrt(i,j,k,ip,3) - p
1374 qrt(i,j,k,ir,3) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1375 qrt(i,j,k,ip,3) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1378 drx = sqrt(sqrt(r*rho_t(i+1,j,k)*rho_t(i,j-1,k)*rho_t(i+1,j-1,k))) - r
1379 dpx = sqrt(sqrt(p* p_t(i+1,j,k)* p_t(i,j-1,k)* p_t(i+1,j-1,k))) - p
1380 dry = qrb(i,j,k,ir,3) - r
1381 dpy = qrb(i,j,k,ip,3) - p
1382 qrb(i,j,k,ir,3) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1383 qrb(i,j,k,ip,3) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1386 drx = sqrt(sqrt(r*rho_t(i-1,j,k)*rho_t(i,j+1,k)*rho_t(i-1,j+1,k))) - r
1387 dpx = sqrt(sqrt(p* p_t(i-1,j,k)* p_t(i,j+1,k)* p_t(i-1,j+1,k))) - p
1388 dry = qlt(i,j,k,ir,3) - r
1389 dpy = qlt(i,j,k,ip,3) - p
1390 qlt(i,j,k,ir,3) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1391 qlt(i,j,k,ip,3) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1394 drx = sqrt(sqrt(r*rho_t(i-1,j,k)*rho_t(i,j-1,k)*rho_t(i-1,j-1,k))) - r
1395 dpx = sqrt(sqrt(p* p_t(i-1,j,k)* p_t(i,j-1,k)* p_t(i-1,j-1,k))) - p
1396 dry = qlb(i,j,k,ir,3) - r
1397 dpy = qlb(i,j,k,ip,3) - p
1398 qlb(i,j,k,ir,3) = merge(drx, dry, abs(drx) <= abs(dry) .or. dry+r<=0.) + r
1399 qlb(i,j,k,ip,3) = merge(dpx, dpy, abs(dpx) <= abs(dpy) .or. dpy+p<=0.) + p
1406 call trace%end (itimer(1), detailed_timer=detailed_timer)
1409 FUNCTION xdn(a)
RESULT (b)
1410 real,
dimension(:,:,:):: a
1411 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
1415 b(i,j,k)=half*(a(i-1,j,k)+a(i,j,k))
1421 FUNCTION ydn(a)
RESULT (b)
1422 real,
dimension(:,:,:):: a
1423 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
1427 b(i,j,k)=half*(a(i,j-1,k)+a(i,j,k))
1433 FUNCTION zdn(a)
RESULT (b)
1434 real,
dimension(:,:,:):: a
1435 real,
dimension(size(a,1),size(a,2),size(a,3)):: b
1439 b(i,j,k)=half*(a(i,j,k-1)+a(i,j,k))
1445 END SUBROUTINE trace3d
1450 subroutine cmpflxm(self, nx, ln, lt1, lt2, bn, bt1, bt2, fl)
1452 real(dp),
dimension(:,:,:,:):: fl
1453 integer:: nx, ln, lt1, lt2, bn, bt1, bt2
1455 integer:: i, j, k, lb(3), ub(3), xdim, l(3), idim
1456 real(dp),
dimension(nx,self%nv):: qleft, qright, fgdnv
1457 real(dp),
dimension(nx):: bn_mean
1459 integer,
save:: itimer=0
1461 associate(qp=>self%qp, qm=>self%qm)
1463 call trace%begin (
'mhd_t%cmpflxm', itimer=itimer, detailed_timer=detailed_timer)
1467 if (verbose>0) print
'(a,2(2x,3i4))',
'mhd_t%cmpflxm: lb, ub=', lb, ub
1477 bn_mean(1:n) = half*(qm(lb(1)-1:ub(1)-1,j,k,bn,xdim)+qp(lb(1):ub(1),j,k,bn,xdim))
1479 qleft(1:n,1) = qm(lb(1)-1:ub(1)-1,j,k,1 ,xdim)
1480 qleft(1:n,2) = qm(lb(1)-1:ub(1)-1,j,k,5 ,xdim)
1481 qleft(1:n,3) = qm(lb(1)-1:ub(1)-1,j,k,ln ,xdim)
1482 qleft(1:n,4) = bn_mean(1:n)
1483 qleft(1:n,5) = qm(lb(1)-1:ub(1)-1,j,k,lt1,xdim)
1484 qleft(1:n,6) = qm(lb(1)-1:ub(1)-1,j,k,bt1,xdim)
1485 qleft(1:n,7) = qm(lb(1)-1:ub(1)-1,j,k,lt2,xdim)
1486 qleft(1:n,8) = qm(lb(1)-1:ub(1)-1,j,k,bt2,xdim)
1488 qright(1:n,1) = qp(lb(1):ub(1),j,k,1 ,xdim)
1489 qright(1:n,2) = qp(lb(1):ub(1),j,k,5 ,xdim)
1490 qright(1:n,3) = qp(lb(1):ub(1),j,k,ln ,xdim)
1491 qright(1:n,4) = bn_mean(1:n)
1492 qright(1:n,5) = qp(lb(1):ub(1),j,k,lt1,xdim)
1493 qright(1:n,6) = qp(lb(1):ub(1),j,k,bt1,xdim)
1494 qright(1:n,7) = qp(lb(1):ub(1),j,k,lt2,xdim)
1495 qright(1:n,8) = qp(lb(1):ub(1),j,k,bt2,xdim)
1496 call solver (self,qleft,qright,fgdnv,n)
1498 fl(lb(1):ub(1),j,k,1 ) = fgdnv(1:n,1)
1499 fl(lb(1):ub(1),j,k,5 ) = fgdnv(1:n,2)
1500 fl(lb(1):ub(1),j,k,ln ) = fgdnv(1:n,3)
1501 fl(lb(1):ub(1),j,k,bn ) = fgdnv(1:n,4)
1502 fl(lb(1):ub(1),j,k,lt1) = fgdnv(1:n,5)
1503 fl(lb(1):ub(1),j,k,bt1) = fgdnv(1:n,6)
1504 fl(lb(1):ub(1),j,k,lt2) = fgdnv(1:n,7)
1505 fl(lb(1):ub(1),j,k,bt2) = fgdnv(1:n,8)
1508 else if (xdim == 2)
then 1513 bn_mean(1:n) = half*(qm(i,lb(2)-1:ub(2)-1,k,bn,xdim)+qp(i,lb(2):ub(2),k,bn,xdim))
1515 qleft(1:n,1) = qm(i,lb(2)-1:ub(2)-1,k,1 ,xdim)
1516 qleft(1:n,2) = qm(i,lb(2)-1:ub(2)-1,k,5 ,xdim)
1517 qleft(1:n,3) = qm(i,lb(2)-1:ub(2)-1,k,ln ,xdim)
1518 qleft(1:n,4) = bn_mean(1:n)
1519 qleft(1:n,5) = qm(i,lb(2)-1:ub(2)-1,k,lt1,xdim)
1520 qleft(1:n,6) = qm(i,lb(2)-1:ub(2)-1,k,bt1,xdim)
1521 qleft(1:n,7) = qm(i,lb(2)-1:ub(2)-1,k,lt2,xdim)
1522 qleft(1:n,8) = qm(i,lb(2)-1:ub(2)-1,k,bt2,xdim)
1524 qright(1:n,1) = qp(i,lb(2):ub(2),k,1 ,xdim)
1525 qright(1:n,2) = qp(i,lb(2):ub(2),k,5 ,xdim)
1526 qright(1:n,3) = qp(i,lb(2):ub(2),k,ln ,xdim)
1527 qright(1:n,4) = bn_mean(1:n)
1528 qright(1:n,5) = qp(i,lb(2):ub(2),k,lt1,xdim)
1529 qright(1:n,6) = qp(i,lb(2):ub(2),k,bt1,xdim)
1530 qright(1:n,7) = qp(i,lb(2):ub(2),k,lt2,xdim)
1531 qright(1:n,8) = qp(i,lb(2):ub(2),k,bt2,xdim)
1532 call solver (self,qleft,qright,fgdnv,n)
1534 fl(i,lb(2):ub(2),k,1 ) = fgdnv(1:n,1)
1535 fl(i,lb(2):ub(2),k,5 ) = fgdnv(1:n,2)
1536 fl(i,lb(2):ub(2),k,ln ) = fgdnv(1:n,3)
1537 fl(i,lb(2):ub(2),k,bn ) = fgdnv(1:n,4)
1538 fl(i,lb(2):ub(2),k,lt1) = fgdnv(1:n,5)
1539 fl(i,lb(2):ub(2),k,bt1) = fgdnv(1:n,6)
1540 fl(i,lb(2):ub(2),k,lt2) = fgdnv(1:n,7)
1541 fl(i,lb(2):ub(2),k,bt2) = fgdnv(1:n,8)
1544 else if (xdim == 3)
then 1549 bn_mean(1:n) = half*(qm(i,j,lb(3)-1:ub(3)-1,bn,xdim)+qp(i,j,lb(3):ub(3),bn,xdim))
1551 qleft(1:n,1) = qm(i,j,lb(3)-1:ub(3)-1,1 ,xdim)
1552 qleft(1:n,2) = qm(i,j,lb(3)-1:ub(3)-1,5 ,xdim)
1553 qleft(1:n,3) = qm(i,j,lb(3)-1:ub(3)-1,ln ,xdim)
1554 qleft(1:n,4) = bn_mean(1:n)
1555 qleft(1:n,5) = qm(i,j,lb(3)-1:ub(3)-1,lt1,xdim)
1556 qleft(1:n,6) = qm(i,j,lb(3)-1:ub(3)-1,bt1,xdim)
1557 qleft(1:n,7) = qm(i,j,lb(3)-1:ub(3)-1,lt2,xdim)
1558 qleft(1:n,8) = qm(i,j,lb(3)-1:ub(3)-1,bt2,xdim)
1560 qright(1:n,1) = qp(i,j,lb(3):ub(3),1 ,xdim)
1561 qright(1:n,2) = qp(i,j,lb(3):ub(3),5 ,xdim)
1562 qright(1:n,3) = qp(i,j,lb(3):ub(3),ln ,xdim)
1563 qright(1:n,4) = bn_mean(1:n)
1564 qright(1:n,5) = qp(i,j,lb(3):ub(3),lt1,xdim)
1565 qright(1:n,6) = qp(i,j,lb(3):ub(3),bt1,xdim)
1566 qright(1:n,7) = qp(i,j,lb(3):ub(3),lt2,xdim)
1567 qright(1:n,8) = qp(i,j,lb(3):ub(3),bt2,xdim)
1568 call solver (self,qleft,qright,fgdnv,n)
1570 fl(i,j,lb(3):ub(3),1 ) = fgdnv(1:n,1)
1571 fl(i,j,lb(3):ub(3),5 ) = fgdnv(1:n,2)
1572 fl(i,j,lb(3):ub(3),ln ) = fgdnv(1:n,3)
1573 fl(i,j,lb(3):ub(3),bn ) = fgdnv(1:n,4)
1574 fl(i,j,lb(3):ub(3),lt1) = fgdnv(1:n,5)
1575 fl(i,j,lb(3):ub(3),bt1) = fgdnv(1:n,6)
1576 fl(i,j,lb(3):ub(3),lt2) = fgdnv(1:n,7)
1577 fl(i,j,lb(3):ub(3),bt2) = fgdnv(1:n,8)
1584 call trace%end (itimer, detailed_timer=detailed_timer)
1585 end subroutine cmpflxm
1590 SUBROUTINE cmp_mag_flx (self, qRT, qRB, qLT, qLB, lp1 ,lp2 ,lor ,bp1 , bp2, bor, emf, ngrid, nvar)
1595 integer:: lp1, lp2, lor, bp1, bp2, bor, ngrid, nvar
1596 integer:: irt, jrt, krt, irb, jrb, krb, ilt, jlt, klt, ilb, jlb, klb
1597 real(dp),
DIMENSION(:,:,:):: emf
1598 real(dp),
DIMENSION(:,:,:,:,:):: qrt, qlt, qrb, qlb
1600 integer:: i, j, k, n, xdim, l, u
1601 real(dp),
DIMENSION(1:ngrid,nvar)::qll,qrl,qlr,qrr
1603 integer,
save:: itimer=0
1605 call trace%begin (
'mhd_t%cmp_mag_flx', itimer=itimer, detailed_timer=detailed_timer)
1614 smallp = smallr*smallc**2
1622 irt = 1; jrt = 1; krt = 0
1623 irb = 1; jrb = 0; krb = 0
1624 ilt = 0; jlt = 1; klt = 0
1625 ilb = 0; jlb = 0; klb = 0
1626 else if (xdim == 2)
then 1627 irt = 1; jrt = 0; krt = 1
1628 irb = 0; jrb = 0; krb = 1
1629 ilt = 1; jlt = 0; klt = 0
1630 ilb = 0; jlb = 0; klb = 0
1631 else if (xdim == 1)
then 1632 irt = 0; jrt = 1; krt = 1
1633 irb = 0; jrb = 1; krb = 0
1634 ilt = 0; jlt = 0; klt = 1
1635 ilb = 0; jlb = 0; klb = 0
1639 do k = self%mesh(3)%li,self%mesh(3)%uo
1640 do j = self%mesh(2)%li,self%mesh(2)%uo
1644 qll(i,1) = max(qrt(i-irt, j-jrt, k-krt, 1, xdim), smallr)
1645 qrl(i,1) = max(qlt(i-ilt, j-jlt, k-klt, 1, xdim), smallr)
1646 qlr(i,1) = max(qrb(i-irb, j-jrb, k-krb, 1, xdim), smallr)
1647 qrr(i,1) = max(qlb(i-ilb, j-jlb, k-klb, 1, xdim), smallr)
1648 qll(i,2) = max(qrt(i-irt, j-jrt, k-krt, 5, xdim), smallp)
1649 qrl(i,2) = max(qlt(i-ilt, j-jlt, k-klt, 5, xdim), smallp)
1650 qlr(i,2) = max(qrb(i-irb, j-jrb, k-krb, 5, xdim), smallp)
1651 qrr(i,2) = max(qlb(i-ilb, j-jlb, k-klb, 5, xdim), smallp)
1652 qll(i,3) = qrt(i-irt, j-jrt, k-krt, lp1, xdim)
1653 qrl(i,3) = qlt(i-ilt, j-jlt, k-klt, lp1, xdim)
1654 qlr(i,3) = qrb(i-irb, j-jrb, k-krb, lp1, xdim)
1655 qrr(i,3) = qlb(i-ilb, j-jlb, k-klb, lp1, xdim)
1656 qll(i,4) = qrt(i-irt, j-jrt, k-krt, lp2, xdim)
1657 qrl(i,4) = qlt(i-ilt, j-jlt, k-klt, lp2, xdim)
1658 qlr(i,4) = qrb(i-irb, j-jrb, k-krb, lp2, xdim)
1659 qrr(i,4) = qlb(i-ilb, j-jlb, k-klb, lp2, xdim)
1660 qll(i,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1661 qrl(i,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1662 qlr(i,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1663 qrr(i,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1664 qll(i,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1665 qrl(i,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1666 qlr(i,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1667 qrr(i,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1668 qll(i,5) = qrt(i-irt, j-jrt, k-krt, lor, xdim)
1669 qrl(i,5) = qlt(i-ilt, j-jlt, k-klt, lor, xdim)
1670 qlr(i,5) = qrb(i-irb, j-jrb, k-krb, lor, xdim)
1671 qrr(i,5) = qlb(i-ilb, j-jlb, k-klb, lor, xdim)
1672 qll(i,8) = qrt(i-irt, j-jrt, k-krt, bor, xdim)
1673 qrl(i,8) = qlt(i-ilt, j-jlt, k-klt, bor, xdim)
1674 qlr(i,8) = qrb(i-irb, j-jrb, k-krb, bor, xdim)
1675 qrr(i,8) = qlb(i-ilb, j-jlb, k-klb, bor, xdim)
1677 call solver_2d (self,qll,qrl,qlr,qrr,emf,j,k,l,u,nvar,xdim,detailed_timer=detailed_timer)
1680 else if (xdim == 2)
then 1681 do i = self%mesh(1)%li,self%mesh(1)%uo
1682 do k = self%mesh(3)%li,self%mesh(3)%uo
1686 qll(j,1) = max(qrt(i-irt, j-jrt, k-krt, 1, xdim), smallr)
1687 qrl(j,1) = max(qlt(i-ilt, j-jlt, k-klt, 1, xdim), smallr)
1688 qlr(j,1) = max(qrb(i-irb, j-jrb, k-krb, 1, xdim), smallr)
1689 qrr(j,1) = max(qlb(i-ilb, j-jlb, k-klb, 1, xdim), smallr)
1690 qll(j,2) = max(qrt(i-irt, j-jrt, k-krt, 5, xdim), smallp)
1691 qrl(j,2) = max(qlt(i-ilt, j-jlt, k-klt, 5, xdim), smallp)
1692 qlr(j,2) = max(qrb(i-irb, j-jrb, k-krb, 5, xdim), smallp)
1693 qrr(j,2) = max(qlb(i-ilb, j-jlb, k-klb, 5, xdim), smallp)
1694 qll(j,3) = qrt(i-irt, j-jrt, k-krt, lp1, xdim)
1695 qrl(j,3) = qlt(i-ilt, j-jlt, k-klt, lp1, xdim)
1696 qlr(j,3) = qrb(i-irb, j-jrb, k-krb, lp1, xdim)
1697 qrr(j,3) = qlb(i-ilb, j-jlb, k-klb, lp1, xdim)
1698 qll(j,4) = qrt(i-irt, j-jrt, k-krt, lp2, xdim)
1699 qrl(j,4) = qlt(i-ilt, j-jlt, k-klt, lp2, xdim)
1700 qlr(j,4) = qrb(i-irb, j-jrb, k-krb, lp2, xdim)
1701 qrr(j,4) = qlb(i-ilb, j-jlb, k-klb, lp2, xdim)
1702 qll(j,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1703 qrl(j,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1704 qlr(j,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1705 qrr(j,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1706 qll(j,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1707 qrl(j,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1708 qlr(j,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1709 qrr(j,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1710 qll(j,5) = qrt(i-irt, j-jrt, k-krt, lor, xdim)
1711 qrl(j,5) = qlt(i-ilt, j-jlt, k-klt, lor, xdim)
1712 qlr(j,5) = qrb(i-irb, j-jrb, k-krb, lor, xdim)
1713 qrr(j,5) = qlb(i-ilb, j-jlb, k-klb, lor, xdim)
1714 qll(j,8) = qrt(i-irt, j-jrt, k-krt, bor, xdim)
1715 qrl(j,8) = qlt(i-ilt, j-jlt, k-klt, bor, xdim)
1716 qlr(j,8) = qrb(i-irb, j-jrb, k-krb, bor, xdim)
1717 qrr(j,8) = qlb(i-ilb, j-jlb, k-klb, bor, xdim)
1719 call solver_2d (self,qll,qrl,qlr,qrr,emf,k,i,l,u,nvar,xdim,detailed_timer=detailed_timer)
1722 else if (xdim == 3)
then 1723 do j = self%mesh(2)%li,self%mesh(2)%uo
1724 do i = self%mesh(1)%li,self%mesh(1)%uo
1728 qll(k,1) = max(qrt(i-irt, j-jrt, k-krt, 1, xdim), smallr)
1729 qrl(k,1) = max(qlt(i-ilt, j-jlt, k-klt, 1, xdim), smallr)
1730 qlr(k,1) = max(qrb(i-irb, j-jrb, k-krb, 1, xdim), smallr)
1731 qrr(k,1) = max(qlb(i-ilb, j-jlb, k-klb, 1, xdim), smallr)
1732 qll(k,2) = max(qrt(i-irt, j-jrt, k-krt, 5, xdim), smallp)
1733 qrl(k,2) = max(qlt(i-ilt, j-jlt, k-klt, 5, xdim), smallp)
1734 qlr(k,2) = max(qrb(i-irb, j-jrb, k-krb, 5, xdim), smallp)
1735 qrr(k,2) = max(qlb(i-ilb, j-jlb, k-klb, 5, xdim), smallp)
1736 qll(k,3) = qrt(i-irt, j-jrt, k-krt, lp1, xdim)
1737 qrl(k,3) = qlt(i-ilt, j-jlt, k-klt, lp1, xdim)
1738 qlr(k,3) = qrb(i-irb, j-jrb, k-krb, lp1, xdim)
1739 qrr(k,3) = qlb(i-ilb, j-jlb, k-klb, lp1, xdim)
1740 qll(k,4) = qrt(i-irt, j-jrt, k-krt, lp2, xdim)
1741 qrl(k,4) = qlt(i-ilt, j-jlt, k-klt, lp2, xdim)
1742 qlr(k,4) = qrb(i-irb, j-jrb, k-krb, lp2, xdim)
1743 qrr(k,4) = qlb(i-ilb, j-jlb, k-klb, lp2, xdim)
1744 qll(k,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1745 qrl(k,6) = half*(qrt(i-irt, j-jrt, k-krt, bp1, xdim) + qlt(i-ilt, j-jlt, k-klt, bp1, xdim))
1746 qlr(k,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1747 qrr(k,6) = half*(qrb(i-irb, j-jrb, k-krb, bp1, xdim) + qlb(i-ilb, j-jlb, k-klb, bp1, xdim))
1748 qll(k,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1749 qrl(k,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1750 qlr(k,7) = half*(qrt(i-irt, j-jrt, k-krt, bp2, xdim) + qrb(i-irb, j-jrb, k-krb, bp2, xdim))
1751 qrr(k,7) = half*(qlt(i-ilt, j-jlt, k-klt, bp2, xdim) + qlb(i-ilb, j-jlb, k-klb, bp2, xdim))
1752 qll(k,5) = qrt(i-irt, j-jrt, k-krt, lor, xdim)
1753 qrl(k,5) = qlt(i-ilt, j-jlt, k-klt, lor, xdim)
1754 qlr(k,5) = qrb(i-irb, j-jrb, k-krb, lor, xdim)
1755 qrr(k,5) = qlb(i-ilb, j-jlb, k-klb, lor, xdim)
1756 qll(k,8) = qrt(i-irt, j-jrt, k-krt, bor, xdim)
1757 qrl(k,8) = qlt(i-ilt, j-jlt, k-klt, bor, xdim)
1758 qlr(k,8) = qrb(i-irb, j-jrb, k-krb, bor, xdim)
1759 qrr(k,8) = qlb(i-ilb, j-jlb, k-klb, bor, xdim)
1761 call solver_2d (self,qll,qrl,qlr,qrr,emf,i,j,l,u,nvar,xdim,detailed_timer=detailed_timer)
1765 call trace%end (itimer, detailed_timer=detailed_timer)
1767 END SUBROUTINE cmp_mag_flx
1772 SUBROUTINE solver (self, qleft, qright, fgdnv, ub)
1774 real(dp),
dimension(:,:):: qleft,qright,fgdnv
1777 if (self%riemann ==
'hlld')
then 1778 call rieman%hlld (qleft, qright, fgdnv, ub)
1779 else if (self%riemann ==
'llf')
then 1780 call rieman%llf (qleft, qright, fgdnv, ub)
1782 write(*,*)
'unknown Riemann solver' 1785 END SUBROUTINE solver
1790 SUBROUTINE solver_2d (self,qLL,qRL,qLR,qRR,emf,j,k,l,u,nvar,xdim,detailed_timer)
1792 real(DP),
dimension(:,:),
intent(in):: qll,qrl,qlr,qrr
1793 real(DP),
dimension(:,:,:):: emf
1794 integer,
intent(in):: j,k,l,u,nvar,xdim
1795 logical,
optional:: detailed_timer
1797 if (self%riemann ==
'hlld')
then 1798 call rieman%hlld_2d (qll,qrl,qlr,qrr,emf,j,k,l,u,nvar,xdim,detailed_timer=detailed_timer)
1799 else if (self%riemann ==
'llf')
then 1800 call rieman%llf_2d (qll,qrl,qlr,qrr,emf,j,k,l,u,nvar,xdim,detailed_timer=detailed_timer)
1802 write(*,*)
'unknown Riemann 2D solver' 1805 END SUBROUTINE solver_2d
1808 FUNCTION up(f,i)
RESULT (g)
1809 real,
dimension(:,:,:),
intent(in):: f
1810 real,
dimension(size(f,1),size(f,2),size(f,3)):: g
1813 g = 0.5*(cshift(f,1,i)+f)
1817 FUNCTION gas_pressure (self)
RESULT (pg)
1819 real,
dimension(self%gn(1),self%gn(2),self%gn(3)):: pg
1820 integer :: l(3), u(3)
1822 associate(d => self%mem(:,:,:,self%idx%d ,self%it,1), &
1823 px => self%mem(:,:,:,self%idx%px,self%it,1), &
1824 py => self%mem(:,:,:,self%idx%py,self%it,1), &
1825 pz => self%mem(:,:,:,self%idx%pz,self%it,1), &
1826 bx => self%mem(:,:,:,self%idx%bx,self%it,1), &
1827 by => self%mem(:,:,:,self%idx%by,self%it,1), &
1828 bz => self%mem(:,:,:,self%idx%bz,self%it,1), &
1829 e => self%mem(:,:,:,self%idx%e ,self%it,1))
1830 if (isothermal)
then 1831 pg = d*self%csound**2
1834 (e-0.5*(up(bx,1)**2+up(by,2)**2+up(bz,3)**2+(px**2+py**2+pz**2)/d))
1837 END FUNCTION gas_pressure
1842 SUBROUTINE gravitational_force (self)
1844 real,
dimension(:,:,:),
pointer:: phi, d, et, gradphi1, gradphi2, gradphi3
1845 real,
dimension(:,:,:,:),
pointer:: p
1847 real,
dimension(3):: wtp1, wtm1
1850 associate(l => self%li, u => self%ui)
1851 wtp1 = 0.50 * 4.0 / 3.0 / self%ds
1852 wtm1 = 0.25 * 1.0 / 3.0 / self%ds
1853 phi => self%mem(:,:,:,self%idx%phi,self%it,1)
1854 d => self%mem(:,:,:,self%idx%d,self%it,1)
1855 et => self%mem(:,:,:,self%idx%e,self%new,1)
1856 p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%new,1)
1857 allocate (gradphi1(self%gn(1),self%gn(2),self%gn(3)), &
1858 gradphi2(self%gn(1),self%gn(2),self%gn(3)), &
1859 gradphi3(self%gn(1),self%gn(2),self%gn(3)))
1865 etemp = et(i,j,k) - 0.5 * p(i,j,k,1)**2 / d(i,j,k) &
1866 - 0.5 * p(i,j,k,2)**2 / d(i,j,k) &
1867 - 0.5 * p(i,j,k,3)**2 / d(i,j,k)
1869 gradphi1(i,j,k) = wtp1(1) * (phi(i+1,j,k) - phi(i-1,j,k)) &
1870 - wtm1(1) * (phi(i+2,j,k) - phi(i-2,j,k))
1871 gradphi2(i,j,k) = wtp1(2) * (phi(i,j+1,k) - phi(i,j-1,k)) &
1872 - wtm1(2) * (phi(i,j+2,k) - phi(i,j-2,k))
1873 gradphi3(i,j,k) = wtp1(3) * (phi(i,j,k+1) - phi(i,j,k-1)) &
1874 - wtm1(3) * (phi(i,j,k+2) - phi(i,j,k-2))
1876 p(i,j,k,1) = p(i,j,k,1) - self%dtime * d(i,j,k) * gradphi1(i,j,k)
1877 p(i,j,k,2) = p(i,j,k,2) - self%dtime * d(i,j,k) * gradphi2(i,j,k)
1878 p(i,j,k,3) = p(i,j,k,3) - self%dtime * d(i,j,k) * gradphi3(i,j,k)
1880 et(i,j,k) = etemp + 0.5 * p(i,j,k,1)**2 / d(i,j,k) &
1881 + 0.5 * p(i,j,k,2)**2 / d(i,j,k) &
1882 + 0.5 * p(i,j,k,3)**2 / d(i,j,k)
1887 deallocate (gradphi1, gradphi2, gradphi3)
1889 END SUBROUTINE gravitational_force
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
RAMSES Godunov solvers, use of guard zones; specifically in HLLD.
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
Module with list handling for generic class task_t objects.
The lock module uses nested locks, to allow versatile use of locks, where a procedure may want to mak...