36 real,
dimension(:,:,:),
pointer :: d, dddt, e, dedt
37 real,
dimension(:,:,:,:),
pointer :: p, dpdt, b, dbdt
38 real,
dimension(:,:,:),
pointer :: lnd, ee, pg=>null(), tt=>null()
39 real,
dimension(:,:,:),
allocatable:: qrad
40 procedure(divb_clean1),
pointer :: divb_clean => null()
45 logical:: first_time=.true.
50 procedure:: pre_update
53 procedure:: gas_pressure
54 procedure:: gas_pressure_sub
56 logical,
save:: unsigned=.true.
57 logical,
save:: do_test=.false.
58 logical,
save:: do_smooth=.false.
59 logical,
save:: do_maxwell = .true.
60 logical,
save:: do_temperature = .false.
61 integer,
save:: verbose=0
62 integer,
save:: divb_cleaner=2
69 SUBROUTINE pre_init (self)
70 class(
mhd_t),
target:: self
73 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
74 real(8),
save:: gamma=1.4_8
75 logical,
save:: first_time=.true.
76 logical,
save:: mhd=.false.
77 character(len=16),
save:: eos
78 namelist /stagger_params/ nu, courant, gamma, csound, cdtd, hardwire, eos, &
79 do_smooth, unsigned, do_maxwell, do_temperature, verbose, divb_cleaner, &
82 call trace%begin(
'mhd_t%pre_init')
89 read (io%input, stagger_params)
90 if (io%master)
write (*, stagger_params)
94 if (timestep%time_order <= 0)
call mpi%abort(
"The Stagger solvers require `time_order` > 0. Abort!")
102 self%courant = courant
108 self%initial%mhd = self%mhd
109 call self%initial%init (self%kind,
real(self%gamma))
110 self%mhd = self%initial%mhd
112 self%kind =
'stagger2e_mhd_patch' 114 self%kind =
'stagger2e_hd_patch' 116 if (self%nv == 0)
then 123 self%nv = self%nv + 5
124 if (self%mhd) self%nv = self%nv + 3
126 select case(divb_cleaner)
128 self%divb_clean => divb_clean1
130 self%divb_clean => divb_clean2
133 call self%idx%init (5, self%mhd)
137 call self%extras_t%pre_init
139 END SUBROUTINE pre_init
145 SUBROUTINE init (self)
149 call trace%begin(
'mhd_t%init')
154 call self%patch_t%init
156 self%mesh(i)%h(self%idx%px+i-1) = -0.5
157 if (self%mhd) self%mesh(i)%h(self%idx%bx+i-1) = -0.5
159 if (verbose>1 .and. self%track)
then 161 print
'("h:",8f5.2)', self%mesh(i)%h
164 call self%patch_t%init
165 select case(divb_cleaner)
167 self%divb_clean => divb_clean1
169 self%divb_clean => divb_clean2
171 call io%abort(
'mhd_mod:: divb cleaner method not understood')
173 call self%gpatch_t%init
174 self%unsigned(self%idx%d) = unsigned
175 self%pervolume(self%idx%e) = unsigned
179 call self%extras_t%init
184 SUBROUTINE pde (self)
188 real,
dimension(:,:,:),
pointer:: d, e, dddt, dedt
189 real,
dimension(:,:,:),
pointer:: dpdtx, dpdty, dpdtz, dbdtx, dbdty, dbdtz
190 real,
dimension(:,:,:),
pointer:: bx, by, bz, px, py, pz, phi, qr
191 real,
dimension(:,:,:),
pointer:: ux, uy, uz, pg
192 real,
dimension(:,:,:,:),
pointer:: u, p, b, dpdt, dbdt
193 real,
dimension(:,:,:),
pointer:: ex, ey, ez, jx, jy, jz
194 real,
dimension(:,:,:),
pointer:: fdx, fdy, fdz, ldx, ldy, ldz, ddx, ddy, ddz
195 real,
dimension(:,:,:),
pointer:: fxx, fxy
196 real,
dimension(:,:,:),
pointer:: fyy, fyz
197 real,
dimension(:,:,:),
pointer:: fzx, fzz
198 real,
dimension(:,:,:),
allocatable:: du, ss, cs, pa, pp, pb, q, uu, divu
199 real,
dimension(:,:,:),
allocatable,
target:: ee, lnd
200 real,
dimension(:,:,:),
allocatable:: txx, txy
201 real,
dimension(:,:,:),
allocatable:: tyy, tyz
202 real,
dimension(:,:,:),
allocatable:: tzx, tzz
203 real,
dimension(:,:,:,:),
allocatable,
target:: emf, j
204 real,
dimension(:,:,:,:),
allocatable,
target:: fx, fy, fz, fd, ld, dd, fm
205 integer :: omp_get_num_threads, omp_get_thread_num
206 logical,
save:: first_time=.true.
208 real(8):: ds(3), dsmax
210 integer,
save:: itimer=0
212 call trace%begin(
'mhd_t%pde',itimer=itimer)
213 associate(i => self%idx)
214 self%d => self%mem(:,:,:,i%d ,self%it,1); self%dddt => self%mem(:,:,:,i%d ,self%it,2)
215 self%e => self%mem(:,:,:,i%e ,self%it,1); self%dedt => self%mem(:,:,:,i%e ,self%it,2)
216 self%p => self%mem(:,:,:,i%px:i%pz,self%it,1); self%dpdt => self%mem(:,:,:,i%px:i%pz,self%it,2)
217 d=>self%d; dddt=>self%dddt; e=>self%e; dedt=>self%dedt;
218 p=>self%p; dpdt=>self%dpdt
219 px=>p(:,:,:,1); py=>p(:,:,:,2); pz=>p(:,:,:,3);
220 dpdtx=>dpdt(:,:,:,1); dpdty=>dpdt(:,:,:,2); dpdtz=>dpdt(:,:,:,3);
222 self%B => self%mem(:,:,:,i%bx:i%bz,self%it,1); self%dBdt => self%mem(:,:,:,i%bx:i%bz,self%it,2)
223 b=>self%B; dbdt=>self%dBdt
224 bx=>b(:,:,:,1); by=>b(:,:,:,2); bz=>b(:,:,:,3);
225 dbdtx=>dbdt(:,:,:,1); dbdty=>dbdt(:,:,:,2); dbdtz=>dbdt(:,:,:,3);
228 call allocate_vectors_a (self%gn, ld, dd, emf, j, fd, fx, fy, fz, fm)
229 call allocate_scalars_a (self%gn, lnd, pa, du, ee, cs, pp, pb, q, uu, divu)
230 if (self%do_pebbles)
then 231 u => self%vgas(:,:,:,:,self%it)
232 pg => self%pressure(:,:,:,self%it)
233 self%density => self%mem(:,:,:,1,:,1)
235 call allocate_scalars (self%gn, pg)
236 call allocate_vectors (self%gn, u)
238 ux=>u(:,:,:,1); uy=>u(:,:,:,2); uz=>u(:,:,:,3)
239 ex=>emf(:,:,:,1); ey=>emf(:,:,:,2); ez=>emf(:,:,:,3)
240 jx=>j(:,:,:,1); jy=>j(:,:,:,2); jz=>j(:,:,:,3)
241 fxx=>fx(:,:,:,1); fxy=>fx(:,:,:,2)
242 fyy=>fy(:,:,:,2); fyz=>fy(:,:,:,3)
243 fzx=>fz(:,:,:,1); fzz=>fz(:,:,:,3)
244 fdx=>fd(:,:,:,1); fdy=>fd(:,:,:,2); fdz=>fd(:,:,:,3)
245 ldx=>ld(:,:,:,1); ldy=>ld(:,:,:,2); ldz=>ld(:,:,:,3)
246 ddx=>dd(:,:,:,1); ddy=>dd(:,:,:,2); ddz=>dd(:,:,:,3)
262 call self%gas_pressure_sub (d, self%lnd, self%ee, pg)
268 call allocate_scalars_a (self%gn, txx, tyy, tzz)
278 dsmax = maxval(ds,self%n > 2)
279 du = -dsmax*min(divu,0.0)
281 pa = self%nu(2)*d*xyzsm(du**2)
283 pa = self%nu(2)*d*du**2
291 uu = sqrt(fxx+fyy+fzz)
301 cs = sqrt((self%gamma*pg+2.0*pb + pa)/d)
303 u_max = self%fmaxval(cs,outer=.true.)
304 self%u_max = max(self%u_max,u_max)
305 print 1,
'u_max(ca) ', u_max, self%u_max
306 1
format(1x,a,1p,2g12.3)
311 du = (dsmax*self%nu(1))*(cs + 2.*uu)
314 fxx = pp + d*(fxx - du*txx) - ex
315 fyy = pp + d*(fyy - du*tyy) - ey
316 fzz = pp + d*(fzz - du*tzz) - ez
318 fxx = pp + d*(fxx - du*txx)
319 fyy = pp + d*(fyy - du*tyy)
320 fzz = pp + d*(fzz - du*tzz)
324 cs = sqrt((self%gamma*pg+pa)/d)
326 u_max = self%fmaxval(cs,outer=.true.)
327 self%u_max = max(self%u_max,u_max)
328 print 1,
'u_max(cs) ', u_max, self%u_max
331 du = (dsmax*self%nu(1))*(cs + 2.*uu)
336 fxx = pp + d*(fxx - du*txx)
337 fyy = pp + d*(fyy - du*tyy)
338 fzz = pp + d*(fzz - du*tzz)
340 u_max = self%fmaxval(cs+uu,outer=.true.)
341 self%u_max = max(self%u_max,u_max)
342 if (verbose>2) print 1,
'u_max(cs+u)', u_max, self%u_max
347 q = d*du*(txx**2 + tyy**2 + tzz**2)
353 fd = -(self%nu(4)*dd*
down(du))*
ddown(ds,lnd)
356 u_max = self%cdtd*dsmax*self%fmaxval(abs(dddt/d),outer=.true.)
357 self%u_max = max(self%u_max,u_max)
358 if (verbose>2) print 1,
'u_max(dddt)', u_max, self%u_max
363 fxy = ydn(ux)*xdn(uy)
364 fyz = zdn(uy)*ydn(uz)
365 fzx = xdn(uz)*zdn(ux)
369 call allocate_scalars_a (self%gn, txy, tyz, tzx)
370 txy = ddxdn(ds,uy)+ddydn(ds,ux)
371 tyz = ddydn(ds,uz)+ddzdn(ds,uy)
372 tzx = ddzdn(ds,ux)+ddxdn(ds,uz)
373 fxy = fxy - self%nu(3)*0.5*xdn1(ydn1(du))*txy
374 fyz = fyz - self%nu(3)*0.5*ydn1(zdn1(du))*tyz
375 fzx = fzx - self%nu(3)*0.5*zdn1(xdn1(du))*tzx
376 q = q + self%nu(3)*0.5*d*du*(xup(yup(txy))**2 + &
379 call deallocate_scalars_a (txy, tyz, tzx)
383 fxy =
exp(ydn(ldx))*fxy
384 fyz =
exp(zdn(ldy))*fyz
385 fzx =
exp(xdn(ldz))*fzx
389 if (self%mhd.and.do_maxwell)
then 390 fxy = fxy - ydn(bx)*xdn(by)
391 fyz = fyz - zdn(by)*ydn(bz)
392 fzx = fzx - xdn(bz)*zdn(bx)
399 if (self%nu(4) > 0.0)
then 404 fxx = fxx - pa*ddxup(ds,ldx)*xup(ux)
405 fyy = fyy - pa*ddyup(ds,ldy)*yup(uy)
406 fzz = fzz - pa*ddzup(ds,ldz)*zup(uz)
407 fxy = fxy + (xdn1(fdy)*ydn(ux)+ydn1(fdx)*xdn(uy))
408 fyz = fyz + (ydn1(fdz)*zdn(uy)+zdn1(fdy)*ydn(uz))
409 fzx = fzx + (zdn1(fdx)*xdn(uz)+xdn1(fdz)*zdn(ux))
415 dpdtx = - ddxdn1(ds,fxx) - ddyup1(ds,fxy) - ddzup1(ds,fzx)
416 dpdty = - ddxup1(ds,fxy) - ddydn1(ds,fyy) - ddzup1(ds,fyz)
417 dpdtz = - ddxup1(ds,fzx) - ddyup1(ds,fyz) - ddzdn1(ds,fzz)
422 if (
allocated(self%force_per_unit_mass))
then 423 dpdt = dpdt + dd*self%force_per_unit_mass
425 if (
allocated(self%force_per_unit_volume))
then 426 dpdt = dpdt + self%force_per_unit_volume
431 if (self%gamma == 1d0)
then 436 afm = self%faver(fm(:,:,:,1))
438 dedt = -
div(ds,fm) - pg*divu + q
439 if (verbose==1) print
'(a,10f12.6)',
'average(Q):', &
440 self%time, self%time*self%faver(q), &
441 self%faver(e-d**self%gamma/(self%gamma-1.0)), &
442 self%faver(- pg*divu), &
443 afm, self%faver(fm(:,:,:,1))
446 call deallocate_scalars_a (txx, tyy, tzz)
455 self%maxJ=maxval(abs(j))
456 emf = self%nu(6)*edge(du)*j
457 q = yup(zup(jx*ex)) + zup(xup(jy*ey)) + xup(yup(jz*ez))
458 if (self%gamma/=1.0) &
463 if (.not.do_maxwell)
then 468 dpdtx = dpdtx + zup(jy)*xdn(fzz) - yup(jz)*xdn(fyy)
469 dpdty = dpdty + xup(jz)*ydn(fxx) - zup(jx)*ydn(fzz)
470 dpdtz = dpdtz + yup(jx)*zdn(fyy) - xup(jy)*zdn(fxx)
475 ex = ex + ydn(uz)*zdn(by) - zdn(uy)*ydn(bz)
476 ey = ey + zdn(ux)*xdn(bz) - xdn(uz)*zdn(bx)
477 ez = ez + xdn(uy)*ydn(bx) - ydn(ux)*xdn(by)
478 call non_ideal%update (self%gn, minval(ds), d, jx, jy, jz, bx, by, bz, q, emf, self%u_max)
479 dbdt = -curl_up(ds,emf)
484 if (
allocated(self%heating_per_unit_mass))
then 485 dedt = dedt + d*self%heating_per_unit_mass
487 if (
allocated(self%heating_per_unit_volume))
then 488 dedt = dedt + self%heating_per_unit_volume
493 if (self%gamma /= 1d0)
then 494 u_max = self%cdtd*dsmax*self%fmaxval(abs((dedt-ee*dddt)/d),outer=.true.)
495 self%u_max = max(self%u_max,u_max)
496 if (verbose>2) print 1,
'u_max(dedt)', u_max, self%u_max
499 call deallocate_vectors_a (ld, dd, emf, j, fd, fx, fy, fz, fm)
500 call deallocate_scalars_a (lnd, pa, du, ee, cs, pp, pb, q, uu, divu)
501 if (.not.self%do_pebbles)
then 502 call deallocate_vectors (u)
503 call deallocate_scalars (pg)
505 call trace%end (itimer)
506 if (mpi%master .and. first_time)
then 507 print *, stagger%count,
' stagger calls', stagger%flops/product(
real(self%n)),
' flops per interior point' 517 SUBROUTINE pre_update (self)
519 real,
dimension(:,:,:),
pointer:: d, e, lnd, ee
524 call trace%begin(
'mhd_t%pre_update')
525 if (do_temperature)
then 526 if (.not.
associated(self%tt)) &
527 allocate (self%tt(self%gn(1),self%gn(2),self%gn(3)))
528 allocate ( lnd(self%gn(1),self%gn(2),self%gn(3)))
529 allocate ( ee(self%gn(1),self%gn(2),self%gn(3)))
530 d => self%mem(:,:,:,self%idx%d,self%it,1)
531 e => self%mem(:,:,:,self%idx%e,self%it,1)
534 call eos%lookup_table (shape(self%tt), lnd=lnd, ee=ee, tt=self%tt)
537 write(io%output,*) self%id,
'mhd_t%pre_update: tt =', &
538 minval(self%tt), maxval(self%tt)
539 call self%aux%register (
'tt', self%tt)
544 call self%extras_t%pre_update
546 END SUBROUTINE pre_update
552 SUBROUTINE update (self)
554 integer,
save:: itimer=0
556 call trace%begin(
'mhd_t%update',itimer=itimer)
563 if (self%gamma==1d0)
then 564 self%mem(:,:,:,self%idx%e,self%it,1) = &
565 self%mem(:,:,:,self%idx%d,self%it,1)*self%csound**2
567 if (self%mhd)
call self%divB_clean
572 call self%courant_condition
574 call self%check_density(
'before timestep')
575 call timestep%update (self%id, self%iit, self%dt, self%mem, self%mesh, self%lock)
576 call self%check_density(
' after timestep')
577 call self%counter_update
579 call trace%end (itimer)
580 END SUBROUTINE update
593 SUBROUTINE divb_clean1(self)
596 real,
dimension(:,:,:,:),
pointer:: b
597 real,
dimension(:,:,:) ,
pointer:: bx, by, bz, divb
598 real,
dimension(:,:,:,:) ,
allocatable:: dndivb
599 real,
parameter:: eps=0.11
601 b => self%mem(:,:,:,self%idx%bx:self%idx%bz,self%it,1)
605 call allocate_scalars (self%gn, divb)
606 call allocate_vectors_a (self%gn, dndivb)
607 divb = ddup(bx,1) + ddup(by,2) + ddup(bz,3)
608 bx = bx + eps*dddn(divb,1)
609 by = by + eps*dddn(divb,2)
610 bz = bz + eps*dddn(divb,3)
611 call deallocate_vectors_a (dndivb)
612 call deallocate_scalars(divb)
616 function ddup(f,i)
result (g)
617 real:: f(:,:,:), g(size(f,1),size(f,2),size(f,3))
620 g = cshift(f,1,i) - f
623 function dddn(f,i)
result (g)
624 real:: f(:,:,:), g(size(f,1),size(f,2),size(f,3))
627 g = f - cshift(f,-1,i)
629 END SUBROUTINE divb_clean1
635 SUBROUTINE divb_clean2 (self)
638 real,
dimension(:,:,:),
pointer:: bx, by, bz
639 real,
dimension(:,:,:),
allocatable:: divb
640 real,
parameter:: eps=0.11
641 integer:: l(3), u(3), ix, iy, iz, i1, j1, k1
644 call trace%begin(
'mhd_t%divb_clean2', itimer=itimer)
645 bx => self%mem(:,:,:,self%idx%bx,self%it,1)
646 by => self%mem(:,:,:,self%idx%by,self%it,1)
647 bz => self%mem(:,:,:,self%idx%bz,self%it,1)
650 i1 = merge(1, 0, self%n(1) > 1)
651 j1 = merge(1, 0, self%n(2) > 1)
652 k1 = merge(1, 0, self%n(3) > 1)
654 do ix=self%mesh(1)%lo,self%mesh(1)%lb,-1
655 bx(ix ,l(2):u(2),l(3):u(3)) = &
656 bx(ix+i1,l(2):u(2),l(3):u(3)) + ( &
657 by(ix,l(2)+j1:u(2)+j1,l(3) :u(3) ) - by(ix,l(2):u(2),l(3):u(3)) + &
658 bz(ix,l(2) :u(2) ,l(3)+k1:u(3)+k1) - bz(ix,l(2):u(2),l(3):u(3)))
660 do ix=self%mesh(1)%ui,self%mesh(1)%ub-1
661 bx(ix+i1,l(2):u(2),l(3):u(3)) = &
662 bx(ix ,l(2):u(2),l(3):u(3)) - ( &
663 by(ix,l(2)+j1:u(2)+j1,l(3) :u(3) ) - by(ix,l(2):u(2),l(3):u(3)) + &
664 bz(ix,l(2) :u(2) ,l(3)+k1:u(3)+k1) - bz(ix,l(2):u(2),l(3):u(3)))
667 do iy=self%mesh(2)%lo,self%mesh(2)%lb,-1
668 by(l(1):u(1),iy ,l(3):u(3)) = &
669 by(l(1):u(1),iy+j1,l(3):u(3)) + ( &
670 bx(l(1)+i1:u(1)+i1,iy,l(3) :u(3) ) - bx(l(1):u(1),iy,l(3):u(3)) + &
671 bz(l(1) :u(1) ,iy,l(3)+k1:u(3)+k1) - bz(l(1):u(1),iy,l(3):u(3)))
673 do iy=self%mesh(2)%ui,self%mesh(2)%ub-1
674 by(l(1):u(1),iy+j1,l(3):u(3)) = &
675 by(l(1):u(1),iy ,l(3):u(3)) - ( &
676 bx(l(1)+i1:u(1)+i1,iy,l(3) :u(3) ) - bx(l(1):u(1),iy,l(3):u(3)) + &
677 bz(l(1) :u(1) ,iy,l(3)+k1:u(3)+k1) - bz(l(1):u(1),iy,l(3):u(3)))
680 do iz=self%mesh(3)%lo,self%mesh(3)%lb,-1
681 bz(l(1):u(1),l(2):u(2),iz ) = &
682 bz(l(1):u(1),l(2):u(2),iz+k1) + ( &
683 bx(l(1)+i1:u(1)+i1,l(2) :u(2) ,iz) - bx(l(1):u(1),l(2):u(2),iz) + &
684 by(l(1) :u(1) ,l(2)+j1:u(2)+j1,iz) - by(l(1):u(1),l(2):u(2),iz))
686 do iz=self%mesh(3)%ui,self%mesh(3)%ub-1
687 bz(l(1):u(1),l(2):u(2),iz+k1) = &
688 bz(l(1):u(1),l(2):u(2),iz ) - ( &
689 bx(l(1)+i1:u(1)+i1,l(2) :u(2) ,iz) - bx(l(1):u(1),l(2):u(2),iz) + &
690 by(l(1) :u(1) ,l(2)+j1:u(2)+j1,iz) - by(l(1):u(1),l(2):u(2),iz))
692 call trace%end (itimer)
693 END SUBROUTINE divb_clean2
698 SUBROUTINE output (self)
700 character(len=64):: filename
704 call trace%begin (
'mhd_t%output')
705 out_next = self%out_next
706 call self%extras_t%output
707 if (self%time >= out_next .and. io%do_output .and. io%do_legacy)
then 709 write (filename,
'(a,i4.4,"_",i4.4,".txt")') trim(io%outputname), &
711 inquire(file=filename, exist=exists)
713 open (io%data_unit, file=trim(filename), form=
'formatted', status=
'old', &
715 write (io%data_unit,
'(8f10.4)') self%nu
721 END SUBROUTINE output
724 FUNCTION gas_pressure (self)
result(pg)
726 real,
dimension(self%gn(1),self%gn(2),self%gn(3)):: pg
727 real,
dimension(:,:,:),
pointer:: d, lnd, e, ee
728 integer,
save:: nprint=1
730 call trace%begin (
'mhd_t%gas_pressure_sub')
731 d => self%mem(:,:,:,self%idx%d ,self%it,1)
732 e => self%mem(:,:,:,self%idx%e ,self%it,1)
733 call allocate_scalars(self%gn,ee)
735 if (trim(self%eos)==
'ideal')
then 736 if (self%gamma==1d0)
then 737 pg = d*self%csound**2
739 pg = d*ee*(self%gamma-1)
742 if (do_temperature)
then 743 call eos%lookup_table (shape(d), d=d, ee=ee, pg=pg, tt=self%tt)
745 call eos%lookup_table (shape(d), d=d, ee=ee, pg=pg)
747 if (verbose > 1 .or. nprint > 0)
then 749 write(io%output,*) self%id,
'lookup_table: lnd =', minval(lnd), maxval(lnd)
750 write(io%output,*) self%id,
'lookup_table: ee =', minval(ee) , maxval(ee)
751 write(io%output,*) self%id,
'lookup_table: pg =', minval(pg) , maxval(pg)
752 if (do_temperature) &
753 write(io%output,*) self%id,
'lookup_table: tt =', minval(self%tt), &
757 call deallocate_scalars(ee)
759 END FUNCTION gas_pressure
762 SUBROUTINE gas_pressure_sub (self, d, lnd, ee, pg)
764 real,
dimension(:,:,:):: pg
765 real,
dimension(:,:,:),
pointer:: d, lnd, ee
766 integer,
save:: nprint=1
768 call trace%begin (
'mhd_t%gas_pressure_sub')
769 if (trim(self%eos)==
'ideal')
then 770 if (self%gamma==1d0)
then 771 pg = d*self%csound**2
773 pg = d*ee*(self%gamma-1)
776 if (do_temperature)
then 777 call eos%lookup_table (shape(lnd), lnd=lnd, ee=ee, pg=pg, tt=self%tt)
779 call eos%lookup_table (shape(lnd), lnd=lnd, ee=ee, pg=pg)
781 if (verbose > 1 .or. nprint > 0)
then 783 write(io%output,*) self%id,
'lookup_table: lnd =', minval(lnd), maxval(lnd)
784 write(io%output,*) self%id,
'lookup_table: ee =', minval(ee) , maxval(ee)
785 write(io%output,*) self%id,
'lookup_table: pg =', minval(pg) , maxval(pg)
786 if (do_temperature) &
787 write(io%output,*) self%id,
'lookup_table: tt =', minval(self%tt), &
792 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.
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...
Module with which one can register any number of pointers to real or integer arrays, and then output the contents of the links later.
6th order stagger operators, with self-test procedure