MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_calculate_xw.t
Go to the documentation of this file.
1!> Handles computations for coordinates and variables in output
3 implicit none
4
5contains
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 .and. allocated(iw_mag)) 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 double precision, dimension(ix^S,ndim) :: xC
162 double precision, dimension(ix^S,nw+nwauxio) :: wC
163 double precision, dimension(ix^S,ndim) :: x_TMP
164 double precision, dimension(ix^S,nw+nwauxio) :: w_TMP
165 integer :: ix^L
166
167 double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
168 double precision, dimension(ndim,ndim) :: normal
169 integer, dimension(nw) :: vectoriw
170 integer :: ix^D, idim, iw, ivector, iw0
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 logical, save:: first=.true.
309 character(len=name_len):: wname
310 character(len=std_len):: aux_variable_names
311 character(len=std_len):: scanstring
312
313 ! in case additional variables are computed and stored for output
314 if (nwauxio>0) then
315 if (.not. associated(usr_add_aux_names)) then
316 call mpistop("usr_add_aux_names not defined")
317 else
318 call usr_add_aux_names(aux_variable_names)
319 end if
320 end if
321
322 ! --- part to provide variable names from prim_wnames/cons_wnames strings
323 if(saveprim) then
324 scanstring=trim(aux_variable_names)
325 wnamei(1:nw)=prim_wnames(1:nw)
326 else
327 scanstring=trim(aux_variable_names)
328 wnamei(1:nw)=cons_wnames(1:nw)
329 endif
330
331 space_position=index(scanstring,' ')
332 do iw=nw+1,nw+nwauxio
333 do while (space_position==1)
334 scanstring=scanstring(2:)
335 space_position=index(scanstring,' ')
336 enddo
337 wname=scanstring(:space_position-1)
338 scanstring=scanstring(space_position+1:)
339 space_position=index(scanstring,' ')
340
341 ! fill all names, even those that we will not write away from the first nw
342 wnamei(iw)=trim(wname)
343 enddo
344 ! --- end of part to provide variable names
345
346 select case (coordinate)
347 case( spherical )
348 xandwnamei(1)="r";{^nooned xandwnamei(2)="Theta"};{^ifthreed xandwnamei(3)="Phi"}
349 case( cylindrical )
350 xandwnamei(1)="R";
351 {^nooned
352 if( phi_ == 2 )then
353 xandwnamei(2)="Phi"
354 else
355 xandwnamei(2)="Z"
356 endif}
357 {^ifthreed
358 if( phi_ == 2 )then
359 xandwnamei(3)="Z"
360 else
361 xandwnamei(3)="Phi"
362 endif}
363 case default
364 xandwnamei(1)="X";{^nooned xandwnamei(2)="Y"};{^ifthreed xandwnamei(3)="Z"}
365 end select
366
367 xandwnamei(ndim+1:ndim+nw+nwauxio)=wnamei(1:nw+nwauxio)
368
369 ! in outfilehead, collect the dimensional names, and all output variable names
370 ! first all dimensions
371 write(outfilehead,'(a)') trim(xandwnamei(1))
372 {^nooned
373 do iw=2,ndim
374 wname=xandwnamei(iw)
375 write(outfilehead,'(a)')outfilehead(1:len_trim(outfilehead))//" "//trim(wname)
376 enddo
377 }
378 ! then all nw variables, with w_write control for inclusion
379 do iw=ndim+1,ndim+nw
380 wname=xandwnamei(iw)
381 if(w_write(iw-ndim)) then
382 write(outfilehead,'(a)')outfilehead(1:len_trim(outfilehead))//" "//trim(wname)
383 endif
384 enddo
385 ! then all nwauxio variables
386 if(nwauxio>0) then
387 do iw=ndim+nw+1,ndim+nw+nwauxio
388 wname=xandwnamei(iw)
389 write(outfilehead,'(a)')outfilehead(1:len_trim(outfilehead))//" "//trim(wname)
390 enddo
391 endif
392
393 if(first.and.mype==0)then
394 print*,'-------------------------------------------------------------------------------'
395 write(unitterm,*)'Saving visual data. Coordinate directions and variable names are:'
396 ind=0
397 do iw=1,ndim
398 ind=ind+1
399 print *,ind,xandwnamei(iw)
400 enddo
401 do iw=1+ndim,nw+ndim
402 if(w_write(iw-ndim)) then
403 ind=ind+1
404 write(*,*) ind,wnamei(iw-ndim)
405 end if
406 end do
407 do iw=ndim+nw+1,ndim+nw+nwauxio
408 ind=ind+1
409 print *,ind,wnamei(iw-ndim)
410 enddo
411 write(unitterm,*)'time =', global_time
412 print*,'-------------------------------------------------------------------------------'
413 first=.false.
414 endif
415
416 end subroutine getheadernames
417
418 !> computes cell corner (xC) and cell center (xCC) coordinates
419 subroutine calc_x(igrid,xC,xCC)
421
422 integer, intent(in) :: igrid
423 double precision, intent(out) :: xC(ixMlo^D-1:ixMhi^D,ndim)
424 double precision, intent(out) :: xCC(ixMlo^D:ixMhi^D,ndim)
425
426 integer :: ixC^L, idims, level, ix
427
428 level=node(plevel_,igrid)
429
430 ! coordinates of cell centers
431 xcc(ixm^t,:)=ps(igrid)%x(ixm^t,:)
432
433 ! coordinates of cell corners
434 ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
435 if(slab_uniform)then
436 do idims=1,ndim
437 xc(ixc^s,idims)=ps(igrid)%x(ixc^s,idims)+0.5d0*dx(idims,level)
438 end do
439 else
440 ! for any non-cartesian or stretched coordinate (allow multiple stretched directions)
441 {do ix=ixcmin^d,ixcmax^d
442 xc(ix^d%ixC^s,^d)=ps(igrid)%x(ix^d%ixC^s,^d)+0.5d0*ps(igrid)%dx(ix^d%ixC^s,^d)
443 end do\}
444 endif
445
446 end subroutine calc_x
447
448end module mod_calculate_xw
Handles computations for coordinates and variables in output.
subroutine to_cartesian(x_tmp, w_tmp, ixl, xc, wc)
convert to cartesian coordinates and vector components
subroutine getheadernames(wnamei, xandwnamei, outfilehead)
get all variables names
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 calc_x(igrid, xc, xcc)
computes cell corner (xC) and cell center (xCC) coordinates
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module with geometry-related routines (e.g., divergence, curl)
Definition mod_geometry.t:2
integer coordinate
Definition mod_geometry.t:7
integer, parameter spherical
integer, parameter cartesian
Definition mod_geometry.t:8
integer, parameter cylindrical
integer, parameter cartesian_expansion
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.
double precision length_convert_factor
Conversion factor for length unit.
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
if true write the w variable in output
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.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer, dimension(:,:), allocatable node
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:55
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition mod_physics.t:50
logical phys_energy
Solve energy equation or not.
Definition mod_physics.t:35
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