MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_calculate_xw.t
Go to the documentation of this file.
1 !> Handles computations for coordinates and variables in output
3  implicit none
4 
5 contains
6 
7  !> Compute both corner as well as cell-centered values for output
8  subroutine calc_grid(qunit,igrid,xC,xCC,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP,normconv,&
9  ixC^L,ixCC^L,first)
10  ! this subroutine computes both corner as well as cell-centered values
11  ! it handles how we do the center to corner averaging, as well as
12  ! whether we switch to cartesian or want primitive or conservative output,
13  ! handling the addition of B0 in B0+B1 cases, ...
14  !
15  ! the normconv is passed on to usr_aux_output for extending with
16  ! possible normalization values for the nw+1:nw+nwauxio entries
17 
20  use mod_limiter
22  use mod_comm_lib, only: mpistop
23 
24  integer, intent(in) :: qunit, igrid
25  double precision, intent(in), dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
26  double precision, intent(in), dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
27  integer :: ixC^L,ixCC^L
28  logical, intent(in) :: first
29 
30  double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
31  double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
32  double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
33  double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
34  double precision,dimension(0:nw+nwauxio),intent(out) :: normconv
35 
36  double precision :: ldw(ixG^T), dwC(ixG^T)
37  double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC
38  double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC
39  double precision, dimension(ixG^T,1:nw+nwauxio) :: w
40  integer :: nxCC^D,idims,jxC^L,iwe
41  integer :: nx^D, nxC^D, ix^D, ix, iw, level, idir
42  logical, save :: subfirst=.true.
43 
44  ! initialize w
45  w=0.d0
46 
47  ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d; ! Corner indices
48  ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d; ! Center indices
49 
50  nx^d=ixmhi^d-ixmlo^d+1;
51  level=node(plevel_,igrid)
52 
53  normconv(0) = length_convert_factor
54  normconv(1:nw) = w_convert_factor
55  block=>ps(igrid)
56  w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
57 
58  if (nwextra>0) then
59  ! here we actually fill the ghost layers for the nwextra variables using
60  ! continuous extrapolation (as these values do not exist normally in ghost cells)
61  do idims=1,ndim
62  select case(idims)
63  {case(^d)
64  jxcmin^dd=ixghi^d+1-nghostcells^d%jxCmin^dd=ixglo^dd;
65  jxcmax^dd=ixghi^dd;
66  do ix^d=jxcmin^d,jxcmax^d
67  w(ix^d^d%jxC^s,nw-nwextra+1:nw) = w(jxcmin^d-1^d%jxC^s,nw-nwextra+1:nw)
68  end do
69  jxcmin^dd=ixglo^dd;
70  jxcmax^dd=ixglo^d-1+nghostcells^d%jxCmax^dd=ixghi^dd;
71  do ix^d=jxcmin^d,jxcmax^d
72  w(ix^d^d%jxC^s,nw-nwextra+1:nw) = w(jxcmax^d+1^d%jxC^s,nw-nwextra+1:nw)
73  end do \}
74  end select
75  end do
76  end if
77 
78  ! next lines needed when usr_aux_output uses gradients
79  ! and later on when dwlimiter2 is used
80  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
81  if(nwauxio>0)then
82  ! auxiliary io variables can be computed and added by user
83  ! next few lines ensure correct usage of routines like divvector etc
84  ! default (no) normalization for auxiliary variables
85  normconv(nw+1:nw+nwauxio)=one
86 
87  if (.not. associated(usr_aux_output)) then
88  call mpistop("usr_aux_output not defined")
89  else
90  call usr_aux_output(ixg^ll,ixm^ll^ladd1,w,ps(igrid)%x,normconv)
91  end if
92  endif
93 
94  ! In case primitives to be saved: use primitive subroutine
95  ! extra layer around mesh only needed when storing corner values and averaging
96  if(saveprim.and.first) call phys_to_primitive(ixg^ll,ixm^ll^ladd1,w(ixg^t,1:nw),ps(igrid)%x)
97 
98  if(b0field) then
99  ! B0+B1 split handled here
100  if(.not.saveprim.and.phys_energy) then
101  w(ixg^t,iw_e)=w(ixg^t,iw_e)+0.5d0*sum(ps(igrid)%B0(ixg^t,:,0)**2,dim=ndim+1) &
102  + sum(w(ixg^t,iw_mag(:))*ps(igrid)%B0(ixg^t,:,0),dim=ndim+1)
103  end if
104  w(ixg^t,iw_mag(:))=w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
105  end if
106  ! compute the cell-center values for w first
107  ! cell center values obtained from mere copy
108  wcc(ixcc^s,:)=w(ixcc^s,:)
109 
110  ! compute the corner values for w now by averaging
111 
112  if(slab_uniform) then
113  ! for slab symmetry: no geometrical info required
114  do iw=1,nw+nwauxio
115  {do ix^db=ixcmin^db,ixcmax^db\}
116  wc(ix^d,iw)=sum(w(ix^d:ix^d+1,iw))/dble(2**ndim)
117  {end do\}
118  end do
119  else
120  do iw=1,nw+nwauxio
121  {do ix^db=ixcmin^db,ixcmax^db\}
122  wc(ix^d,iw)=sum(w(ix^d:ix^d+1,iw)*ps(igrid)%dvolume(ix^d:ix^d+1)) &
123  /sum(ps(igrid)%dvolume(ix^d:ix^d+1))
124  {end do\}
125  end do
126  endif
127 
128  if(nocartesian) then
129  ! keep the coordinate and vector components
130  xc_tmp(ixc^s,1:ndim) = xc(ixc^s,1:ndim)
131  wc_tmp(ixc^s,1:nw+nwauxio) = wc(ixc^s,1:nw+nwauxio)
132  xcc_tmp(ixcc^s,1:ndim) = xcc(ixcc^s,1:ndim)
133  wcc_tmp(ixcc^s,1:nw+nwauxio) = wcc(ixcc^s,1:nw+nwauxio)
134  else
135  ! do all conversions to cartesian coordinates and vector components
136  ! start for the corner values
137  call to_cartesian(xc_tmp,wc_tmp,ixc^l,xc,wc)
138  ! then cell center values
139  call to_cartesian(xcc_tmp,wcc_tmp,ixcc^l,xcc,wcc)
140  endif
141 
142  ! Warning: differentiate between idl/idlCC/tecplot/tecplotCC/vtu(B)/vtu(B)CC
143  if(nwaux>0 .and. mype==0 .and. first.and.subfirst) then
144  ! when corner values are computed and auxiliaries present: warn!
145  if(convert_type=='idl'.or.convert_type=='tecplot' &
146  .or.convert_type=='vtu'.or.convert_type=='vtuB') &
147  write(*,*) 'Warning: also averaged auxiliaries within calc_grid'
148  subfirst=.false.
149  endif
150 
151  end subroutine calc_grid
152 
153  !> convert to cartesian coordinates and vector components
154  subroutine to_cartesian(x_TMP,w_TMP,ix^L,xC,wC)
155  ! conversion of coordinate and vector components from cylindrical/spherical
156  ! to cartesian coordinates and components done here
157  ! Also: nullifying values lower than smalldouble
159  use mod_geometry
160 
161  integer :: ix^L, ix^D, idim, iw, ivector, iw0
162  integer, dimension(nw) :: vectoriw
163  double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
164  double precision, dimension(ndim,ndim) :: normal
165 
166  double precision, dimension(ix^S,ndim) :: xC
167  double precision, dimension(ix^S,nw+nwauxio) :: wC
168 
169  double precision, dimension(ix^S,ndim) :: x_TMP
170  double precision, dimension(ix^S,nw+nwauxio) :: w_TMP
171 
172  iw0=0
173  vectoriw=-1
174  if(nvector>0) then
175  do ivector=1,nvector
176  do idim=1,ndim
177  vectoriw(iw_vector(ivector)+idim)=iw_vector(ivector)
178  end do
179  end do
180  endif
181  x_tec=0.d0
182  {do ix^db=ixmin^db,ixmax^db\}
183  select case (coordinate)
185  x_tec(1:ndim)=xc(ix^d,1:ndim)
186  w_tec(1:nw+nwauxio)=wc(ix^d,1:nw+nwauxio)
187  case (cylindrical)
188  {^ifoned
189  x_tec(1)=xc(ix^d,1)}
190  {^iftwod
191  select case (phi_)
192  case (2)
193  x_tec(1)=xc(ix^d,1)*cos(xc(ix^d,2))
194  x_tec(2)=xc(ix^d,1)*sin(xc(ix^d,2))
195  case default
196  x_tec(1)=xc(ix^d,1)
197  x_tec(2)=xc(ix^d,2)
198  end select}
199  {^ifthreed
200  x_tec(1)=xc(ix^d,1)*cos(xc(ix^d,phi_))
201  x_tec(2)=xc(ix^d,1)*sin(xc(ix^d,phi_))
202  x_tec(3)=xc(ix^d,z_)}
203 
204  if (nvector>0) then
205  {^ifoned normal(1,1)=one}
206 
207  {^iftwod
208  select case (phi_)
209  case (2)
210  normal(1,1)=cos(xc(ix^d,2))
211  normal(1,2)=-sin(xc(ix^d,2))
212  normal(2,1)=sin(xc(ix^d,2))
213  normal(2,2)=cos(xc(ix^d,2))
214  case default
215  normal(1,1)=one
216  normal(1,2)=zero
217  normal(2,1)=zero
218  normal(2,2)=one
219  end select}
220 
221  {^ifthreed
222  normal(1,1)=cos(xc(ix^d,phi_))
223  normal(1,phi_)=-sin(xc(ix^d,phi_))
224  normal(1,z_)=zero
225 
226  normal(2,1)=sin(xc(ix^d,phi_))
227  normal(2,phi_)=cos(xc(ix^d,phi_))
228  normal(2,z_)=zero
229 
230  normal(3,1)=zero
231  normal(3,phi_)=zero
232  normal(3,z_)=one}
233  end if
234  do iw=1,nw+nwauxio
235  if (iw<=nw) iw0=vectoriw(iw)
236  if (iw0>=0.and.iw<=iw0+ndim.and.iw<=nw) then
237  idim=iw-iw0
238  w_tec(iw0+idim)={^d&wc(ix^dd,iw0+^d)*normal(idim,^d)+}
239  else
240  w_tec(iw)=wc(ix^d,iw)
241  end if
242  end do
243  case (spherical)
244  x_tec(1)=xc(ix^d,1){^nooned*sin(xc(ix^d,2))}{^ifthreed*cos(xc(ix^d,3))}
245  {^iftwod
246  x_tec(2)=xc(ix^d,1)*cos(xc(ix^d,2))}
247  {^ifthreed
248  x_tec(2)=xc(ix^d,1)*sin(xc(ix^d,2))*sin(xc(ix^d,3))
249  x_tec(3)=xc(ix^d,1)*cos(xc(ix^d,2))}
250 
251  if (nvector>0) then
252  {^ifoned normal(1,1)=one}
253  {^nooned
254  normal(1,1)=sin(xc(ix^d,2)){^ifthreed*cos(xc(ix^d,3))}
255  normal(1,2)=cos(xc(ix^d,2)){^ifthreed*cos(xc(ix^d,3))
256  normal(1,3)=-sin(xc(ix^d,3))}}
257 
258  {^iftwod
259  normal(2,1)=cos(xc(ix^d,2))
260  normal(2,2)=-sin(xc(ix^d,2))}
261  {^ifthreed
262  normal(2,1)=sin(xc(ix^d,2))*sin(xc(ix^d,3))
263  normal(2,2)=cos(xc(ix^d,2))*sin(xc(ix^d,3))
264  normal(2,3)=cos(xc(ix^d,3))
265 
266  normal(3,1)=cos(xc(ix^d,2))
267  normal(3,2)=-sin(xc(ix^d,2))
268  normal(3,3)=zero}
269  end if
270  do iw=1,nw+nwauxio
271  if (iw<=nw) iw0=vectoriw(iw)
272  if (iw0>=0.and.iw<=iw0+ndim.and.iw<=nw) then
273  idim=iw-iw0
274  w_tec(iw0+idim)={^d&wc(ix^dd,iw0+^d)*normal(idim,^d)+}
275  else
276  w_tec(iw)=wc(ix^d,iw)
277  end if
278  end do
279  case default
280  write(*,*) "No converter for coordinate=",coordinate
281  end select
282  x_tmp(ix^d,1:ndim)=x_tec(1:ndim)
283  w_tmp(ix^d,1:nw+nwauxio)=w_tec(1:nw+nwauxio)
284  ! Be aware that small values are nullified here!!!
285  where(dabs(w_tmp(ix^d,1:nw+nwauxio))<smalldouble)
286  w_tmp(ix^d,1:nw+nwauxio)=zero
287  endwhere
288  {end do\}
289 
290  end subroutine to_cartesian
291 
292  !> get all variables names
293  subroutine getheadernames(wnamei,xandwnamei,outfilehead)
294  ! this collects all variables names in the wnamei character array, getting the info from
295  ! the prim_wnames/cons_wnames strings (depending on saveprim). It combines this info with names
296  ! for the dimensional directions in the xandwnamei array. In the outfilehead, it collects
297  ! the dimensional names, and only those names from the nw variables for output (through w_write)
298  ! together with the added names for nwauxio variables
301  use mod_geometry
302  use mod_comm_lib, only: mpistop
303 
304  character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
305  character(len=1024) :: outfilehead
306 
307  integer:: space_position,iw,ind
308  character(len=name_len):: wname
309  character(len=std_len):: aux_variable_names
310  character(len=std_len):: scanstring
311 
312  logical, save:: first=.true.
313 
314  ! in case additional variables are computed and stored for output
315  if (nwauxio>0) then
316  if (.not. associated(usr_add_aux_names)) then
317  call mpistop("usr_add_aux_names not defined")
318  else
319  call usr_add_aux_names(aux_variable_names)
320  end if
321  end if
322 
323  ! --- part to provide variable names from prim_wnames/cons_wnames strings
324  if(saveprim) then
325  scanstring=trim(aux_variable_names)
326  wnamei(1:nw)=prim_wnames(1:nw)
327  else
328  scanstring=trim(aux_variable_names)
329  wnamei(1:nw)=cons_wnames(1:nw)
330  endif
331 
332  space_position=index(scanstring,' ')
333  do iw=nw+1,nw+nwauxio
334  do while (space_position==1)
335  scanstring=scanstring(2:)
336  space_position=index(scanstring,' ')
337  enddo
338  wname=scanstring(:space_position-1)
339  scanstring=scanstring(space_position+1:)
340  space_position=index(scanstring,' ')
341 
342  ! fill all names, even those that we will not write away from the first nw
343  wnamei(iw)=trim(wname)
344  enddo
345  ! --- end of part to provide variable names
346 
347  select case (coordinate)
348  case( spherical )
349  xandwnamei(1)="r";{^nooned xandwnamei(2)="Theta"};{^ifthreed xandwnamei(3)="Phi"}
350  case( cylindrical )
351  xandwnamei(1)="R";
352  {^nooned
353  if( phi_ == 2 )then
354  xandwnamei(2)="Phi"
355  else
356  xandwnamei(2)="Z"
357  endif}
358  {^ifthreed
359  if( phi_ == 2 )then
360  xandwnamei(3)="Z"
361  else
362  xandwnamei(3)="Phi"
363  endif}
364  case default
365  xandwnamei(1)="X";{^nooned xandwnamei(2)="Y"};{^ifthreed xandwnamei(3)="Z"}
366  end select
367 
368  xandwnamei(ndim+1:ndim+nw+nwauxio)=wnamei(1:nw+nwauxio)
369 
370  ! in outfilehead, collect the dimensional names, and all output variable names
371  ! first all dimensions
372  write(outfilehead,'(a)') trim(xandwnamei(1))
373  {^nooned
374  do iw=2,ndim
375  wname=xandwnamei(iw)
376  write(outfilehead,'(a)')outfilehead(1:len_trim(outfilehead))//" "//trim(wname)
377  enddo
378  }
379  ! then all nw variables, with w_write control for inclusion
380  do iw=ndim+1,ndim+nw
381  wname=xandwnamei(iw)
382  if(w_write(iw-ndim)) then
383  write(outfilehead,'(a)')outfilehead(1:len_trim(outfilehead))//" "//trim(wname)
384  endif
385  enddo
386  ! then all nwauxio variables
387  if(nwauxio>0) then
388  do iw=ndim+nw+1,ndim+nw+nwauxio
389  wname=xandwnamei(iw)
390  write(outfilehead,'(a)')outfilehead(1:len_trim(outfilehead))//" "//trim(wname)
391  enddo
392  endif
393 
394  if(first.and.mype==0)then
395  print*,'-------------------------------------------------------------------------------'
396  write(unitterm,*)'Saving visual data. Coordinate directions and variable names are:'
397  ind=0
398  do iw=1,ndim
399  ind=ind+1
400  print *,ind,xandwnamei(iw)
401  enddo
402  do iw=1+ndim,nw+ndim
403  if(w_write(iw-ndim)) then
404  ind=ind+1
405  write(*,*) ind,wnamei(iw-ndim)
406  end if
407  end do
408  do iw=ndim+nw+1,ndim+nw+nwauxio
409  ind=ind+1
410  print *,ind,wnamei(iw-ndim)
411  enddo
412  write(unitterm,*)'time =', global_time
413  print*,'-------------------------------------------------------------------------------'
414  first=.false.
415  endif
416 
417  end subroutine getheadernames
418 
419  !> computes cell corner (xC) and cell center (xCC) coordinates
420  subroutine calc_x(igrid,xC,xCC)
422 
423  integer, intent(in) :: igrid
424  double precision, intent(out) :: xC(ixMlo^D-1:ixMhi^D,ndim)
425  double precision, intent(out) :: xCC(ixMlo^D:ixMhi^D,ndim)
426 
427  integer :: ixC^L, idims, level, ix
428 
429  level=node(plevel_,igrid)
430 
431  ! coordinates of cell centers
432  xcc(ixm^t,:)=ps(igrid)%x(ixm^t,:)
433 
434  ! coordinates of cell corners
435  ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
436  if(slab_uniform)then
437  do idims=1,ndim
438  xc(ixc^s,idims)=ps(igrid)%x(ixc^s,idims)+0.5d0*dx(idims,level)
439  end do
440  else
441  ! for any non-cartesian or stretched coordinate (allow multiple stretched directions)
442  {do ix=ixcmin^d,ixcmax^d
443  xc(ix^d%ixC^s,^d)=ps(igrid)%x(ix^d%ixC^s,^d)+0.5d0*ps(igrid)%dx(ix^d%ixC^s,^d)
444  end do\}
445  endif
446 
447  end subroutine calc_x
448 
449 end module mod_calculate_xw
Handles computations for coordinates and variables in output.
subroutine calc_grid(qunit, igrid, xC, xCC, xC_TMP, xCC_TMP, wC_TMP, wCC_TMP, normconv, ixCL, ixCCL, first)
Compute both corner as well as cell-centered values for output.
subroutine getheadernames(wnamei, xandwnamei, outfilehead)
get all variables names
subroutine calc_x(igrid, xC, xCC)
computes cell corner (xC) and cell center (xCC) coordinates
subroutine to_cartesian(x_TMP, w_TMP, ixL, xC, wC)
convert to cartesian coordinates and vector components
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
integer coordinate
Definition: mod_geometry.t:7
integer, parameter spherical
Definition: mod_geometry.t:11
integer, parameter cartesian
Definition: mod_geometry.t:8
integer, parameter cylindrical
Definition: mod_geometry.t:10
integer, parameter cartesian_expansion
Definition: mod_geometry.t:12
integer, parameter cartesian_stretched
Definition: mod_geometry.t:9
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision, dimension(:), allocatable w_convert_factor
Conversion factors the primitive variables.
type(state), pointer block
Block pointer for using one block and its previous state.
integer ixghi
Upper index of grid block arrays.
double precision global_time
The global simulation time.
logical saveprim
If true, convert from conservative to primitive variables in output.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer mype
The rank of the current MPI task.
integer, parameter plevel_
double precision length_convert_factor
integer ixm
the mesh range of a physical block without ghost cells
integer nwauxio
Number of auxiliary variables that are only included in output.
integer, parameter unitterm
Unit for standard output.
logical, dimension(:), allocatable w_write
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer, dimension(:,:), allocatable node
double precision, dimension(ndim) dxlevel
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Module with slope/flux limiters.
Definition: mod_limiter.t:2
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:56
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition: mod_physics.t:16
logical phys_energy
Solve energy equation or not.
Definition: mod_physics.t:27
Module with all the methods that users can customize in AMRVAC.
procedure(aux_output), pointer usr_aux_output
procedure(add_aux_names), pointer usr_add_aux_names