MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_random_heating.t
Go to the documentation of this file.
3  use mod_comm_lib, only: mpistop
4  double precision, allocatable :: va1(:),ta1(:),tb1(:),va2(:),ta2(:),tb2(:)
5  !> This module is for imposing randomized heating at x-axis for 2.5D simulation, for details please see Zhou et al. (2020, Nature Astronomy,4,994) and
6  !> Li et al. (2022,ApJ,926,216).
7  !> Users can define the parameters "trelax,vtim,vlim,periods,variation,num,mnx" in the mod_usr.t file.
8  !> Example for using this module please see tests/mhd/coronal_rain_2.5D
9  !trelax is the relaxization time for the simulation, if the relaxization is done seperately, then trelax=0.d0
10  !periods=300,variation=75 !The time durations of each heating episode is typically 300s +/- 75s
11  !vtim = 100 !100 heating episodes, e.g., 500 min. Enough for the simulation
12  !num = 1000 !1000 waves for each episode
13  !si=5./6. !spectral index
14  !mnx could be chosen as twice of the spatial resolution in the simulation
15  !vlim could be chosen as the same with mnx
16 
17 contains
18 
19  subroutine getlqgrid(lQgrid,ixI^L,ixO^L,qt,trelax,x,mode,vtim,vlim)
21  integer, intent(in) :: ixI^L, ixO^L, mode, vtim, vlim
22  double precision, intent(in) :: qt, trelax, x(ixI^S,1:ndim)
23  double precision, intent(inout) :: lQgrid(ixI^S)
24 
25  double precision :: tl1,tl2,tl3,tl4,tl5,tl6,t2,stt
26  double precision :: tr1,tr2,tr3,tr4,tr5,tr6
27  double precision :: va0
28  integer :: Bp,Bt1,Bt2,ix^D,i
29 
30  t2 = qt - trelax
31  va0=0.4d0
32  stt = (xprobmax1- xprobmin1) / (dble(vlim) - 1.d0)
33 
34  select case(mode)
35  ! mode 1 for gaussian, mode 2 for sinc
36  case(1)
37  if(qt .lt. trelax) then
38  lqgrid(ixo^s)=0.d0
39  else
40  do i = 4, vtim
41  if(qt .lt. trelax + tb1(i)) then
42  tl1 = dexp(-(t2 - tb1(i-3))**2/(0.5d0*ta1(i-3))**2)
43  tl2 = dexp(-(t2 - tb1(i-2))**2/(0.5d0*ta1(i-2))**2)
44  tl3 = dexp(-(t2 - tb1(i-1))**2/(0.5d0*ta1(i-1))**2)
45  tl4 = dexp(-(t2 - tb1(i))**2/(0.5d0*ta1(i))**2)
46  tl5 = dexp(-(t2 - tb1(i+1))**2/(0.5d0*ta1(i+1))**2)
47  tl6 = dexp(-(t2 - tb1(i+2))**2/(0.5d0*ta1(i+2))**2)
48  {do ix^db = ixomin^db,ixomax^db\}
49  bp = floor((x(ix^d,1) - xprobmin1) / stt) + 1
50  if (bp .lt. 1) bp=1
51  if(x(ix^d,1) .lt. 0.d0) then
52  lqgrid(ix^d) = tl1 * (va1(bp + (i-4)*vlim)*(one-va0)+va0) &
53  + tl2 * (va1(bp + (i-3)*vlim)*(one-va0)+va0) &
54  + tl3 * (va1(bp + (i-2)*vlim)*(one-va0)+va0) &
55  + tl4 * (va1(bp + (i-1)*vlim)*(one-va0)+va0) &
56  + tl5 * (va1(bp + (i)*vlim)*(one-va0)+va0) &
57  + tl6 * (va1(bp + (i+1)*vlim)*(one-va0)+va0)
58  endif
59  {enddo\}
60  exit
61  endif
62  enddo
63 
64  do i = 4, vtim
65  if(qt .lt. trelax + tb2(i)) then
66  tr1 = dexp(-(t2 - tb2(i-3))**2/(0.5d0*ta2(i-3))**2)
67  tr2 = dexp(-(t2 - tb2(i-2))**2/(0.5d0*ta2(i-2))**2)
68  tr3 = dexp(-(t2 - tb2(i-1))**2/(0.5d0*ta2(i-1))**2)
69  tr4 = dexp(-(t2 - tb2(i))**2/(0.5d0*ta2(i))**2)
70  tr5 = dexp(-(t2 - tb2(i+1))**2/(0.5d0*ta2(i+1))**2)
71  tr6 = dexp(-(t2 - tb2(i+2))**2/(0.5d0*ta2(i+2))**2)
72  {do ix^db = ixomin^db,ixomax^db\}
73  bp = floor((x(ix^d,1) - xprobmin1) / stt) + 1
74  if (bp .lt. 1) bp=1
75  if(x(ix^d,1) .ge. 0.d0) then
76  lqgrid(ix^d) = tr1 * (va2(bp + (i-4)*vlim)*(one-va0)+va0) &
77  + tr2 * (va2(bp + (i-3)*vlim)*(one-va0)+va0) &
78  + tr3 * (va2(bp + (i-2)*vlim)*(one-va0)+va0) &
79  + tr4 * (va2(bp + (i-1)*vlim)*(one-va0)+va0) &
80  + tr5 * (va2(bp + (i)*vlim)*(one-va0)+va0) &
81  + tr6 * (va2(bp + (i+1)*vlim)*(one-va0)+va0)
82  endif
83  {enddo\}
84  exit
85  endif
86  enddo
87  endif
88  case(2)
89  if(qt .lt. trelax) then
90  tr1 = 0.d0
91  tr3 = 0.d0
92  bt1 = 0
93  else if(qt .lt. trelax + tb1(1)) then
94  tr1 = dsin(t2 * dpi / tb1(2))
95  tr3 = 0.d0
96  bt1 = 0
97  else if(qt .lt. trelax + tb1(2)) then
98  tr1 = dsin(t2 * dpi / tb1(2))
99  tr3 = dsin((t2 - tb1(1)) * dpi / (tb1(3) - tb1(1)))
100  bt1 = 0
101  else
102  do i = 2, vtim
103  if(t2 .lt. tb1(i)) then
104  bt1 = i - 1
105  exit
106  endif
107  enddo
108  tr1 = dsin((t2 - tb1(bt1-1)) * dpi / (tb1(bt1+1) - tb1(bt1-1)))
109  tr3 = dsin((t2 - tb1(bt1)) * dpi / (tb1(bt1+2) - tb1(bt1)))
110  endif
111  if(qt .lt. trelax) then
112  tr2 = 0.d0
113  tr4 = 0.d0
114  bt2 = 0
115  else if(qt .lt. trelax + tb2(1)) then
116  tr2 = dsin(t2 * dpi / tb2(2))
117  tr4 = 0.d0
118  bt2 = 0
119  else if(qt .lt. trelax + tb2(2)) then
120  tr1 = dsin(t2 * dpi / tb2(2))
121  tr3 = dsin((t2 - tb2(1)) * dpi / (tb2(3) - tb2(1)))
122  bt2 = 0
123  else
124  do i = 2, vtim
125  if(t2 .lt. tb2(i)) then
126  bt2 = i - 1
127  exit
128  endif
129  enddo
130  tr2 = dsin((t2 - tb2(bt2-1)) * dpi / (tb2(bt2+1) - tb2(bt2-1)))
131  tr4 = dsin((t2 - tb2(bt2)) * dpi / (tb2(bt2+2) - tb2(bt2)))
132  endif
133  {do ix^db = ixomin^db,ixomax^db\}
134  bp = floor((x(ix^d,1) - xprobmin1) / stt) + 1
135  if(x(ix^d,1) .lt. 0.d0) then
136  lqgrid(ix^d) = tr1 * (va1(bp + bt1*vlim)*(one-va0)+va0) + tr3 * (va1(bp+(bt1+1)*vlim)*(one-va0)+va0)
137  else
138  lqgrid(ix^d) = tr2 * (va2(bp + bt2*vlim)*(one-va0)+va0) + tr4 * (va2(bp+(bt2+1)*vlim)*(one-va0)+va0)
139  endif
140  {enddo\}
141  case default
142  call mpistop('unknow mode')
143  end select
144 
145  end subroutine getlqgrid
146 
147  subroutine generatetv(vtim,vlim,periods,variation,mnx,num,si)
148  integer, intent(in) :: vtim,vlim,mnx,num
149  double precision, intent(in) :: periods,variation,si
150  integer :: i,vx,j
151  character(len=100) :: filename
152  logical :: alive
153 
154  ! generate random T, typically 300s +/- 75s
155 
156  allocate(ta1(vtim)) ! left footpoint
157  allocate(tb1(vtim))
158  allocate(ta2(vtim)) ! right footpoint
159  allocate(tb2(vtim))
160 
161  if (ndim/=2) call mpistop("Randomized heating for 2D, 2.5D simulations only!")
162  if (mype==0) then
163  write(filename,"(a)") "LeftT.dat"
164  inquire(file=filename,exist=alive)
165  if(alive) then
166  open(unit=21,file=filename,status='old')
167  do j=1,vtim
168  read(21,*) tb1(j)
169  enddo
170  close(21)
171  else
172  call randomt(tb1,vtim,periods,variation,filename)
173  endif
174 
175  write(filename,"(a)") "RightT.dat"
176  inquire(file=filename,exist=alive)
177  if(alive) then
178  open(unit=21,file=filename,status='old')
179  do j=1,vtim
180  read(21,*) tb2(j)
181  enddo
182  close(21)
183  else
184  call randomt(tb2,vtim,periods,variation,filename)
185  endif
186 
187  ta1 = tb1
188  ta2 = tb2
189  do i = 2, vtim
190  ta1(i) = tb1(i) - tb1(i-1)
191  ta2(i) = tb2(i) - tb2(i-1)
192  enddo
193  endif
194  call mpi_barrier(icomm,ierrmpi)
195  if (npe>1) then
196  call mpi_bcast(tb1,vtim,mpi_double_precision,0,icomm,ierrmpi)
197  call mpi_bcast(ta1,vtim,mpi_double_precision,0,icomm,ierrmpi)
198  call mpi_bcast(tb2,vtim,mpi_double_precision,0,icomm,ierrmpi)
199  call mpi_bcast(ta2,vtim,mpi_double_precision,0,icomm,ierrmpi)
200  endif
201 
202  ! spatial distribution
203 
204  vx=vtim*vlim
205  allocate(va1(vx))
206  allocate(va2(vx))
207  if (mype==0) then
208  write(filename,"(a)") "LeftV.dat"
209  inquire(file=filename,exist=alive)
210  if(alive) then
211  open(unit=22,file=filename,status='old')
212  do j=1,vx
213  read(22,*) va1(j)
214  enddo
215  close(22)
216  else
217  call randomv(va1,vtim,vlim,mnx,num,si,filename)
218  endif
219 
220  write(filename,"(a)") "RightV.dat"
221  inquire(file=filename,exist=alive)
222  if(alive) then
223  open(unit=22,file=filename,status='old')
224  do j=1,vx
225  read(22,*) va2(j)
226  enddo
227  close(22)
228  else
229  call randomv(va2,vtim,vlim,mnx,num,si,filename)
230  endif
231  endif
232  if (npe>1) then
233  call mpi_bcast(va1,vx,mpi_double_precision,0,icomm,ierrmpi)
234  call mpi_bcast(va2,vx,mpi_double_precision,0,icomm,ierrmpi)
235  endif
236 
237  end subroutine generatetv
238 
239 
240  subroutine randomt(tb,vtim,periods,variation,filename)
241  integer, intent(in) :: vtim
242  double precision, intent(in) :: periods,variation
243  double precision, intent(inout) :: tb(:)
244  character(len=100), intent(in) :: filename
245  integer :: i
246  double precision :: tt1,tt,mm
247 
248  ! generate random T, typically 300s +/- 75s
249 
250  tt = 0.d0
251  do i = 1, vtim
252  call random_number(mm)
253  tt1 = periods + variation * 2. * (mm - 0.5)
254  tt = tt + tt1
255  tb(i) = tt / unit_time
256  end do
257 
258  open(unit=21,file=filename,form='formatted',status='new')
259  write(21,'(es12.4)') tb
260  close(21)
261 
262  end subroutine randomt
263 
264  subroutine randomv(va,vtim,vlim,mnx,num,si,filename)
265  integer, intent(in) :: vtim,vlim,mnx,num
266  double precision, intent(in) :: si
267  double precision, intent(inout) :: va(:)
268  character(len=100), intent(in) :: filename
269  double precision, allocatable :: vn(:,:),vm(:)
270  double precision :: lambda,c,lambda_min,lambda_max,rms,ll
271  integer :: i
272  integer :: j, kk, lth = 1
273 
274  allocate(vn(num, mnx))
275  allocate(vm(mnx))
276  do kk = 1, vtim
277  lambda_min = lth * 1.d0 / (mnx * 1.d0)
278  lambda_max = lth * 2.d0
279  do i = 1, num
280  lambda = i * 1.d0 / num * (lambda_max - lambda_min) + lambda_min
281  c = lambda ** si
282  call random_number(ll)
283  do j = 1, mnx
284  vn(i,j) = c*sin(2.*dpi*(j*1.d0/(mnx-1)*lth)/lambda+2.*dpi*2.*(ll-0.5))
285  end do
286  end do
287  rms = 0.d0
288  vm = 0.d0
289  do j = 1, mnx
290  do i = 1, num
291  vm(j) = vm(j) + vn(i, j)
292  end do
293  rms = rms + vm(j)**2
294  end do
295  vm = vm / dsqrt(rms/mnx) ! normalization
296  do j = 1, vlim
297  va(vlim*(kk-1)+j) = vm(j)**2
298  end do
299  end do
300 
301  open(unit=22,file=filename,form='formatted',status='new')
302  write(22,'(es12.4)') va
303  close(22)
304 
305  end subroutine randomv
306 
307 end module mod_random_heating
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision unit_time
Physical scaling factor for time.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
integer ierrmpi
A global MPI error return code.
integer npe
The number of MPI tasks.
double precision, dimension(:), allocatable va2
double precision, dimension(:), allocatable tb1
double precision, dimension(:), allocatable ta2
double precision, dimension(:), allocatable ta1
double precision, dimension(:), allocatable va1
subroutine generatetv(vtim, vlim, periods, variation, mnx, num, si)
subroutine getlqgrid(lQgrid, ixIL, ixOL, qt, trelax, x, mode, vtim, vlim)
This module is for imposing randomized heating at x-axis for 2.5D simulation, for details please see ...
subroutine randomv(va, vtim, vlim, mnx, num, si, filename)
double precision, dimension(:), allocatable tb2
subroutine randomt(tb, vtim, periods, variation, filename)