34 real,
dimension(:,:,:),
pointer:: d, dddt, s, dsdt, pg
35 real,
dimension(:,:,:,:),
pointer:: p, dpdt, b, dbdt
39 logical:: first_time=.true.
40 logical:: do_force=.true.
48 procedure:: gas_pressure
49 procedure:: gas_pressure_sub
52 logical,
save:: flop_count=.false.
53 logical,
save:: do_2nd_div=.true.
54 logical,
save:: do_smooth=.false.
55 logical,
save:: do_maxwell=.true.
62 SUBROUTINE init (self)
65 real,
save:: nu(6)=[0.1,1.0,0.0,0.5,0.5,0.5], csound=1.0, cdtd=0.5, courant=0.2
66 real(8),
save:: gamma=1.4_8
67 logical,
save:: first_time=.true.
68 character(len=16),
save:: eos
69 namelist /stagger_params/ nu, courant, gamma, csound, cdtd, hardwire, eos, &
70 do_smooth, do_2nd_div, flop_count, verbose, do_maxwell
72 call trace%begin(
'mhd_t%init')
77 read (io%input, stagger_params)
78 if (io%master)
write (*, stagger_params)
82 if (timestep%time_order <= 0)
call mpi%abort(
"The Stagger solvers require `time_order` > 0. Abort!")
89 self%courant = courant
95 call self%idx%init (5, self%mhd)
96 self%initial%mhd = self%mhd
97 call self%initial%init (self%kind)
98 self%mhd = self%initial%mhd
100 self%kind =
'stagger2_mhd_patch' 102 self%kind =
'stagger2_hd_patch' 117 call self%patch_t%init
118 call self%idx%init (5, self%mhd)
121 self%mesh(i)%h(self%idx%px+i-1) = -0.5
122 if (self%mhd) self%mesh(i)%h(self%idx%bx+i-1) = -0.5
124 if (io%verbose>1 .and. self%track)
then 126 print
'("h:",8f5.2)', self%mesh(i)%h
129 call self%gpatch_t%init
130 self%unsigned(self%idx%d) = .true.
131 self%pervolume(self%idx%s) = .true.
132 call self%extras_t%init
147 SUBROUTINE test (self)
149 integer :: ix, iy, iz, m(3), n(3), n0(3), i1
150 real :: fx, fy, fz, eps
152 class(
mhd_t),
pointer:: tmp
156 call trace%begin(
'mhd_t%test')
160 call tmp%idx%init (5, tmp%mhd)
161 call tmp%patch_t%init
163 if (any(tmp%n /= n0))
then 165 print *,
'mhd_t%test not possible with n, no_mans_land =', tmp%n, tmp%no_mans_land
174 do iz=tmp%mesh(3)%lb,tmp%mesh(3)%ub
175 fz = sin(4.*cgs%pi*tmp%mesh(3)%r(iz)/tmp%size(3))
176 do iy=tmp%mesh(2)%lb,tmp%mesh(2)%ub
177 fy = sin(4.*cgs%pi*tmp%mesh(2)%r(iy)/tmp%size(2))
178 do ix=tmp%mesh(1)%lb,tmp%mesh(1)%ub
179 fx = sin(4.*cgs%pi*tmp%mesh(1)%r(ix)/tmp%size(1))
180 tmp%mem(ix,iy,iz,tmp%idx%d,tmp%it,1) = 1.0+0.1*(fx+fy+fz)
181 tmp%mem(ix,iy,iz,tmp%idx%s,tmp%it,1) = 0.1*(fx+fy+fz)
182 tmp%mem(ix,iy,iz,tmp%idx%px,tmp%it,1) = fx
183 tmp%mem(ix,iy,iz,tmp%idx%py,tmp%it,1) = fy
184 tmp%mem(ix,iy,iz,tmp%idx%pz,tmp%it,1) = fz
186 tmp%mem(ix,iy,iz,tmp%idx%bx,tmp%it,1) = fy
187 tmp%mem(ix,iy,iz,tmp%idx%by,tmp%it,1) = fz
188 tmp%mem(ix,iy,iz,tmp%idx%bz,tmp%it,1) = fx
197 tmp%do_force = .false.
199 tmp%do_force = .true.
202 eps = 3e-6/minval(tmp%ds)
204 print *,
'minval(ds), eps =', minval(tmp%ds), eps
205 n = nint(tmp%size*0.5_8/tmp%ds)
206 do ix=tmp%mesh(1)%li,tmp%mesh(1)%ui
207 i1 = mod(ix-tmp%mesh(1)%li+n(1),2*n(1)) + tmp%mesh(1)%li
208 ok = all(abs(tmp%mem(ix,m(2),m(3),1:5,tmp%it,2)- &
209 tmp%mem(i1,m(2),m(3),1:5,tmp%it,2)) < eps)
210 allok = allok .and. ok
211 if (io%verbose>0 .or. io_unit%do_validate .or. .not.ok) &
212 print 1,ix,tmp%mem(ix,m(2),m(3),1:5,tmp%it,2), &
213 tmp%mem(i1,m(2),m(3),1:5,tmp%it,2)
214 1
format(
"mhd_t%test:",i4,1p,2(2x,5e14.5))
216 do iy=tmp%mesh(2)%li,tmp%mesh(2)%ui
217 i1 = mod(iy-tmp%mesh(2)%li+n(2),2*n(2)) + tmp%mesh(2)%li
218 ok = all(abs(tmp%mem(m(1),iy,m(3),1:5,tmp%it,2)- &
219 tmp%mem(m(1),i1,m(3),1:5,tmp%it,2)) < eps)
220 allok = allok .and. ok
221 if (io%verbose>0 .or. io_unit%do_validate .or. .not.ok) &
222 print 1,iy,tmp%mem(m(1),iy,m(3),1:5,tmp%it,2), &
223 tmp%mem(m(1),i1,m(3),1:5,tmp%it,2)
225 do iz=tmp%mesh(3)%li,tmp%mesh(3)%ui
226 i1 = mod(iz-tmp%mesh(3)%li+n(3),2*n(3)) + tmp%mesh(3)%li
227 ok = all(abs(tmp%mem(m(1),m(2),iz,1:5,tmp%it,2)- &
228 tmp%mem(m(1),m(2),i1,1:5,tmp%it,2)) < eps)
229 allok = allok .and. ok
230 if (io%verbose>0 .or. io_unit%do_validate .or. .not.ok) &
231 print 1,iz,tmp%mem(m(1),m(2),iz,1:5,tmp%it,2), &
232 tmp%mem(m(1),m(2),i1,1:5,tmp%it,2)
235 print *,
'mhd_t%test passed' 238 print *,
'mhd_t%test failed' 240 call io%abort(
'mhd_t%test')
248 SUBROUTINE pde (self)
251 real,
dimension(:,:,:),
pointer:: d, s, dddt, dsdt
252 real,
dimension(:,:,:),
pointer:: dpdtx, dpdty, dpdtz, dbdtx, dbdty, dbdtz
253 real,
dimension(:,:,:),
pointer:: bx, by, bz, px, py, pz, phi
254 real,
dimension(:,:,:),
pointer:: ux, uy, uz, pg
255 real,
dimension(:,:,:,:),
pointer:: u, p, b, dpdt, dbdt
256 real,
dimension(:,:,:),
pointer:: ex, ey, ez, jx, jy, jz
257 real,
dimension(:,:,:),
pointer:: fdx, fdy, fdz, ldx, ldy, ldz, ddx, ddy, ddz
258 real,
dimension(:,:,:),
pointer:: fxx, fxy
259 real,
dimension(:,:,:),
pointer:: fyy, fyz
260 real,
dimension(:,:,:),
pointer:: fzx, fzz
261 real,
dimension(:,:,:),
allocatable:: du, lnd, ss, cs, pa, pp, pb, q, uu
262 real,
dimension(:,:,:),
allocatable:: txx, txy
263 real,
dimension(:,:,:),
allocatable:: tyy, tyz
264 real,
dimension(:,:,:),
allocatable:: tzx, tzz
265 real,
dimension(:,:,:,:),
allocatable,
target:: emf, j
266 real,
dimension(:,:,:,:),
allocatable,
target:: fx, fy, fz, fd, ld, dd, fm
267 integer :: omp_get_num_threads, omp_get_thread_num
268 logical,
save:: first_time=.true.
270 real(8):: ds(3), dsmax
271 integer,
save:: itimer=0
274 call trace%begin(
'mhd_t%pde',itimer=itimer)
275 associate(i => self%idx)
276 self%d => self%mem(:,:,:,i%d ,self%it,1); self%dddt => self%mem(:,:,:,i%d ,self%it,2)
277 self%s => self%mem(:,:,:,i%s ,self%it,1); self%dsdt => self%mem(:,:,:,i%s ,self%it,2)
278 self%p => self%mem(:,:,:,i%px:i%pz,self%it,1); self%dpdt => self%mem(:,:,:,i%px:i%pz,self%it,2)
279 d=>self%d; dddt=>self%dddt; s=>self%s; dsdt=>self%dsdt;
280 p=>self%p; dpdt=>self%dpdt
281 px=>p(:,:,:,1); py=>p(:,:,:,2); pz=>p(:,:,:,3);
282 dpdtx=>dpdt(:,:,:,1); dpdty=>dpdt(:,:,:,2); dpdtz=>dpdt(:,:,:,3);
284 self%B => self%mem(:,:,:,i%bx:i%bz,self%it,1); self%dBdt => self%mem(:,:,:,i%bx:i%bz,self%it,2)
285 b=>self%B; dbdt=>self%dBdt
286 bx=>b(:,:,:,1); by=>b(:,:,:,2); bz=>b(:,:,:,3);
287 dbdtx=>dbdt(:,:,:,1); dbdty=>dbdt(:,:,:,2); dbdtz=>dbdt(:,:,:,3);
290 call allocate_vectors_a (self%gn, ld, dd, emf, j, fd, fx, fy, fz, fm)
291 call allocate_scalars_a (self%gn, lnd, pa, du, ss, cs, pp, pb, q)
292 if (self%do_pebbles)
then 293 u => self%vgas(:,:,:,:,self%it)
294 pg => self%pressure(:,:,:,self%it)
295 self%density => self%mem(:,:,:,1,:,1)
297 call allocate_scalars (self%gn, pg)
298 call allocate_vectors (self%gn, u)
300 ux=>u(:,:,:,1); uy=>u(:,:,:,2); uz=>u(:,:,:,3)
301 ex=>emf(:,:,:,1); ey=>emf(:,:,:,2); ez=>emf(:,:,:,3)
302 jx=>j(:,:,:,1); jy=>j(:,:,:,2); jz=>j(:,:,:,3)
303 fxx=>fx(:,:,:,1); fxy=>fx(:,:,:,2)
304 fyy=>fy(:,:,:,2); fyz=>fy(:,:,:,3)
305 fzx=>fz(:,:,:,1); fzz=>fz(:,:,:,3)
306 fdx=>fd(:,:,:,1); fdy=>fd(:,:,:,2); fdz=>fd(:,:,:,3)
307 ldx=>ld(:,:,:,1); ldy=>ld(:,:,:,2); ldz=>ld(:,:,:,3)
308 ddx=>dd(:,:,:,1); ddy=>dd(:,:,:,2); ddz=>dd(:,:,:,3)
321 call self%gas_pressure_sub (d, s, ss, pg)
327 call allocate_scalars_a (self%gn, txx, tyy, tzz)
337 dsmax = maxval(ds,self%n > 2)
338 du = -dsmax*min(du,0.0)
344 call allocate_scalars_a (self%gn, uu)
345 uu = sqrt(fxx+fyy+fzz)
348 pa = self%nu(2)*d*stagger%xyzsm(du)**2
350 pa = self%nu(2)*d*du**2
362 cs = sqrt((self%gamma*pg+2.0*pb + pa)/d)
364 u_max = self%fmaxval(cs,outer=.true.)
365 self%u_max = max(self%u_max,u_max)
366 print 1, self%id,
'u_max(ca) ', u_max, self%u_max
367 1
format(i6,2x,a,1p,2g12.3)
373 du = (dsmax*self%nu(1))*(cs + 2.*uu)
376 fxx = pp + d*(fxx - du*txx) - ex
377 fyy = pp + d*(fyy - du*tyy) - ey
378 fzz = pp + d*(fzz - du*tzz) - ez
380 fxx = pp + d*(fxx - du*txx)
381 fyy = pp + d*(fyy - du*tyy)
382 fzz = pp + d*(fzz - du*tzz)
389 cs = sqrt((self%gamma*pg+pa)/d)
391 u_max = self%fmaxval(cs,outer=.true.)
392 self%u_max = max(self%u_max,u_max)
393 print 1, self%id,
'u_max(cs) ', u_max, self%u_max
402 du = (dsmax*self%nu(1))*(cs + 2.*uu)
403 fxx = pp + d*(fxx - du*txx)
404 fyy = pp + d*(fyy - du*tyy)
405 fzz = pp + d*(fzz - du*tzz)
408 print *, stagger%count,
' stagger calls', stagger%flops/product(
real(self%n)),
' flops per interior point' 409 u_max = self%fmaxval(cs + uu,outer=.true.)
410 self%u_max = max(self%u_max,u_max)
411 if (verbose>2) print 1, self%id,
'u_max(cs+u) ', u_max, self%u_max
416 q = d*du*(txx**2 + tyy**2 + tzz**2)
417 call deallocate_scalars_a (txx, tyy, tzz, uu)
418 if (verbose==1) print *,
'average(Q):', self%faver(q)
422 fd = -(self%nu(4)*dd*
down(du))*
ddown(ds,lnd)
425 u_max = self%cdtd*dsmax*self%fmaxval(abs(dddt/d),outer=.true.)
426 self%u_max = max(self%u_max,u_max)
432 fxy = ydn(ux)*xdn(uy)
433 fyz = zdn(uy)*ydn(uz)
434 fzx = xdn(uz)*zdn(ux)
440 call allocate_scalars_a (self%gn, txy, tyz, tzx)
441 txy = ddxdn(ds,uy)+ddydn(ds,ux)
442 tyz = ddydn(ds,uz)+ddzdn(ds,uy)
443 tzx = ddzdn(ds,ux)+ddxdn(ds,uz)
444 fxy = fxy - self%nu(3)*0.5*xdn1(ydn1(du))*txy
445 fyz = fyz - self%nu(3)*0.5*ydn1(zdn1(du))*tyz
446 fzx = fzx - self%nu(3)*0.5*zdn1(xdn1(du))*tzx
447 q = q + self%nu(3)*.25*d*du*(xup(yup(txy))**2 + &
450 call deallocate_scalars_a (txy, tyz, tzx)
454 fxy =
exp(ydn(ldx))*fxy
455 fyz =
exp(zdn(ldy))*fyz
456 fzx =
exp(xdn(ldz))*fzx
460 if (self%mhd.and.do_maxwell)
then 461 fxy = fxy - ydn(bx)*xdn(by)
462 fyz = fyz - zdn(by)*ydn(bz)
463 fzx = fzx - xdn(bz)*zdn(bx)
470 if (self%nu(4) > 0.0)
then 475 fxx = fxx - pa*ddxup(ds,ldx)*xup(ux)
476 fyy = fyy - pa*ddyup(ds,ldy)*yup(uy)
477 fzz = fzz - pa*ddzup(ds,ldz)*zup(uz)
478 fxy = fxy + (xdn1(fdy)*ydn(ux)+ydn1(fdx)*xdn(uy))
479 fyz = fyz + (ydn1(fdz)*zdn(uy)+zdn1(fdy)*ydn(uz))
480 fzx = fzx + (zdn1(fdx)*xdn(uz)+xdn1(fdz)*zdn(ux))
487 dpdtx = - ddxdn1(ds,fxx) - ddyup1(ds,fxy) - ddzup1(ds,fzx)
488 dpdty = - ddxup1(ds,fxy) - ddydn1(ds,fyy) - ddzup1(ds,fyz)
489 dpdtz = - ddxup1(ds,fzx) - ddyup1(ds,fyz) - ddzdn1(ds,fzz)
491 dpdtx = - ddxdn(ds,fxx) - ddyup(ds,fxy) - ddzup(ds,fzx)
492 dpdty = - ddxup(ds,fxy) - ddydn(ds,fyy) - ddzup(ds,fyz)
493 dpdtz = - ddxup(ds,fzx) - ddyup(ds,fyz) - ddzdn(ds,fzz)
499 if (self%do_force)
then 500 if (
allocated(self%force_per_unit_mass))
then 501 dpdt = dpdt + dd*self%force_per_unit_mass
503 if (
allocated(self%force_per_unit_volume))
then 504 dpdt = dpdt + self%force_per_unit_volume
512 if (self%gamma==1d0)
then 522 dsdt = -
div(ds,fm) + q*d/pg
532 emf = self%nu(6)*
edge(du)*j
533 q = yup(zup(jx*ex)) + zup(xup(jy*ey)) + xup(yup(jz*ez))
537 if (.not.do_maxwell)
then 543 dpdtx = dpdtx + zup(jy)*xdn(fzz) - yup(jz)*xdn(fyy)
544 dpdty = dpdty + xup(jz)*ydn(fxx) - zup(jx)*ydn(fzz)
545 dpdtz = dpdtz + yup(jx)*zdn(fyy) - xup(jy)*zdn(fxx)
550 ex = ex + ydn(uz)*zdn(by) - zdn(uy)*ydn(bz)
551 ey = ey + zdn(ux)*xdn(bz) - xdn(uz)*zdn(bx)
552 ez = ez + xdn(uy)*ydn(bx) - ydn(ux)*xdn(by)
553 call non_ideal%update (self%gn, minval(ds), d, jx, jy, jz, bx, by, bz, q, emf, self%u_max)
554 if (self%gamma==1d0)
then 564 u_max = self%cdtd*dsmax*self%fmaxval(abs((dsdt-s*dddt/d)/d),outer=.true.)
565 self%u_max = max(self%u_max,u_max)
566 if (verbose>2) print 1, self%id,
'u_max(dsdt) ', u_max, self%u_max
568 call deallocate_vectors_a (ld, dd, emf, j, fd, fx, fy, fz, fm)
569 call deallocate_scalars_a (lnd, pa, du, ss, cs, pp, pb, q)
570 if (.not.self%do_pebbles)
then 571 call deallocate_vectors (u)
572 call deallocate_scalars (pg)
574 call trace%end (itimer)
575 if (mpi%master .and. (first_time .or. flop_count))
then 576 print *, stagger%count,
' stagger calls', stagger%flops/product(
real(self%n)),
' flops per interior point' 587 SUBROUTINE update (self)
589 integer,
save:: itimer=0
591 call trace%begin(
'mhd_t%update',itimer=itimer)
599 call self%courant_condition
601 call self%check_density(
'before timestep')
602 call timestep%update (self%id, self%iit, self%dt, self%mem, self%mesh, self%lock)
603 call self%check_density(
' after timestep')
605 call self%counter_update
606 call trace%end (itimer)
607 END SUBROUTINE update
612 SUBROUTINE output (self)
614 character(len=64):: filename
618 out_next = self%out_next
619 call trace%begin(
'mhd_t%output')
620 call self%extras_t%output
621 if (self%time >= self%out_next .and. io%do_output .and. io%do_legacy)
then 623 write (filename,
'(a,i4.4,"_",i4.4,".txt")') trim(io%outputname), &
625 inquire(file=filename, exist=exists)
627 open (io%data_unit, file=trim(filename), form=
'formatted', status=
'old', &
629 write (io%data_unit,
'(8f10.4)') self%nu
635 END SUBROUTINE output
638 FUNCTION gas_pressure (self)
RESULT (pg)
640 real,
dimension(:,:,:),
pointer:: d, s
641 real,
dimension(self%gn(1),self%gn(2),self%gn(3)):: pg
642 integer,
save:: nprint=3, n(3)
644 if (trim(self%eos)==
'ideal')
then 645 d => self%mem(:,:,:,self%idx%d,self%it,1)
646 s => self%mem(:,:,:,self%idx%s,self%it,1)
647 if (self%gamma==1d0)
then 649 pg = d*self%csound**2
651 s => self%mem(:,:,:,self%idx%s,self%it,1)
652 pg = d**self%gamma*exp(s/d*(self%gamma-1))
658 call eos%lookup (n, x=s, y=d)
660 call eos%lookup (n, x=s, y=d, pg=pg)
663 if (io%master .and. nprint>0)
then 665 print *,
'mhd_t%gas_pressure: min, max =', self%fminval(pg), self%fmaxval(pg)
667 END FUNCTION gas_pressure
670 SUBROUTINE gas_pressure_sub (self, d, s, ss, pg)
672 real,
dimension(:,:,:),
pointer:: d, s, pg
673 real,
dimension(:,:,:):: ss
675 integer,
save:: nprint=3
677 if (trim(self%eos)==
'ideal')
then 678 if (self%gamma==1d0)
then 680 pg = d*self%csound**2
682 pg = exp(self%gamma*log(d))*exp(ss*(self%gamma-1.0))
686 call eos%lookup (shape(d), x=s, y=d, pg=pg)
688 if (io%master .and. nprint>0)
then 690 print *,
'mhd_t%gas_pressure: min, max =', self%fminval(pg), self%fmaxval(pg)
692 END SUBROUTINE gas_pressure_sub
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
Define code units, in terms of (here) CGS units.
RAMSES Godunov solvers, use of guard zones; specifically in HLLD.
Equation of state module for any sort of tables, provided by a reader.
Fundamental constants in CGS and SI units.
Simple initical condition module for tests.
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
6th order stagger operators, with self-test procedure