39 character(len=10):: riemann=
'hllc' 40 real,
pointer,
dimension(:,:,:,:):: uin
41 real,
pointer,
dimension(:,:,:,:):: unew
42 real,
pointer,
dimension(:,:,:,:):: q
43 real,
pointer,
dimension(:,:,:,:,:):: qp, qm
44 real,
pointer,
dimension(:,:,:,:,:):: dq
45 real,
pointer,
dimension(:,:,:,:,:):: fl
46 real,
pointer,
dimension(:,:,:):: c, s
51 procedure:: gas_pressure
53 integer,
save:: courant_type=0
54 logical,
save:: debug=.false.
55 logical,
save:: first_time=.true.
56 logical,
save:: detailed_timer=.false.
57 logical,
save:: unsigned=.true.
58 character(len=16):: eqn=
'mhd' 64 SUBROUTINE init (self)
66 character(len=10):: riemann=
'hllc' 67 real,
save:: csound=1.
69 namelist /ramses_params/ gamma, riemann, slope_type, courant_factor, smallr, &
70 smallc, isothermal, csound, detailed_timer, unsigned, courant_type, debug
72 character(len=120):: id = &
73 '$Id: 60e649508a40479a148c99021f42dfa3fe535b77 $ solvers/ramses/hydro/hd_mod.f90' 75 call trace%begin(
'hd_mod%init')
76 call trace%print_id (id)
80 self%kind =
'ramses_hd_patch' 81 if (self%nv == 0) self%nv = 5
82 call self%idx%init (5, self%mhd)
83 call self%patch_t%init
88 read (io_unit%input, ramses_params)
89 if (mpi%rank==0)
write (*, ramses_params)
96 self%riemann = riemann
97 self%courant = courant_factor
100 self%staggered = .false.
104 self%unsigned(self%idx%d) = unsigned
105 self%unsigned(self%idx%e) = unsigned
106 self%pervolume(self%idx%px) = unsigned
107 self%pervolume(self%idx%py) = unsigned
108 self%pervolume(self%idx%pz) = unsigned
109 call self%gpatch_t%init
110 call self%extras_t%init
117 SUBROUTINE dealloc (self)
120 call trace%begin (
'hd_t%dealloc')
121 call self%extras_t%dealloc
123 END SUBROUTINE dealloc
127 SUBROUTINE update (self)
129 integer:: n(3), m, i, l(3), u(3)
130 integer,
save:: itimer=0
131 real,
dimension(:,:,:,:),
pointer:: p
132 real,
dimension(:,:,:),
pointer:: d, e, ux=>null(), uy=>null(), uz=>null()
133 real,
dimension(:,:,:),
pointer:: tmp
136 call trace%begin (
'hd_t%update', itimer=itimer)
137 if (io%smallr /= 0.0) smallr = io%smallr
138 if (io%courant /= 0.0) self%courant = io%courant
149 self%uin => self%mem(:,:,:,:,self%it,1)
150 self%unew => self%mem(:,:,:,:,self%new,1)
153 allocate (self%q (n(1),n(2),n(3),self%nv))
154 allocate (self%qp(n(1),n(2),n(3),self%nv,ndim))
155 allocate (self%qm(n(1),n(2),n(3),self%nv,ndim))
156 allocate (self%dq(n(1),n(2),n(3),self%nv,ndim))
157 allocate (self%fl(n(1),n(2),n(3),self%nv,ndim))
158 allocate (self%c(n(1),n(2),n(3)))
159 allocate (self%s(n(1),n(2),n(3)))
161 if (detailed_timer)
call trace_end (itimer)
163 if (self%id==229.or.self%id==231)
then 164 print *,self%id,
'mk1:px', self%mem(19-2:19+2,19,19,self%idx%px,self%it,1), self%idx%px
165 print *,self%id,
'mk1: d', self%mem(19-2:19+2,19,19,self%idx%d ,self%it,1)
166 print *,self%id,
'mk1: e', self%mem(19-2:19+2,19,19,self%idx%e ,self%it,1), self%idx%e
168 if (self%id==230)
then 169 print *,self%id,
'mk1:px', self%mem( m-7:m-3 ,19,19,self%idx%px,self%it,1), self%idx%px
170 print *,self%id,
'mk1: d', self%mem( m-7:m-3 ,19,19,self%idx%d ,self%it,1)
171 print *,self%id,
'mk1: e', self%mem( m-7:m-3 ,19,19,self%idx%e ,self%it,1), self%idx%e
174 d => self%mem(:,:,:, self%idx%d ,self%it,1)
175 e => self%mem(:,:,:, self%idx%e ,self%it,1)
176 p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%it,1)
182 self%mem(:,:,:,self%idx%e,self%new,1) = &
183 self%mem(:,:,:,self%idx%d,self%new,1)*self%csound**2
184 if (detailed_timer)
call trace_begin (
'hd_t%update', itimer=itimer)
188 p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%new,1)
189 if (
allocated(self%force_per_unit_mass))
then 191 p(:,:,:,i) = p(:,:,:,i) + self%dtime*self%force_per_unit_mass(:,:,:,i)*d
194 if (
allocated(self%force_per_unit_volume))
then 195 p = p + self%dtime*self%force_per_unit_volume
197 call validate%check (self%link, p,
'after force')
198 if (omp_lock%tasks)
then 199 call self%lock%unset (
'update')
204 call self%counter_update
205 deallocate (self%q,self%qp,self%qm,self%dq,self%fl,self%c,self%s)
206 call trace%end (itimer)
207 END SUBROUTINE update
212 SUBROUTINE p2u(self, U, it)
214 real,
dimension(:,:,:,:),
pointer:: u
219 associate(d => self%mem(:,:,:,self%idx%d, it,1), &
220 p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
222 u(:,:,:,i) = p(:,:,:,i)/d
230 SUBROUTINE u2p(self, U, it)
232 real,
dimension(:,:,:,:),
pointer:: u
237 associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
238 p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
240 p(:,:,:,i)= u(:,:,:,i)*d
248 SUBROUTINE e2e_th(self, E_th, it)
250 real,
dimension(:,:,:),
pointer:: e_th
253 real,
dimension(:,:,:),
pointer:: k
256 associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
257 e => self%mem(:,:,:,self%idx%e ,it,1), &
258 p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
259 allocate(k(
size(d,1),
size(d,2),
size(d,3)))
268 END SUBROUTINE e2e_th
273 SUBROUTINE e2s(self, s, it)
275 real,
dimension(:,:,:),
pointer:: s
278 real,
dimension(:,:,:),
pointer:: k
285 associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
286 e => self%mem(:,:,:,self%idx%e ,it,1), &
287 p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
288 allocate(k(
size(d,1),
size(d,2),
size(d,3)))
294 s = log((e-k)*(gamma-1d0)/d**gamma)/(gamma-1d0)
302 SUBROUTINE s2e (self, s, it)
304 real,
dimension(:,:,:),
pointer:: s
307 real,
dimension(:,:,:),
pointer:: k
310 associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
311 e => self%mem(:,:,:,self%idx%e ,it,1), &
312 p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
317 allocate(k(
size(d,1),
size(d,2),
size(d,3)))
323 e = d**gamma*exp(s*(gamma-1d0)) + k
331 SUBROUTINE e_th2e(self, E_th, it)
334 real,
dimension(:,:,:),
pointer:: e_th
336 real,
dimension(:,:,:),
pointer:: k
339 associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
340 e => self%mem(:,:,:,self%idx%e ,it,1), &
341 p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
342 allocate(k(
size(d,1),
size(d,2),
size(d,3)))
351 END SUBROUTINE e_th2e
354 subroutine godunov (self)
364 integer::i0,j0,k0,i,j,k,idim,ivar,li(3),ui(3)
365 integer,
save:: itimer=0
366 associate(unew=>self%unew, uin=>self%uin, fl=>self%fl)
370 if (.not.detailed_timer)
call trace%begin (
'hd_t%godunov')
373 call check_small (self, uin)
378 if (detailed_timer)
call trace%begin (
'hd_t%godunov', itimer=itimer)
379 if (omp_lock%tasks)
then 380 call self%lock%set (
'update')
385 i0 = merge(1,0,idim==1)
386 j0 = merge(1,0,idim==2)
387 k0 = merge(1,0,idim==3)
392 unew(i,j,k,ivar) = unew(i,j,k,ivar)+ &
393 & (fl(i ,j ,k ,ivar,idim) &
394 & -fl(i+i0,j+j0,k+k0,ivar,idim))
400 if (detailed_timer)
then 401 call trace%end (itimer)
406 end subroutine godunov
410 subroutine check_small (self, uu)
412 real,
dimension(:,:,:,:):: uu
413 integer:: i,j,k,lo(3),uo(3),n_smalld,n_smallp
414 real(dp):: c,d,u,v,w,p,ekin
415 integer,
save:: nprint=10
416 integer,
save:: itimer=0
420 if (detailed_timer)
call trace%begin (
'hd_t%check_small', itimer=itimer)
429 if (uu(i,j,k,self%idx%d) < smallr)
then 433 write (io_unit%log,
'(i7,3i5,1p,2e12.3)') self%id,i,j,k,uu(i,j,k,self%idx%d),smallr
435 n_smalld = n_smalld+1
436 uu(i,j,k,self%idx%d) = smallr
438 d = uu(i,j,k,self%idx%d)
440 u = uu(i,j,k,self%idx%px)*c
441 v = uu(i,j,k,self%idx%py)*c
442 w = uu(i,j,k,self%idx%pz)*c
443 ekin = 0.5*d*(u**2+v**2+w**2)
445 uu(i,j,k,self%idx%e) = max(uu(i,j,k,self%idx%e), smallr*smallc**2)
447 p = (uu(i,j,k,self%idx%e)-ekin)*(gamma-one)
448 if (p < smallr*smallc**2)
then 450 uu(i,j,k,self%idx%e) = p/(gamma-one)+ekin
451 n_smallp = n_smallp+1
457 if (n_smalld+n_smallp > 0)
then 458 write (io_unit%log,*)
'small: id,nr,np', self%id,n_smalld,n_smallp
460 if (panic .and. self%n_dump > 0)
then 462 self%n_dump = self%n_dump-1
463 write (io_unit%dump) self%id, self%gn, self%nv, self%nt, self%nw, self%istep
464 write (io_unit%dump) self%time, self%dtime
469 write (io_unit%dump) self%mem(:,:,:,:,jt,iw)
472 write (io_unit%mpi,*)
'density < SMALLR, dumped', self%id, self%time
475 if (detailed_timer)
call trace%end (itimer)
476 end subroutine check_small
482 subroutine unsplit(self)
486 call trace%begin (
'hd_t%unsplit')
488 call self%courant_condition (detailed_timer)
493 end subroutine unsplit
499 subroutine ctoprim (self)
501 real(dp):: smallp, dtxhalf
502 integer,
save:: itimer=0
504 real,
dimension(:,:,:),
pointer:: d, e, pg
505 associate(uin=>self%uin, q=>self%q, c=>self%c, s=>self%s)
507 if (detailed_timer)
then 508 call trace%begin (
'hd_mod::ctoprim', itimer=itimer)
510 call trace%begin (
'hd_mod::ctoprim')
512 smallp = smallr*smallc**2/gamma
513 q(:,:,:,1) = max(uin(:,:,:,self%idx%d),smallr)
515 q(:,:,:,2) = uin(:,:,:,self%idx%px)*s
516 q(:,:,:,3) = uin(:,:,:,self%idx%py)*s
517 q(:,:,:,4) = uin(:,:,:,self%idx%pz)*s
518 s = q(:,:,:,2)**2+q(:,:,:,3)**2+q(:,:,:,4)**2
522 s = max((uin(:,:,:,self%idx%e)-half*q(:,:,:,1)*s(:,:,:)),smallp)
527 pg = d*self%csound**2
533 if (self%do_pebbles)
then 534 self%vgas(:,:,:,1:3,self%it) = q(:,:,:,2:4)
535 self%pressure(:,:,:, self%it) = pg
537 c(:,:,:)=sqrt(gamma*q(:,:,:,5)/q(:,:,:,1))
538 if (courant_type==1)
then 539 c = c + max(abs(q(:,:,:,2)), abs(q(:,:,:,3)), abs(q(:,:,:,4)))
541 c = c + (abs(q(:,:,:,2)) + abs(q(:,:,:,3)) + abs(q(:,:,:,4)))/3.0
545 self%u_max = maxval(c(l(1):u(1),l(2):u(2),l(3):u(3)))
546 if (detailed_timer)
then 547 call trace%end (itimer)
552 end subroutine ctoprim
557 subroutine uslope(self)
560 integer::i, j, k, n, l(3), u(3), islope_type
561 real(dp)::dsgn, dlim, dcen, dlft, drgt, slop, c1, c2
562 real(dp)::dfll,dflm,dflr,dfml,dfmm,dfmr,dfrl,dfrm,dfrr
563 real(dp)::dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl
564 real(dp)::dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm
565 real(dp)::dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr
566 real(dp)::vmin,vmax,dfx,dfy,dfz,dff,idff,xslope_type
567 integer,
save:: itimer=0
568 associate(q=>self%q, dq=>self%dq)
570 if (detailed_timer)
then 571 call trace%begin (
'hd_mod::uslope', itimer=itimer)
573 call trace%begin (
'hd_mod::uslope')
576 islope_type = int(slope_type)
579 if(islope_type==-1)
then 584 dq(i,j,k,n,1) = half*(q(i+1,j,k,n) - q(i-1,j,k,n))
585 dq(i,j,k,n,2) = half*(q(i,j+1,k,n) - q(i,j-1,k,n))
586 dq(i,j,k,n,3) = half*(q(i,j,k+1,n) - q(i,j,k-1,n))
591 else if(islope_type==0)
then 593 if (detailed_timer)
then 594 call trace%end (itimer)
599 else if (islope_type==1)
then 605 dlft = q(i ,j,k,n) - q(i-1,j,k,n)
606 drgt = q(i+1,j,k,n) - q(i ,j,k,n)
607 dq(i,j,k,n,1) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
609 dlft = q(i,j ,k,n) - q(i,j-1,k,n)
610 drgt = q(i,j+1,k,n) - q(i,j ,k,n)
611 dq(i,j,k,n,2) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
613 dlft = q(i,j,k ,n) - q(i,j,k-1,n)
614 drgt = q(i,j,k+1,n) - q(i,j,k ,n)
615 dq(i,j,k,n,3) = merge(zero,merge(min(dlft,drgt),max(dlft,drgt),dlft>0),(dlft*drgt)<=zero)
620 else if (islope_type==2)
then 628 dlft = c1*(q(i ,j,k,n) - q(i-1,j,k,n))
629 drgt = c1*(q(i+1,j,k,n) - q(i ,j,k,n))
630 dcen = c2*(dlft+drgt)
631 dsgn = sign(one, dcen)
632 slop = min(abs(dlft),abs(drgt))
633 dlim = merge(zero,slop,(dlft*drgt)<=zero)
634 dq(i,j,k,n,1) = dsgn*min(dlim,abs(dcen))
638 dlft = c1*(q(i,j ,k,n) - q(i,j-1,k,n))
639 drgt = c1*(q(i,j+1,k,n) - q(i,j ,k,n))
640 dcen = c2*(dlft+drgt)
641 dsgn = sign(one,dcen)
642 slop = min(abs(dlft),abs(drgt))
643 dlim = merge(zero,slop,(dlft*drgt)<=zero)
644 dq(i,j,k,n,2) = dsgn*min(dlim,abs(dcen))
648 dlft = c1*(q(i,j,k ,n) - q(i,j,k-1,n))
649 drgt = c1*(q(i,j,k+1,n) - q(i,j,k ,n))
650 dcen = c2*(dlft+drgt)
651 dsgn = sign(one,dcen)
652 slop = min(abs(dlft),abs(drgt))
653 dlim = merge(zero,slop,(dlft*drgt)<=zero)
654 dq(i,j,k,n,3) = dsgn*min(dlim,abs(dcen))
659 else if (islope_type==3)
then 660 xslope_type =
real(slope_type - 2_dp,kind=dp)
665 dflll = q(i-1,j-1,k-1,n)-q(i,j,k,n)
666 dflml = q(i-1,j ,k-1,n)-q(i,j,k,n)
667 dflrl = q(i-1,j+1,k-1,n)-q(i,j,k,n)
668 dfmll = q(i ,j-1,k-1,n)-q(i,j,k,n)
669 dfmml = q(i ,j ,k-1,n)-q(i,j,k,n)
670 dfmrl = q(i ,j+1,k-1,n)-q(i,j,k,n)
671 dfrll = q(i+1,j-1,k-1,n)-q(i,j,k,n)
672 dfrml = q(i+1,j ,k-1,n)-q(i,j,k,n)
673 dfrrl = q(i+1,j+1,k-1,n)-q(i,j,k,n)
675 dfllm = q(i-1,j-1,k ,n)-q(i,j,k,n)
676 dflmm = q(i-1,j ,k ,n)-q(i,j,k,n)
677 dflrm = q(i-1,j+1,k ,n)-q(i,j,k,n)
678 dfmlm = q(i ,j-1,k ,n)-q(i,j,k,n)
679 dfmmm = q(i ,j ,k ,n)-q(i,j,k,n)
680 dfmrm = q(i ,j+1,k ,n)-q(i,j,k,n)
681 dfrlm = q(i+1,j-1,k ,n)-q(i,j,k,n)
682 dfrmm = q(i+1,j ,k ,n)-q(i,j,k,n)
683 dfrrm = q(i+1,j+1,k ,n)-q(i,j,k,n)
685 dfllr = q(i-1,j-1,k+1,n)-q(i,j,k,n)
686 dflmr = q(i-1,j ,k+1,n)-q(i,j,k,n)
687 dflrr = q(i-1,j+1,k+1,n)-q(i,j,k,n)
688 dfmlr = q(i ,j-1,k+1,n)-q(i,j,k,n)
689 dfmmr = q(i ,j ,k+1,n)-q(i,j,k,n)
690 dfmrr = q(i ,j+1,k+1,n)-q(i,j,k,n)
691 dfrlr = q(i+1,j-1,k+1,n)-q(i,j,k,n)
692 dfrmr = q(i+1,j ,k+1,n)-q(i,j,k,n)
693 dfrrr = q(i+1,j+1,k+1,n)-q(i,j,k,n)
695 vmin = min(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
696 & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
697 & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
698 vmax = max(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
699 & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
700 & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
702 dfx = half*(q(i+1,j,k,n)-q(i-1,j,k,n))
703 dfy = half*(q(i,j+1,k,n)-q(i,j-1,k,n))
704 dfz = half*(q(i,j,k+1,n)-q(i,j,k-1,n))
705 idff = xslope_type / sqrt(dfx**2+dfy**2+dfz**2+tiny(1.0_dp))
716 slop = merge(min(one,min(abs(vmin),abs(vmax))*idff), one, idff<0.001_dp * huge(1.0_dp))
718 dq(i,j,k,n,1) = dlim*dfx
719 dq(i,j,k,n,2) = dlim*dfy
720 dq(i,j,k,n,3) = dlim*dfz
725 else if (islope_type==4)
then 726 xslope_type =
real(slope_type - 3_dp,kind=dp)
731 real(dp):: adfx(self%nv), adfy(self%nv), adfz(self%nv)
734 dflll = q(i-1,j-1,k-1,n)-q(i,j,k,n)
735 dflml = q(i-1,j ,k-1,n)-q(i,j,k,n)
736 dflrl = q(i-1,j+1,k-1,n)-q(i,j,k,n)
737 dfmll = q(i ,j-1,k-1,n)-q(i,j,k,n)
738 dfmml = q(i ,j ,k-1,n)-q(i,j,k,n)
739 dfmrl = q(i ,j+1,k-1,n)-q(i,j,k,n)
740 dfrll = q(i+1,j-1,k-1,n)-q(i,j,k,n)
741 dfrml = q(i+1,j ,k-1,n)-q(i,j,k,n)
742 dfrrl = q(i+1,j+1,k-1,n)-q(i,j,k,n)
744 dfllm = q(i-1,j-1,k ,n)-q(i,j,k,n)
745 dflmm = q(i-1,j ,k ,n)-q(i,j,k,n)
746 dflrm = q(i-1,j+1,k ,n)-q(i,j,k,n)
747 dfmlm = q(i ,j-1,k ,n)-q(i,j,k,n)
748 dfmmm = q(i ,j ,k ,n)-q(i,j,k,n)
749 dfmrm = q(i ,j+1,k ,n)-q(i,j,k,n)
750 dfrlm = q(i+1,j-1,k ,n)-q(i,j,k,n)
751 dfrmm = q(i+1,j ,k ,n)-q(i,j,k,n)
752 dfrrm = q(i+1,j+1,k ,n)-q(i,j,k,n)
754 dfllr = q(i-1,j-1,k+1,n)-q(i,j,k,n)
755 dflmr = q(i-1,j ,k+1,n)-q(i,j,k,n)
756 dflrr = q(i-1,j+1,k+1,n)-q(i,j,k,n)
757 dfmlr = q(i ,j-1,k+1,n)-q(i,j,k,n)
758 dfmmr = q(i ,j ,k+1,n)-q(i,j,k,n)
759 dfmrr = q(i ,j+1,k+1,n)-q(i,j,k,n)
760 dfrlr = q(i+1,j-1,k+1,n)-q(i,j,k,n)
761 dfrmr = q(i+1,j ,k+1,n)-q(i,j,k,n)
762 dfrrr = q(i+1,j+1,k+1,n)-q(i,j,k,n)
764 vmin = min(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
765 & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
766 & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
767 vmax = max(dflll,dflml,dflrl,dfmll,dfmml,dfmrl,dfrll,dfrml,dfrrl, &
768 & dfllm,dflmm,dflrm,dfmlm,dfmmm,dfmrm,dfrlm,dfrmm,dfrrm, &
769 & dfllr,dflmr,dflrr,dfmlr,dfmmr,dfmrr,dfrlr,dfrmr,dfrrr)
771 dfx = half*(q(i+1,j,k,n)-q(i-1,j,k,n))
772 dfy = half*(q(i,j+1,k,n)-q(i,j-1,k,n))
773 dfz = half*(q(i,j,k+1,n)-q(i,j,k-1,n))
774 idff = xslope_type / sqrt(dfx**2+dfy**2+dfz**2+tiny(1.0_dp))
786 slop = merge(min(one,min(abs(vmin),abs(vmax))*idff), one, idff<0.001_dp * huge(1.0_dp))
792 dlim = min(dlim,slop)
800 dq(i,j,k,n,1) = dlim*adfx(n)
801 dq(i,j,k,n,2) = dlim*adfy(n)
802 dq(i,j,k,n,3) = dlim*adfz(n)
809 write(*,*)
'Unknown slope type' 812 if (detailed_timer)
then 813 call trace%end (itimer)
818 end subroutine uslope
823 subroutine trace3d(self)
825 integer ::i, j, k, lb(3), ub(3)
826 integer ::ir=1, iu=2, iv=3, iw=4, ip=5
827 real(dp)::dtdx, dtdy, dtdz
828 real(dp)::r, u, v, w, p, a
829 real(dp)::drx, dux, dvx, dwx, dpx, dax
830 real(dp)::dry, duy, dvy, dwy, dpy, day
831 real(dp)::drz, duz, dvz, dwz, dpz, daz
832 real(dp)::sr0, su0, sv0, sw0, sp0, sa0
833 real(dp)::ff(self%gn(1),3)
834 integer,
save:: itimer=0
835 associate(q=>self%q, dq=>self%dq, qp=>self%qp, qm=>self%qm)
837 if (detailed_timer)
then 838 call trace%begin (
'hd_mod::trace3d', itimer=itimer)
840 call trace%begin (
'hd_mod::trace3d')
842 dtdx = self%dtime/self%mesh(1)%d
843 dtdy = self%dtime/self%mesh(2)%d
844 dtdz = self%dtime/self%mesh(3)%d
850 if (
allocated(self%force_per_unit_mass))
then 852 ff(i,1) = self%force_per_unit_mass(i,j,k,1)*self%mesh(1)%d
853 ff(i,2) = self%force_per_unit_mass(i,j,k,2)*self%mesh(2)%d
854 ff(i,3) = self%force_per_unit_mass(i,j,k,3)*self%mesh(3)%d
887 sr0 = -u*drx-v*dry-w*drz - (dux+dvy+dwz)*r
888 sp0 = -u*dpx-v*dpy-w*dpz - (dux+dvy+dwz)*gamma*p
889 su0 = -u*dux-v*duy-w*duz - dpx/r + ff(i,1)
890 sv0 = -u*dvx-v*dvy-w*dvz - dpy/r + ff(i,2)
891 sw0 = -u*dwx-v*dwy-w*dwz - dpz/r + ff(i,3)
894 qp(i,j,k,ir,1) = r - half*drx + sr0*dtdx*half
895 qp(i,j,k,ip,1) = p - half*dpx + sp0*dtdx*half
896 qp(i,j,k,iu,1) = u - half*dux + su0*dtdx*half
897 qp(i,j,k,iv,1) = v - half*dvx + sv0*dtdx*half
898 qp(i,j,k,iw,1) = w - half*dwx + sw0*dtdx*half
899 qp(i,j,k,ir,1) = max(smallr, qp(i,j,k,ir,1))
902 qm(i+1,j,k,ir,1) = r + half*drx + sr0*dtdx*half
903 qm(i+1,j,k,ip,1) = p + half*dpx + sp0*dtdx*half
904 qm(i+1,j,k,iu,1) = u + half*dux + su0*dtdx*half
905 qm(i+1,j,k,iv,1) = v + half*dvx + sv0*dtdx*half
906 qm(i+1,j,k,iw,1) = w + half*dwx + sw0*dtdx*half
907 qm(i+1,j,k,ir,1) = max(smallr, qm(i+1,j,k,ir,1))
910 qp(i,j,k,ir,2) = r - half*dry + sr0*dtdy*half
911 qp(i,j,k,ip,2) = p - half*dpy + sp0*dtdy*half
912 qp(i,j,k,iu,2) = u - half*duy + su0*dtdy*half
913 qp(i,j,k,iv,2) = v - half*dvy + sv0*dtdy*half
914 qp(i,j,k,iw,2) = w - half*dwy + sw0*dtdy*half
915 qp(i,j,k,ir,2) = max(smallr, qp(i,j,k,ir,2))
918 qm(i,j+1,k,ir,2) = r + half*dry + sr0*dtdy*half
919 qm(i,j+1,k,ip,2) = p + half*dpy + sp0*dtdy*half
920 qm(i,j+1,k,iu,2) = u + half*duy + su0*dtdy*half
921 qm(i,j+1,k,iv,2) = v + half*dvy + sv0*dtdy*half
922 qm(i,j+1,k,iw,2) = w + half*dwy + sw0*dtdy*half
923 qm(i,j+1,k,ir,2) = max(smallr, qm(i,j+1,k,ir,2))
926 qp(i,j,k,ir,3) = r - half*drz + sr0*dtdz*half
927 qp(i,j,k,ip,3) = p - half*dpz + sp0*dtdz*half
928 qp(i,j,k,iu,3) = u - half*duz + su0*dtdz*half
929 qp(i,j,k,iv,3) = v - half*dvz + sv0*dtdz*half
930 qp(i,j,k,iw,3) = w - half*dwz + sw0*dtdz*half
931 qp(i,j,k,ir,3) = max(smallr, qp(i,j,k,ir,3))
934 qm(i,j,k+1,ir,3) = r + half*drz + sr0*dtdz*half
935 qm(i,j,k+1,ip,3) = p + half*dpz + sp0*dtdz*half
936 qm(i,j,k+1,iu,3) = u + half*duz + su0*dtdz*half
937 qm(i,j,k+1,iv,3) = v + half*dvz + sv0*dtdz*half
938 qm(i,j,k+1,iw,3) = w + half*dwz + sw0*dtdz*half
939 qm(i,j,k+1,ir,3) = max(smallr, qm(i,j,k+1,ir,3))
943 if (detailed_timer)
then 944 call trace%end (itimer)
949 end subroutine trace3d
954 subroutine cmpflxm(self)
957 integer:: i, j, k, lb(3), ub(3)
958 real(dp),
dimension(self%mesh(1)%gn,self%nv):: qleft, qright, fgdnv
961 integer,
save:: itimer=0
962 associate(qp=>self%qp, qm=>self%qm, fl=>self%fl)
964 if (detailed_timer)
then 965 call trace%begin (
'hd_mod::cmpflxm', itimer=itimer)
967 call trace%begin (
'hd_mod::cmpflxm')
970 dtds = self%dtime/self%mesh%d
973 qleft=0.0; qright=0.0
977 qleft(i,1) = qm(i,j,k,1,1); qright(i,1) = qp(i,j,k,1,1)
978 qleft(i,2) = qm(i,j,k,2,1); qright(i,2) = qp(i,j,k,2,1)
979 qleft(i,3) = qm(i,j,k,5,1); qright(i,3) = qp(i,j,k,5,1)
980 qleft(i,4) = qm(i,j,k,3,1); qright(i,4) = qp(i,j,k,3,1)
981 qleft(i,5) = qm(i,j,k,4,1); qright(i,5) = qp(i,j,k,4,1)
983 if (detailed_timer)
call trace_end (itimer)
984 call solver (self,self%mesh(1)%gn,qleft,qright,fgdnv,error,detailed_timer)
985 if (detailed_timer)
call trace_begin (
'hd_mod::cmpflxm', itimer=itimer)
986 if (error)
call error_handler(
'x', j, k, lb(1), ub(1))
988 fl(i,j,k,self%idx%d ,1) = fgdnv(i,1)*dtds(1)
989 fl(i,j,k,self%idx%px,1) = fgdnv(i,2)*dtds(1)
990 fl(i,j,k,self%idx%e ,1) = fgdnv(i,3)*dtds(1)
991 fl(i,j,k,self%idx%py,1) = fgdnv(i,4)*dtds(1)
992 fl(i,j,k,self%idx%pz,1) = fgdnv(i,5)*dtds(1)
996 if (any(self%mesh%gn/=self%mesh(1)%gn))
then 1003 qleft(j,1) = qm(i,j,k,1,2); qright(j,1) = qp(i,j,k,1,2)
1004 qleft(j,2) = qm(i,j,k,3,2); qright(j,2) = qp(i,j,k,3,2)
1005 qleft(j,3) = qm(i,j,k,5,2); qright(j,3) = qp(i,j,k,5,2)
1006 qleft(j,4) = qm(i,j,k,2,2); qright(j,4) = qp(i,j,k,2,2)
1007 qleft(j,5) = qm(i,j,k,4,2); qright(j,5) = qp(i,j,k,4,2)
1009 if (detailed_timer)
call trace_end (itimer)
1010 call solver (self,self%mesh(2)%gn,qleft,qright,fgdnv,error,detailed_timer)
1011 if (detailed_timer)
call trace_begin (
'hd_mod::cmpflxm', itimer=itimer)
1012 if (error)
call error_handler(
'y', i, k, lb(2), ub(2))
1014 fl(i,j,k,self%idx%d ,2) = fgdnv(j,1)*dtds(2)
1015 fl(i,j,k,self%idx%py,2) = fgdnv(j,2)*dtds(2)
1016 fl(i,j,k,self%idx%e ,2) = fgdnv(j,3)*dtds(2)
1017 fl(i,j,k,self%idx%px,2) = fgdnv(j,4)*dtds(2)
1018 fl(i,j,k,self%idx%pz,2) = fgdnv(j,5)*dtds(2)
1025 qleft(k,1) = qm(i,j,k,1,3); qright(k,1) = qp(i,j,k,1,3)
1026 qleft(k,2) = qm(i,j,k,4,3); qright(k,2) = qp(i,j,k,4,3)
1027 qleft(k,3) = qm(i,j,k,5,3); qright(k,3) = qp(i,j,k,5,3)
1028 qleft(k,4) = qm(i,j,k,2,3); qright(k,4) = qp(i,j,k,2,3)
1029 qleft(k,5) = qm(i,j,k,3,3); qright(k,5) = qp(i,j,k,3,3)
1031 if (detailed_timer)
call trace_end (itimer)
1032 call solver (self,self%mesh(3)%gn,qleft,qright,fgdnv,error,detailed_timer)
1033 if (detailed_timer)
call trace_begin (
'hd_mod::cmpflxm', itimer=itimer)
1034 if (error)
call error_handler(
'z', i, j, lb(3), ub(3))
1036 fl(i,j,k,self%idx%d ,3) = fgdnv(k,1)*dtds(3)
1037 fl(i,j,k,self%idx%pz,3) = fgdnv(k,2)*dtds(3)
1038 fl(i,j,k,self%idx%e ,3) = fgdnv(k,3)*dtds(3)
1039 fl(i,j,k,self%idx%px,3) = fgdnv(k,4)*dtds(3)
1040 fl(i,j,k,self%idx%py,3) = fgdnv(k,5)*dtds(3)
1046 if (detailed_timer)
then 1047 call trace%end (itimer)
1052 subroutine cmpflxm_y
1054 real(dp),
dimension(self%mesh(2)%gn,self%nv):: qleft, qright, fgdnv
1055 associate(qp=>self%qp, qm=>self%qm, fl=>self%fl)
1058 dtds = self%dtime/self%mesh%d
1061 qleft=0.0; qright=0.0
1064 qleft(:,1) = qm(i,:,k,1,2); qright(:,1) = qp(i,:,k,1,2)
1065 qleft(:,2) = qm(i,:,k,3,2); qright(:,2) = qp(i,:,k,3,2)
1066 qleft(:,3) = qm(i,:,k,5,2); qright(:,3) = qp(i,:,k,5,2)
1067 qleft(:,4) = qm(i,:,k,2,2); qright(:,4) = qp(i,:,k,2,2)
1068 qleft(:,5) = qm(i,:,k,4,2); qright(:,5) = qp(i,:,k,4,2)
1069 if (detailed_timer)
call trace_end (itimer)
1070 call solver (self,self%mesh(2)%gn,qleft,qright,fgdnv,error,detailed_timer)
1071 if (error)
call error_handler(
'y', i, k, lb(2), ub(2))
1072 if (detailed_timer)
call trace_begin (
'hd_mod::cmpflxm', itimer=itimer)
1073 fl(i,:,k,self%idx%d ,2) = fgdnv(:,1)*dtds(2)
1074 fl(i,:,k,self%idx%py,2) = fgdnv(:,2)*dtds(2)
1075 fl(i,:,k,self%idx%e ,2) = fgdnv(:,3)*dtds(2)
1076 fl(i,:,k,self%idx%px,2) = fgdnv(:,4)*dtds(2)
1077 fl(i,:,k,self%idx%pz,2) = fgdnv(:,5)*dtds(2)
1081 end subroutine cmpflxm_y
1084 subroutine cmpflxm_z
1086 real(dp),
dimension(self%mesh(3)%gn,self%nv):: qleft, qright, fgdnv
1087 associate(qp=>self%qp, qm=>self%qm, fl=>self%fl)
1089 dtds = self%dtime/self%mesh%d
1092 qleft=0.0; qright=0.0
1095 qleft(:,1) = qm(i,j,:,1,3); qright(:,1) = qp(i,j,:,1,3)
1096 qleft(:,2) = qm(i,j,:,4,3); qright(:,2) = qp(i,j,:,4,3)
1097 qleft(:,3) = qm(i,j,:,5,3); qright(:,3) = qp(i,j,:,5,3)
1098 qleft(:,4) = qm(i,j,:,2,3); qright(:,4) = qp(i,j,:,2,3)
1099 qleft(:,5) = qm(i,j,:,3,3); qright(:,5) = qp(i,j,:,3,3)
1100 if (detailed_timer)
call trace_end (itimer)
1101 call solver (self,self%mesh(3)%gn,qleft,qright,fgdnv,error,detailed_timer)
1102 if (error)
call error_handler(
'z', i, j, lb(3), ub(3))
1103 if (detailed_timer)
call trace_begin (
'hd_mod::cmpflxm', itimer=itimer)
1104 fl(i,j,:,self%idx%d ,3) = fgdnv(:,1)*dtds(3)
1105 fl(i,j,:,self%idx%pz,3) = fgdnv(:,2)*dtds(3)
1106 fl(i,j,:,self%idx%e ,3) = fgdnv(:,3)*dtds(3)
1107 fl(i,j,:,self%idx%px,3) = fgdnv(:,4)*dtds(3)
1108 fl(i,j,:,self%idx%py,3) = fgdnv(:,5)*dtds(3)
1112 end subroutine cmpflxm_z
1115 subroutine error_handler (label, i1, i2, lb, ub)
1116 class(*),
pointer:: link
1117 character(len=*):: label
1118 integer:: i1, i2, lb, ub
1119 class(link_t),
pointer:: nbor
1121 write (io_unit%log,*)
'solver error: id, i[12] =', self%id, i1, i2
1126 write (io_unit%log,
'(a,1p,2(5e12.3,2x))') &
1127 trim(label)//
'solver error:', qleft(i,1:5), qright(i,1:5)
1139 do while (
associated(nbor))
1140 write (io_unit%log,
'(a,i9,f10.3,3x,3f6.2)') &
1141 trim(label)//
'solver error: nbor, time diff', &
1142 nbor%task%id, (nbor%task%time-self%time)/self%dtime, &
1143 (nbor%task%position-self%position)/self%size
1146 call mpi%abort (
'HLLC')
1148 end subroutine cmpflxm
1153 subroutine solver (self,ngrid,qleft,qright,fgdnv,error,detailed_timer)
1156 real(dp),
dimension(:,:):: qleft,qright,fgdnv
1157 logical:: error, detailed_timer
1159 if(self%riemann.eq.
'acoustic')
then 1160 call rieman%acoustic(qleft,qright,fgdnv,ngrid,self%idx%e)
1161 else if (self%riemann.eq.
'exact')
then 1162 call rieman%approx (qleft,qright,fgdnv,ngrid,self%idx%e)
1163 else if (self%riemann.eq.
'llf')
then 1164 call rieman%llf (qleft,qright,fgdnv,ngrid,self%idx%e)
1165 else if (self%riemann.eq.
'hllc')
then 1166 call rieman%hllc (qleft,qright,fgdnv,ngrid,self%idx%e,isothermal,error,detailed_timer)
1167 else if (self%riemann.eq.
'hllc_a')
then 1168 call rieman%hllc_a (qleft,qright,fgdnv,ngrid,self%idx%e)
1169 else if (self%riemann.eq.
'hllc_v')
then 1170 call rieman%hllc_v (qleft,qright,fgdnv,ngrid,self%idx%e)
1171 else if (self%riemann.eq.
'hll')
then 1172 call rieman%hll (qleft,qright,fgdnv,ngrid,self%idx%e)
1174 write(*,*)
'unknown Riemann solver' 1177 end subroutine solver
1181 FUNCTION gas_pressure (self)
RESULT (pg)
1183 real,
dimension(self%gn(1),self%gn(2),self%gn(3)):: pg
1185 associate(d => self%mem(:,:,:,self%idx%d ,self%it,1), &
1186 px => self%mem(:,:,:,self%idx%px,self%it,1), &
1187 py => self%mem(:,:,:,self%idx%py,self%it,1), &
1188 pz => self%mem(:,:,:,self%idx%pz,self%it,1), &
1189 e => self%mem(:,:,:,self%idx%e ,self%it,1))
1190 if (isothermal)
then 1191 pg = d*self%csound**2
1193 pg = (gamma-1d0)*(e-0.5*((px**2+py**2+pz**2)/d))
1196 END FUNCTION gas_pressure
1200 SUBROUTINE gravitational_force (self)
1202 real,
dimension(:,:,:),
pointer:: phi, d, et, gradphi1, gradphi2, gradphi3
1203 real,
dimension(:,:,:,:),
pointer:: p
1205 real,
dimension(3):: wtp1, wtm1
1208 associate(l => self%li, u => self%ui)
1209 wtp1 = 0.50 * 4.0 / 3.0 / self%ds
1210 wtm1 = 0.25 * 1.0 / 3.0 / self%ds
1211 phi => self%mem(:,:,:,self%idx%phi,self%it,1)
1212 d => self%mem(:,:,:,self%idx%d,self%it,1)
1213 et => self%mem(:,:,:,self%idx%e,self%new,1)
1214 p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%new,1)
1215 allocate (gradphi1(self%gn(1),self%gn(2),self%gn(3)), &
1216 gradphi2(self%gn(1),self%gn(2),self%gn(3)), &
1217 gradphi3(self%gn(1),self%gn(2),self%gn(3)))
1223 etemp = et(i,j,k) - 0.5 * p(i,j,k,1)**2 / d(i,j,k) &
1224 - 0.5 * p(i,j,k,2)**2 / d(i,j,k) &
1225 - 0.5 * p(i,j,k,3)**2 / d(i,j,k)
1227 gradphi1(i,j,k) = wtp1(1) * (phi(i+1,j,k) - phi(i-1,j,k)) &
1228 - wtm1(1) * (phi(i+2,j,k) - phi(i-2,j,k))
1229 gradphi2(i,j,k) = wtp1(2) * (phi(i,j+1,k) - phi(i,j-1,k)) &
1230 - wtm1(2) * (phi(i,j+2,k) - phi(i,j-2,k))
1231 gradphi3(i,j,k) = wtp1(3) * (phi(i,j,k+1) - phi(i,j,k-1)) &
1232 - wtm1(3) * (phi(i,j,k+2) - phi(i,j,k-2))
1234 p(i,j,k,1) = p(i,j,k,1) - self%dtime * d(i,j,k) * gradphi1(i,j,k)
1235 p(i,j,k,2) = p(i,j,k,2) - self%dtime * d(i,j,k) * gradphi2(i,j,k)
1236 p(i,j,k,3) = p(i,j,k,3) - self%dtime * d(i,j,k) * gradphi3(i,j,k)
1238 et(i,j,k) = etemp + 0.5 * p(i,j,k,1)**2 / d(i,j,k) &
1239 + 0.5 * p(i,j,k,2)**2 / d(i,j,k) &
1240 + 0.5 * p(i,j,k,3)**2 / d(i,j,k)
1245 deallocate (gradphi1, gradphi2, gradphi3)
1247 END SUBROUTINE gravitational_force
Generic validation module. The general idea is to be able to compare two runs at critical points in t...
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Equation of state module for any sort of tables, provided by a reader.
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...