DISPATCH
mesh_mod.f90
1 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
2 !> Template module for mesh
3 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
4 MODULE mesh_mod
5  USE trace_mod
6  USE io_mod
7  USE mpi_mod
8  USE kinds_mod
9  implicit none
10  private
11  !-----------------------------------------------------------------------------
12  ! A generic orthogonal, curvilinear mesh type.
13  !
14  ! Rather than have one type for each cardinal direction, all metric factors
15  ! for all directions are defined here, but only those for the current axis
16  ! are allocated upon initialisation. This way, you can have one type for meshes
17  ! in any direction.
18  !
19  ! Primary scale factors: h2c, h2f, h31c, h31f, h32c, h32f
20  ! `c` means the quantity is located at cell-centre; `f` means the quantity
21  ! is located at the cell-face. Regardless of whether you are using a
22  ! staggered method or not, you will generally need these factors if using
23  ! curvilinear coords.
24  ! Secondary volume and area factors: dvol1c, dvol1f, dvol2c, dvol2f, dvol3c,
25  ! dvol3f, dar1c, dar1f, dar2c, dar2f,
26  ! dar31c, dar31f, dar32c, dar32f
27  !
28  ! Current requirements for curvilinear meshes:
29  ! - The scale factors are either a function of one variable or are separable
30  ! (e.g., h3 = r sin \theta in spherical coords. but can be expressed as
31  ! h31(i) = r(i) and h32(j) = sin theta(j)).
32  ! - The mesh spacing (self%d) is assumed to be uniform in all directions.
33  !-----------------------------------------------------------------------------
34  !-----------------------------------------------------------------------------
35  type, abstract, public :: mesh_t
36  integer :: id=0 ! copy of task id, for info
37  integer :: n = 16 ! number of "active" zones
38  integer :: ng = 2 ! number of ghost zones
39  integer :: nc ! number of cells
40  integer :: lb ! lower boundary of mesh incl. ghost zones
41  integer :: li ! lower boundary of mesh excl. ghost zones
42  integer :: lo
43  integer :: ub ! upper boundary of mesh incl. ghost zones
44  integer :: ui ! upper boundary of mesh excl. ghost zones
45  integer :: uo
46  integer :: gn = 20 ! total number of zones incl. ghost zones (= n + 2 * ng)
47  integer :: idir ! identify which direction the mesh applies to
48  real(8) :: p = 0d0 ! position of patch centre in **Cartesian** coords.
49  real(8) :: s = 0d0 ! spatial extent
50  real(8) :: d = 0d0 ! grid spacing
51  real(8) :: b = 1d0 ! box size
52  real(8) :: o = 0.0 ! position offset in index space
53  real(4) :: lf = 0.0, uf = 0.0 ! interpolation boundary indices (must be floating point)
54  ! the variables `llc_cart` and `llc_nat` describe the same location, but in
55  ! different reference frames.
56  real(8) :: llc_cart = 0.0 ! position of the patch's lower left corner in Cartesian
57  ! coords.; replacement for `p`.
58  real(8) :: llc_nat = 0.0 ! position of the patch's lower left corner in relative,
59  ! native coords.; always = 0 in Cartesian coords.
60  real(8) :: centre_nat ! volume centroid of the mesh in its *own* coords.
61  real(8), dimension(:), pointer :: r => null() ! coordinate relative to patch centre
62  real(8), dimension(:), pointer :: h => null() ! offset for staggered variables;
63  ! should be initialised in `mhd_t%init`
64  real, dimension(:), pointer :: ddn => null() ! borrow for now (FIXME)
65  real, dimension(:), pointer :: dup => null() ! borrow for now (FIXME)
66  real, dimension(:), pointer :: dsup => null() ! borrow for now (FIXME)
67  real, dimension(:), pointer :: dsdn => null() ! borrow for now (FIXME)
68  real, dimension(:), pointer :: b_m => null() ! borrow for now (FIXME) -- bifrost xm/ym/zm
69  real, dimension(:), pointer :: b_mdn => null() ! borrow for now (FIXME) -- bifrost xmdn/ymdn/zmdn
70  real, dimension(:), pointer :: b_1d => null() ! borrow for now (FIXME) -- bifrost dx1d/dy1d/dz1d
71  real, dimension(:), pointer :: b_didup => null() ! borrow for now (FIXME) -- bifrost dxidxup/...
72  real, dimension(:), pointer :: b_diddn => null() ! borrow for now (FIXME) -- bifrost dxidxdn/...
73  real, dimension(:,:), pointer:: zlincoeff => null()! borrow for now (FIXME) -- bifrost zlincoeff
74  real, dimension(:,:), pointer:: dzlincoeff => null()! borrow for now (FIXME) -- bifrost zlincoeff
75  real, dimension(:,:), pointer:: zuincoeff => null()! borrow for now (FIXME) -- bifrost zlincoeff
76  real, dimension(:,:), pointer:: dzuincoeff => null()! borrow for now (FIXME) -- bifrost zlincoeff
77  real(8) :: b_d ! borrow for now (FIXME) -- grid spacing in bifrost
78  real :: b_ad, b_bd, b_cd !borrow for now (FIXME) -- bifrost cdx/bdx/adx...
79  integer(4):: b_s ! borrow for now (FIXME) -- same as li, but in bifrost
80  integer(4):: b_sb ! borrow for now (FIXME) -- same as lb, but in bifrost
81  integer(4):: b_e ! borrow for now (FIXME) -- same as ui, but in bifrost
82  integer(4):: b_eb ! borrow for now (FIXME) -- sane as ub, but in bifrost
83  integer(kind=4) :: b_r, b_rb
84 
85  character(len=16):: type = ''
86  real(8), allocatable :: zeta ! position angle of the phi = 0 axis relative to the x-axis (cyl. and sph. coords.).
87  real(8), allocatable :: xi ! position angle of the theta = 0 relative to the neg. z-axis (sph. coords.).
88  real(8):: origin=0d0
89  !---------------------------------------------------------------------------
90  ! Physical boundary condition info; indices beyond which BCs rule
91  !---------------------------------------------------------------------------
92  logical:: lower_boundary = .false.
93  logical:: upper_boundary = .false.
94  logical:: no_mans_land = .false.
95 
96  procedure(), pointer, nopass :: currenttocartesian => null()
97  ! geometric factors as procedure pointers
98  procedure(func_interface), pointer :: h2c => null()
99  procedure(func_interface), pointer :: h2f => null()
100  procedure(func_interface), pointer :: h31c => null()
101  procedure(func_interface), pointer :: h31f => null()
102  procedure(func_interface), pointer :: h32f => null()
103  procedure(func_interface), pointer :: h32c => null()
104  procedure(func_interface), pointer :: dvol1c => null()
105  procedure(func_interface), pointer :: dvol1f => null()
106  procedure(func_interface), pointer :: dvol2c => null()
107  procedure(func_interface), pointer :: dvol2f => null()
108  procedure(func_interface), pointer :: dvol3c => null()
109  procedure(func_interface), pointer :: dvol3f => null()
110  procedure(func_interface), pointer :: dar1c => null()
111  procedure(func_interface), pointer :: dar1f => null()
112  procedure(func_interface), pointer :: dar2c => null()
113  procedure(func_interface), pointer :: dar2f => null()
114  procedure(func_interface), pointer :: dar31c => null()
115  procedure(func_interface), pointer :: dar31f => null()
116  procedure(func_interface), pointer :: dar32c => null()
117  procedure(func_interface), pointer :: dar32f => null()
118  contains
119  procedure :: init => init_mesh
120  procedure :: deallocate => deallocate_coord
121  procedure :: print => print_mesh
122  procedure(init_geometryinterface), deferred :: init_geometry
123  end type
124  abstract interface
125  subroutine init_geometryinterface(self, idir, origin)
126  import mesh_t
127  class(mesh_t) :: self
128  integer, intent(in) :: idir
129  real(8), intent(in) :: origin
130  end subroutine
131 
132  function func_interface(self, i) result(x)
133  import mesh_t
134  class(mesh_t) :: self
135  integer, intent(in) :: i
136  real(8) :: x
137  end function
138  end interface
139 
140  !-----------------------------------------------------------------------------
141  ! Utility list of coordinate systems available for meshes
142  !-----------------------------------------------------------------------------
143  type :: mesh_type
144  integer :: cartesian = 1
145  integer :: spherical = 2
146  integer :: cylindrical = 3
147  end type
148  type(mesh_type), public :: mesh_types
149 
150  !-----------------------------------------------------------------------------
151  ! Cartesian: x,y,z
152  !-----------------------------------------------------------------------------
153  type, extends(mesh_t), public :: cartesian_mesh
154  contains
155  procedure :: init_geometry => init_cartesian_geometry
156  end type cartesian_mesh
157 
158  !-----------------------------------------------------------------------------
159  ! cylindrical: z, R, Phi
160  ! Assert that, when Phi is 0, the r-coord. and x-coord. overlie one another.
161  ! The z-coords. should always overlap.
162  ! Further, the angle `zeta` specifies the offset between the Phi = 0 and x axes.
163  !-----------------------------------------------------------------------------
164  type, extends(mesh_t), public :: cylindrical_mesh
165  contains
166  procedure :: init_geometry => init_cylindrical_geometry
167  end type cylindrical_mesh
168 
169  !-----------------------------------------------------------------------------
170  ! spherical: r, theta, phi
171  ! Assert that, when theta = 0, the r-coord. and the *negative* z-coord overlie
172  ! one another.
173  ! Further, assert that, when phi = 0 and theta = pi/2, the r-coord. and x-coord.
174  ! overlie one another.
175  ! The angle `zeta` specifies the offset between the phi = 0 and x axes, while
176  ! the angle `xi` specified the offset between the theta=0 and neg. z axes.
177  !-----------------------------------------------------------------------------
178  type, extends(mesh_t), public :: spherical_mesh
179  contains
180  procedure :: init_geometry => init_spherical_geometry
181  end type spherical_mesh
182 
184  module procedure :: cylindricaltocartesian4, cylindricaltocartesian8
185  end interface
186  interface cylindricaltospherical
187  module procedure :: cylindricaltospherical4, cylindricaltospherical8
188  end interface
189  interface sphericaltocartesian
190  module procedure :: sphericaltocartesian4, sphericaltocartesian8
191  end interface
192  interface cartesiantospherical
193  module procedure :: cartesiantospherical4, cartesiantospherical8
194  end interface
195  interface cartesiantocylindrical
196  module procedure :: cartesiantocylindrical4, cartesiantocylindrical8
197  end interface
198  interface sphericaltocylindrical
199  module procedure :: sphericaltocylindrical4, sphericaltocylindrical8
200  end interface
201  interface currenttocartesian
202  module procedure :: currenttocartesian4, currenttocartesian8
203  end interface
204  interface cartesiantocurrent
205  module procedure :: cartesiantocurrent4, cartesiantocurrent8
206  end interface
207  interface round
208  module procedure :: round4, round8
209  end interface
210 
211  real(8), parameter :: pi = 3.14159265358979323846
212 
213 PUBLIC meshfactory, meshrecycler, &
214  cartesiantocylindrical, cartesiantospherical, cylindricaltocartesian, &
215  cylindricaltospherical, sphericaltocartesian, sphericaltocylindrical, &
216  currenttocartesian, cartesiantocurrent, round
217 CONTAINS
218 
219 !===============================================================================
220 !> Create a mesh collection.
221 !> This method creates three orthogonal meshes to form an 3-dimensional grid.
222 !> It is possible for up to 2 of the meshes to have extent one.
223 !===============================================================================
224 SUBROUTINE meshfactory (meshtype, collection, n, ng, s, p, llc_native, nml)
225  class(mesh_t), dimension(:), pointer :: collection
226  integer, intent(in) :: meshtype
227  integer, dimension(3), optional :: ng
228  integer, dimension(3) :: ngi
229  integer, dimension(3), optional :: n
230  integer, dimension(3) :: ni
231  real(8), dimension(3), optional :: s, p
232  real(8), dimension(3) :: si, pi
233  logical, optional :: nml
234  integer :: i
235  ! expected form of `origin` for Cartesian coords.: [0,0,0] (unused)
236  ! for cylindrical coords.: [rmin, 0, zeta]
237  ! for spherical coords.: [rmin, xi, zeta]
238  real(8), dimension(3), optional :: llc_native
239  real(8), dimension(3) :: llc_native_i
240  !.............................................................................
241  ! set reasonable defaults
242  if (present(n)) then
243  ni = n
244  else
245  ni(:) = 16
246  end if
247  if (present(s)) then
248  si = s
249  else
250  si(:) = 1.0
251  end if
252  if (present(p)) then
253  pi = p
254  else
255  pi(:) = 0.0
256  end if
257  if (present(llc_native)) then
258  llc_native_i = llc_native
259  else
260  llc_native_i(:) = 0.0
261  end if
262  if (present(ng)) then
263  ngi = ng
264  else
265  ngi(:) = 2
266  end if
267 
268 
269  if (associated(collection)) call meshrecycler(collection)
270  ! Allocate based on the desired mesh type.
271  ! The grid is always assumed to be 3-D, even if one of the directions is
272  ! collapsed/a symmetry direction (in which case, n = 1).
273  if (meshtype .eq. mesh_types%Cartesian) then
274  allocate(cartesian_mesh :: collection(3))
275  collection%type = 'Cartesian_mesh'
276  else if ( meshtype .eq. mesh_types%spherical ) then
277  allocate(spherical_mesh :: collection(3))
278  collection%type = 'spherical_mesh'
279  else if ( meshtype .eq. mesh_types%cylindrical ) then
280  allocate(cylindrical_mesh :: collection(3))
281  collection%type = 'cylindrical_mesh'
282  else
283  call mpi%abort("Invalid mesh type. Abort!")
284  end if
285  if (present(nml)) then
286  collection%no_mans_land = nml
287  end if
288 
289  do i=1,3
290  ! initialise the mesh
291  collection(i)%idir = i
292  call collection(i)%init(ni(i),ngi(i),si(i),pi(i),llc_native_i(i))
293 
294  ! initialise geometric scale factors (if applicable).
295  call collection(i)%init_geometry(i,llc_native_i(i))
296  end do
297 
298 END SUBROUTINE meshfactory
299 
300 !===============================================================================
301 !> A means to entirely deallocate a collection of meshes
302 !===============================================================================
303 SUBROUTINE meshrecycler (collection)
304  class(mesh_t), dimension(:), pointer :: collection
305  integer :: i
306  !.............................................................................
307  do i=1,size(collection)
308  call collection(i)%deallocate
309  end do
310  deallocate(collection)
311  nullify (collection)
312 END SUBROUTINE meshrecycler
313 
314 !===============================================================================
315 !> The integer mesh is defined by specifying the number of points `n` and the
316 !> number of guard cells `ng`. The real mesh is then fully specified by giving
317 !> the size `s`. The center of the mesh, `p` can also be supplied.
318 !===============================================================================
319 SUBROUTINE init_mesh (self, n, ng, s, p, llc_native)
320  class(mesh_t) :: self
321  integer :: i
322  integer, optional :: n, ng
323  real(8), optional :: s, p
324  real(8), optional :: llc_native
325  real(8):: h
326  !.............................................................................
327  call trace_begin ('mesh%init')
328 
329  if (present(n)) then
330  self%n = n
331  endif
332  if (present(ng)) then
333  self%ng = ng
334  endif
335  if (present(s)) then
336  self%s = s
337  endif
338  if (present(p)) then
339  self%p = p
340  endif
341  if (present(llc_native)) then
342  self%llc_nat = llc_native ! the lower left corner in the current patch coords.
343  endif
344 
345  ! initialise mesh indices
346  self%lb = 1
347  self%li = self%lb + self%ng
348  self%lo = self%li - 1
349  self%ui = self%lo + self%n
350  self%uo = self%ui + 1
351  self%ub = self%ui + self%ng
352  ! account for collapsed dimensions
353  if (self%n == 1) then
354  self%lo = self%li
355  self%uo = self%ui
356  self%ub = self%ui
357  end if
358  self%gn = self%ub
359  allocate (self%r(self%gn))
360  ! initialise mesh floating point attributes
361  self%o = (self%li+self%ui)*0.5
362 
363  ! calculate the volume centroid of the current mesh
364  self%centre_nat = volumecentroid(self)
365  ! cell-centred coordinate relative to the patch *centre*; the absolute position
366  ! of a cell can be recovered using `self%centre_nat + self%r(i)` in the coord. system
367  ! of the patch, and `self%p + self%centre_nat + self%r(i)`.
368  ! in Cartesian coords., self%centre_nat = self%p
369  if (self%no_mans_land) then
370  h = 0.5_8
371  if (self%n == 1) h = 1.0_8 ! account for collapsed dimensions (otherwise `r` ends up non-zero).
372  self%nc = self%ui - self%li + 1
373  else
374  h = 1.0_8
375  self%nc = max(1, self%ui - self%li)
376  end if
377  self%d = self%s/max(1,self%nc)
378  do i=1,self%gn
379  self%r(i) = real(i - h - self%ng,kind=8) * self%d + self%llc_nat - self%centre_nat
380  end do
381  self%lf = -self%s*0.5_8
382  self%uf = +self%s*0.5_8
383  call trace_end
384 END SUBROUTINE init_mesh
385 
386 !===============================================================================
387 !> deallocates the position coordinate of a mesh
388 !===============================================================================
389 SUBROUTINE deallocate_coord (self)
390  class(mesh_t) :: self
391  !.............................................................................
392  if (associated(self%r)) deallocate (self%r)
393  if (associated(self%h)) deallocate (self%h)
394  if (associated(self%ddn)) deallocate (self%ddn)
395  if (associated(self%dup)) deallocate (self%dup)
396 END SUBROUTINE deallocate_coord
397 
398 !===============================================================================
399 !> print mesh information to the stdout
400 !===============================================================================
401 SUBROUTINE print_mesh (self)
402  class(mesh_t):: self
403  integer:: i
404  !.............................................................................
405 
406  if (io%verbose <= 0) return
407  print*,'::mesh%print::'
408  select type(self)
409  type is(cartesian_mesh)
410  print*,'Cartesian', self%idir
411  type is(spherical_mesh)
412  print*,'spherical', self%idir
413  type is(cylindrical_mesh)
414  print*,'cylindrical', self%idir
415  end select
416  print'(a,3i8)', "n, ng, gn: ", self%n, self%ng, self%gn
417  print'(a,6i8)', "lb, lo, li, ub, uo, ui: ", self%lb, self%lo, self%li, self%ub, self%uo, self%ui
418  print'(a,3x,3g16.8)', "lf, uf, p: ",self%lf, self%uf, self%p
419  print'(a,3x,3g16.8)', "s, d, b: ", self%s, self%d, self%b
420  print'(a,3x,4g16.8)', "llc_nat, centre_nat, llc_cart, o: ",self%llc_nat, self%centre_nat, self%llc_cart, self%o
421  do i=1,self%gn
422  print'(i4,f14.8)',i,self%r(i)
423  end do
424 
425  if (allocated(self%zeta)) print*,"phi axis offset: ", self%zeta
426  if (allocated(self%xi)) print*,"theta axis offset: ", self%xi
427 
428  print*,'----------'
429 END SUBROUTINE print_mesh
430 
431 !===============================================================================
432 !> Dummy initialise for a Cartesian mesh: All scale factors are one and all
433 !> volume elements are equal to dx.
434 !===============================================================================
435 SUBROUTINE init_cartesian_geometry (self, idir, origin)
436  class(cartesian_mesh) :: self
437  integer, intent(in) :: idir
438  real(8), intent(in) :: origin
439  integer :: n
440  !.............................................................................
441  ! All the metric factors are 1 in Cartesian coords. and all the volume coords.
442  ! are simply the cell sizes.
443 
444  select case (idir)
445  case (1)
446  self%h2c => func_one
447  self%h2f => func_one
448  self%h31c => func_one
449  self%h31f => func_one
450  self%dvol1c => func_dx
451  self%dvol1f => func_dx
452 
453  self%dar1c => func_one
454  self%dar1f => func_one
455  self%dar2c => func_one
456  self%dar2f => func_one
457  self%dar31c => func_one
458  self%dar31f => func_one
459  case (2)
460  self%h32c => func_one
461  self%h32f => func_one
462  self%dvol2c => func_dx
463  self%dvol2f => func_dx
464 
465  self%dar32c => func_one
466  self%dar32f => func_one
467  case (3)
468  self%dvol3c => func_dx
469  self%dvol3f => func_dx
470  end select
471 
472 END SUBROUTINE init_cartesian_geometry
473 
474 !===============================================================================
475 !> Initalise the geometric factors for a cylindrical mesh.
476 !> A cylindrical mesh is defined as having the form (z,R,\phi).
477 !===============================================================================
478 SUBROUTINE init_cylindrical_geometry (self, idir, origin)
479  class(cylindrical_mesh) :: self
480  integer, intent(in) :: idir
481  real(8), intent(in) :: origin
482  !.............................................................................
483  select case (idir)
484  case (1)
485  self%h2c => func_one
486  self%h2f => func_one
487  self%h31c => func_one
488  self%h31f => func_one
489  self%dvol1c => func_dx
490  self%dvol1f => func_dx
491 
492  self%dar1c => func_one
493  self%dar1f => func_one
494  self%dar2c => func_one
495  self%dar2f => func_one
496  self%dar31c => func_one
497  self%dar31f => func_one
498  case (2)
499  self%h32c => func_radius_c
500  self%h32f => func_radius_f
501  self%dvol2c => func_dr2_c
502  self%dvol2f => func_dr2_f
503 
504  self%dar32c => func_dar32_c
505  self%dar32f => func_dar32_f
506  case (3)
507  self%dvol3c => func_dx
508  self%dvol3f => func_dx
509 
510  ! when phi = 0 and the radial vector overlies the x-axis, then zeta == 0.
511  allocate(self%zeta)
512  self%zeta = origin
513  end select
514 
515 END SUBROUTINE init_cylindrical_geometry
516 
517 !===============================================================================
518 !> Initalise the geometry factors for a spherical mesh.
519 !> A spherical mesh is defined as having the form (r,\theta,\phi).
520 !===============================================================================
521 SUBROUTINE init_spherical_geometry (self, idir, origin)
522  class(spherical_mesh) :: self
523  integer, intent(in) :: idir
524  real(8), intent(in) :: origin
525  !.............................................................................
526 
527  select case (idir)
528  case (1)
529  self%h2c => func_radius_c
530  self%h2f => func_radius_f
531  self%h31c => func_radius_c
532  self%h31f => func_radius_f
533  self%dvol1c => func_dr3_c
534  self%dvol1f => func_dr3_f
535 
536  self%dar1c => func_dar1_c
537  self%dar1f => func_dar1_f
538  self%dar2c => func_dar2_c
539  self%dar2f => func_dar2_f
540  self%dar31c => func_dar31_c
541  self%dar31f => func_dar31_f
542  case (2)
543  self%h32c => func_sintheta_c
544  self%h32f => func_sintheta_f
545  self%dvol2c => func_dcostheta_c
546  self%dvol2f => func_dcostheta_f
547 
548  self%dar32c => func_dar32_c
549  self%dar32f => func_dar32_f
550 
551  ! when phi = theta = 0 and the radial vector overlies the x-axis, then xi == 0.
552  allocate(self%xi)
553  self%xi = origin
554  case (3)
555  self%dvol3c => func_dx
556  self%dvol3f => func_dx
557 
558  ! when phi = theta = 0 and the radial vector overlies the x-axis, then zeta == 0.
559  allocate(self%zeta)
560  self%zeta = origin
561  end select
562 
563 END SUBROUTINE init_spherical_geometry
564 
565 !===============================================================================
566 !> Collection of functions used to calculate geometric factors.
567 !===============================================================================
568 FUNCTION func_one (self, i) RESULT(one)
569  class(mesh_t) :: self
570  integer, intent(in) :: i
571  real(8) :: one
572  !.............................................................................
573  one = 1.0
574 
575 END FUNCTION func_one
576 
577 FUNCTION func_dx (self, i) RESULT(dx)
578  class(mesh_t) :: self
579  integer, intent(in) :: i
580  real(8) :: dx
581  !.............................................................................
582  dx = self%d
583 
584 END FUNCTION func_dx
585 
586 FUNCTION func_position_c (self, i) RESULT(r)
587  class(mesh_t) :: self
588  integer, intent(in) :: i
589  real(8) :: r
590  !.............................................................................
591  r = self%centre_nat + self%r(i)
592 
593 END FUNCTION func_position_c
594 
595 FUNCTION func_position_f (self, i) RESULT(r)
596  class(mesh_t) :: self
597  integer, intent(in) :: i
598  real(8) :: r
599  !.............................................................................
600  if (i /= self%gn .or. self%gn == 1) then
601  r = self%centre_nat + self%r(i) - 0.5 * self%d
602  else
603  r = self%centre_nat + self%r(i-1) + 0.5 * self%d
604  end if
605 
606 END FUNCTION func_position_f
607 
608 FUNCTION func_radius_c (self, i) RESULT(r)
609  class(mesh_t) :: self
610  integer, intent(in) :: i
611  real(8) :: r
612  !.............................................................................
613  r = abs(func_position_c(self, i))
614 
615 END FUNCTION func_radius_c
616 
617 FUNCTION func_radius_f (self, i) RESULT(r)
618  class(mesh_t) :: self
619  integer, intent(in) :: i
620  real(8) :: r
621  !.............................................................................
622  r = abs(func_position_f(self, i))
623 
624 END FUNCTION func_radius_f
625 
626 FUNCTION func_sintheta_c (self, j) RESULT(sintc)
627  class(mesh_t) :: self
628  integer, intent(in) :: j
629  real(8) :: sintc
630  !.............................................................................
631  sintc = abs(sin(func_position_c(self,j)))
632 
633 END FUNCTION func_sintheta_c
634 
635 FUNCTION func_sintheta_f (self, j) RESULT(sintf)
636  class(mesh_t) :: self
637  integer, intent(in) :: j
638  real(8) :: sintf
639  !.............................................................................
640  sintf = abs(sin(func_position_f(self,j)))
641 
642 END FUNCTION func_sintheta_f
643 
644 !===============================================================================
645 !> Differential volume elements for non-Cartesian grids.
646 !****Note** that the 'f' and 'c' subscripts still denote the centring of a given
647 !> quantity, rather than the centring of the quantities that go into the
648 !> calculation. This is in contrast to ZEUS where the subscript ('a' or 'b') is
649 !> is determined by the centring of the inputs.
650 !===============================================================================
651 FUNCTION func_dr2_c (self, j) RESULT(dvol2c)
652  class(mesh_t) :: self
653  integer, intent(in) :: j
654  real(8) :: dvol2c
655  !.............................................................................
656  dvol2c = 0.5 * abs( self%h32f(j+1) * func_position_f(self,j+1) &
657  - self%h32f(j ) * func_position_f(self,j ))
658 
659 END FUNCTION func_dr2_c
660 
661 FUNCTION func_dr2_f (self, j) RESULT(dvol2f)
662  class(mesh_t) :: self
663  integer, intent(in) :: j
664  real(8) :: dvol2f
665  !.............................................................................
666  dvol2f = 0.5 * abs( self%h32c(j ) * func_position_c(self,j ) &
667  - self%h32c(j-1) * func_position_c(self,j-1))
668 
669 END FUNCTION func_dr2_f
670 
671 FUNCTION func_dr3_c (self, i) RESULT(dvol1c)
672  class(mesh_t) :: self
673  integer, intent(in) :: i
674  real(8) :: dvol1c
675  !.............................................................................
676  dvol1c = self%h2f(i+1) * self%h31f(i+1) * func_position_f(self,i+1) / 3.0d0 &
677  - self%h2f(i ) * self%h31f(i ) * func_position_f(self,i ) / 3.0d0
678 
679 END FUNCTION func_dr3_c
680 
681 FUNCTION func_dr3_f (self, i) RESULT(dvol1f)
682  class(mesh_t) :: self
683  integer, intent(in) :: i
684  real(8) :: dvol1f
685  !.............................................................................
686  dvol1f = self%h2c(i ) * self%h31c(i ) * func_position_c(self,i ) / 3.0d0 &
687  - self%h2c(i-1) * self%h31c(i-1) * func_position_c(self,i-1) / 3.0d0
688 
689 END FUNCTION func_dr3_f
690 
691 FUNCTION func_dcostheta_c (self, j) RESULT(dvol2c)
692  class(mesh_t) :: self
693  integer, intent(in) :: j
694  real(8) :: dvol2c
695  !.............................................................................
696  dvol2c = -cos(func_position_f(self,j+1)) + cos(func_position_f(self,j))
697 
698 END FUNCTION func_dcostheta_c
699 
700 FUNCTION func_dcostheta_f (self, j) RESULT(dvol2f)
701  class(mesh_t) :: self
702  integer, intent(in) :: j
703  real(8) :: dvol2f
704  !.............................................................................
705  dvol2f = -cos(func_position_c(self,j)) + cos(func_position_c(self,j-1))
706 END FUNCTION func_dcostheta_f
707 
708 !===============================================================================
709 !> Secondary area factors; primarily used with operators.
710 !===============================================================================
711 FUNCTION func_dar1_c (self, i) RESULT(dar1c)
712  class(mesh_t) :: self
713  integer, intent(in) :: i
714  real(8) :: dar1c
715  !.............................................................................
716  dar1c = self%h2c(i) * self%h31c(i)
717 
718 END FUNCTION func_dar1_c
719 
720 FUNCTION func_dar1_f (self, i) RESULT(dar1f)
721  class(mesh_t) :: self
722  integer, intent(in) :: i
723  real(8) :: dar1f
724  !.............................................................................
725  dar1f = self%h2f(i) * self%h31f(i)
726 
727 END FUNCTION func_dar1_f
728 
729 FUNCTION func_dar2_c (self,i) RESULT(dar2c)
730  class(mesh_t) :: self
731  integer, intent(in) :: i
732  real(8) :: dar2c
733  !.............................................................................
734  dar2c = self%h31c(i) * self%d / self%dvol1c(i)
735 
736 END FUNCTION func_dar2_c
737 
738 FUNCTION func_dar2_f (self,i) RESULT(dar2f)
739  class(mesh_t) :: self
740  integer, intent(in) :: i
741  real(8) :: dar2f
742  !.............................................................................
743  dar2f = self%h31f(i) * self%d / self%dvol1f(i)
744 
745 END FUNCTION func_dar2_f
746 
747 FUNCTION func_dar31_c (self,i) RESULT(dar31c)
748  class(mesh_t) :: self
749  integer, intent(in) :: i
750  real(8) :: dar31c
751  !.............................................................................
752  dar31c = self%h2c(i) * self%d / self%dvol1c(i)
753 
754 END FUNCTION func_dar31_c
755 
756 FUNCTION func_dar31_f (self,i) RESULT(dar31f)
757  class(mesh_t) :: self
758  integer, intent(in) :: i
759  real(8) :: dar31f
760  !.............................................................................
761  dar31f = self%h2f(i) * self%d / self%dvol1f(i)
762 
763 END FUNCTION func_dar31_f
764 
765 FUNCTION func_dar32_c (self, j) RESULT(dar32c)
766  class(mesh_t) :: self
767  integer, intent(in) :: j
768  real(8) :: dar32c
769  !.............................................................................
770  dar32c = self%d / self%dvol2c(j)
771 
772 END FUNCTION func_dar32_c
773 
774 FUNCTION func_dar32_f (self, j) RESULT(dar32f)
775  class(mesh_t) :: self
776  integer, intent(in) :: j
777  real(8) :: dar32f
778  !.............................................................................
779  dar32f = self%d / self%dvol2f(j)
780 
781 END FUNCTION func_dar32_f
782 
783 !===============================================================================
784 !> Mapping/transformation procedures.
785 !> `OCart` is the origin of the *target* coord. system in *Cartesian* coords.
786 !===============================================================================
787 FUNCTION cartesiantocylindrical4 (x, y, z, OCart) result (vec)
788  real(4), intent(in) :: x, y, z
789  real(4), dimension(3), optional:: OCart
790  real(4), dimension(3) :: vec, lOCart
791  !.............................................................................
792  locart = 0.0
793  if (present(ocart)) locart = ocart
794 
795  vec(1) = z - locart(3) ! Z
796  vec(2) = sqrt((x-locart(1))**2 + (y-locart(2))**2) ! R
797  if (locart(1) == 0.0 .and. locart(2) == 0.0) then
798  vec(3) = atan2(y, x) ! Phi
799  else
800  vec(3) = atan2(y-locart(2), x-locart(1)) ! Phi
801  end if
802  if (vec(3) .lt. 0.0) vec(3) = vec(3) + 2.0 * real(pi,kind=4)
803 
804 END FUNCTION cartesiantocylindrical4
805 
806 FUNCTION cartesiantocylindrical8 (x, y, z, OCart) result (vec)
807  real(8), intent(in) :: x, y, z
808  real(8), dimension(3), optional:: OCart
809  real(8), dimension(3) :: vec, lOCart
810  !.............................................................................
811  locart = 0.0d0
812  if (present(ocart)) locart = ocart
813 
814  vec(1) = z - locart(3) ! Z
815  vec(2) = sqrt((x-locart(1))**2 + (y-locart(2))**2) ! R
816  if (locart(1) == 0.0 .and. locart(2) == 0.0) then
817  vec(3) = atan2(y, x) ! Phi
818  else
819  vec(3) = atan2(y-locart(2), x-locart(1)) ! Phi
820  end if
821  if (vec(3) .lt. 0.0d0) vec(3) = vec(3) + 2.0d0 * pi
822 
823 END FUNCTION cartesiantocylindrical8
824 
825 FUNCTION cartesiantospherical4 (x, y, z, OCart) result (vec)
826  real(4), intent(in) :: x, y, z
827  real(4), dimension(3), optional:: OCart
828  real(4), dimension(3) :: vec, lOCart
829  !.............................................................................
830  locart = 0.0
831  if (present(ocart)) locart = ocart
832 
833  vec(1) = sqrt((x-locart(1))**2 + (y-locart(2))**2 + (z-locart(3))**2) ! r
834  vec(2) = pi - acos((z-locart(3)) / vec(1)) ! theta
835  if (locart(1) == 0.0 .and. locart(2) == 0.0) then
836  vec(3) = -atan2(y, x) ! phi
837  else
838  vec(3) = -atan2((y-locart(2)), (x-locart(1))) ! phi
839  end if
840  if (vec(3) .lt. 0.0) vec(3) = vec(3) + 2.0 * real(pi,kind=4)
841 
842 END FUNCTION cartesiantospherical4
843 
844 FUNCTION cartesiantospherical8 (x, y, z, OCart) result (vec)
845  real(8), intent(in) :: x, y, z
846  real(8), dimension(3), optional:: OCart
847  real(8), dimension(3) :: vec, lOCart
848  !.............................................................................
849  locart = 0.0
850  if (present(ocart)) locart = ocart
851 
852  vec(1) = sqrt((x-locart(1))**2 + (y-locart(2))**2 + (z-locart(3))**2) ! r
853  vec(2) = pi - acos((z-locart(3)) / vec(1)) ! theta
854  if (locart(1) == 0.0 .and. locart(2) == 0.0) then
855  vec(3) = -atan2(y, x) ! phi
856  else
857  vec(3) = -atan2((y-locart(2)), (x-locart(1))) ! phi
858  end if
859  if (vec(3) .lt. 0.0d0) vec(3) = vec(3) + 2.0d0 * pi
860 
861 END FUNCTION cartesiantospherical8
862 
863 FUNCTION cylindricaltocartesian4 (Z, R, Phi, OCart) result (vec)
864  real(4), intent(in) :: Z, R, Phi
865  real(4), dimension(3), optional :: OCart
866  real(4), dimension(3) :: vec, lOCart
867  !.............................................................................
868  locart = 0.0
869  if (present(ocart)) locart = ocart
870 
871  vec(1) = r * cos(phi) - locart(1) ! x
872  vec(2) = r * sin(phi) - locart(2) ! y
873  vec(3) = z - locart(3) ! z
874 
875 END FUNCTION cylindricaltocartesian4
876 
877 FUNCTION cylindricaltocartesian8 (Z, R, Phi, OCart) result (vec)
878  real(8), intent(in) :: Z, R, Phi
879  real(8), dimension(3), optional :: OCart
880  real(8), dimension(3) :: vec, lOCart
881  !.............................................................................
882  locart = 0.0d0
883  if (present(ocart)) locart = ocart
884 
885  vec(1) = r * cos(phi) - locart(1) ! x
886  vec(2) = r * sin(phi) - locart(2) ! y
887  vec(3) = z - locart(3) ! z
888 
889 END FUNCTION cylindricaltocartesian8
890 
891 FUNCTION sphericaltocartesian4 (r, theta, phi, OCart) result (vec)
892  real(4), intent(in) :: r, theta, phi
893  real(4), dimension(3), optional :: OCart
894  real(4), dimension(3) :: vec, lOCart
895  !.............................................................................
896  locart = 0.0
897  if (present(ocart)) locart = ocart
898 
899  vec(1) = r * sin(theta) * cos(phi) - locart(1) ! x
900  vec(2) = -r * sin(theta) * sin(phi) - locart(2) ! y
901  vec(3) = -r * cos(theta) - locart(3) ! z
902 
903 END FUNCTION sphericaltocartesian4
904 
905 FUNCTION sphericaltocartesian8 (r, theta, phi, OCart) result (vec)
906  real(8), intent(in) :: r, theta, phi
907  real(8), dimension(3), optional :: OCart
908  real(8), dimension(3) :: vec, lOCart
909  !.............................................................................
910  locart = 0.0d0
911  if (present(ocart)) locart = ocart
912 
913  vec(1) = r * sin(theta) * cos(phi) - locart(1) ! x
914  vec(2) = -r * sin(theta) * sin(phi) - locart(2) ! y
915  vec(3) = -r * cos(theta) - locart(3) ! z
916 
917 END FUNCTION sphericaltocartesian8
918 
919 FUNCTION sphericaltocylindrical4 (r, theta, phi, OCart) result (vec)
920  real(4), intent(in) :: r, theta, phi
921  real(4), dimension(3), optional :: OCart
922  real(4), dimension(3) :: vec, pCart
923  !.............................................................................
924 
925  ! the most reliable way (that I've found so far) to convert between coord.
926  ! systems that are *not guaranteed to have the same origins* is via Cartesian coords.
927  pcart = sphericaltocartesian4(r, theta, phi, ocart)
928  vec = cartesiantocylindrical4(pcart(1), pcart(2), pcart(3), ocart)
929 
930 END FUNCTION sphericaltocylindrical4
931 
932 FUNCTION sphericaltocylindrical8 (r, theta, phi, OCart) result (vec)
933  real(8), intent(in) :: r, theta, phi
934  real(8), dimension(3), optional :: OCart
935  real(8), dimension(3) :: vec, pCart
936  !.............................................................................
937 
938  ! the most reliable way (that I've found so far) to convert between coord.
939  ! systems that are *not guaranteed to have the same origins* is via Cartesian coords.
940  pcart = sphericaltocartesian8(r, theta, phi, ocart)
941  vec = cartesiantocylindrical8(pcart(1), pcart(2), pcart(3), ocart)
942 
943 END FUNCTION sphericaltocylindrical8
944 
945 FUNCTION cylindricaltospherical4 (Z, R, Phi, OCart) result (vec)
946  real(4), intent(in) :: Z, R, Phi
947  real(4), dimension(3), optional :: OCart
948  real(4), dimension(3) :: vec, pCart
949  !.............................................................................
950 
951  pcart = cylindricaltocartesian4(z, r, phi, ocart)
952  vec = cartesiantospherical4(pcart(1), pcart(2), pcart(3), ocart)
953 
954 END FUNCTION cylindricaltospherical4
955 
956 FUNCTION cylindricaltospherical8 (Z, R, Phi, OCart) result (vec)
957  real(8), intent(in) :: Z, R, Phi
958  real(8), dimension(3), optional :: OCart
959  real(8), dimension(3) :: vec, pCart
960  !.............................................................................
961 
962  pcart = cylindricaltocartesian8(z, r, phi, ocart)
963  vec = cartesiantospherical8(pcart(1), pcart(2), pcart(3), ocart)
964 
965 END FUNCTION cylindricaltospherical8
966 
967 !===============================================================================
968 FUNCTION currenttocartesian4 (mtype, x1, x2, x3) result (vec)
969  integer, intent(in) :: mtype
970  real(4), intent(in) :: x1, x2, x3
971  real(4), dimension(3) :: vec
972  !.............................................................................
973  if (mtype == mesh_types%Cartesian) then
974  ! dummy
975  vec(1) = x1
976  vec(2) = x2
977  vec(3) = x3
978  else if (mtype == mesh_types%spherical) then
979  vec = sphericaltocartesian4(x1, x2, x3)
980  else if (mtype == mesh_types%cylindrical) then
981  vec = cylindricaltocartesian4(x1, x2, x3)
982  end if
983 
984 END FUNCTION currenttocartesian4
985 
986 FUNCTION currenttocartesian8 (mtype, x1, x2, x3) result (vec)
987  integer, intent(in) :: mtype
988  real(8), intent(in) :: x1, x2, x3
989  real(8), dimension(3) :: vec
990  !.............................................................................
991  if (mtype == mesh_types%Cartesian) then
992  ! dummy
993  vec(1) = x1
994  vec(2) = x2
995  vec(3) = x3
996  else if (mtype == mesh_types%spherical) then
997  vec = sphericaltocartesian8(x1, x2, x3)
998  else if (mtype == mesh_types%cylindrical) then
999  vec = cylindricaltocartesian8(x1, x2, x3)
1000  end if
1001 
1002 END FUNCTION currenttocartesian8
1003 
1004 FUNCTION cartesiantocurrent4 (mtype, x1, x2, x3) result (vec)
1005  integer, intent(in) :: mtype
1006  real(4), intent(in) :: x1, x2, x3
1007  real(4), dimension(3) :: vec
1008  !.............................................................................
1009  if (mtype == mesh_types%Cartesian) then
1010  ! dummy
1011  vec(1) = x1
1012  vec(2) = x2
1013  vec(3) = x3
1014  else if (mtype == mesh_types%spherical) then
1015  vec = cartesiantospherical4(x1, x2, x3)
1016  else if (mtype == mesh_types%cylindrical) then
1017  vec = cartesiantocylindrical4(x1, x2, x3)
1018  end if
1019 
1020 END FUNCTION cartesiantocurrent4
1021 
1022 FUNCTION cartesiantocurrent8 (mtype, x1, x2, x3) result (vec)
1023  integer, intent(in) :: mtype
1024  real(8), intent(in) :: x1, x2, x3
1025  real(8), dimension(3) :: vec
1026  !.............................................................................
1027  if (mtype == mesh_types%Cartesian) then
1028  ! dummy
1029  vec(1) = x1
1030  vec(2) = x2
1031  vec(3) = x3
1032  else if (mtype == mesh_types%spherical) then
1033  vec = cartesiantospherical8(x1, x2, x3)
1034  else if (mtype == mesh_types%cylindrical) then
1035  vec = cartesiantocylindrical8(x1, x2, x3)
1036  end if
1037 
1038 END FUNCTION cartesiantocurrent8
1039 
1040 !===============================================================================
1041 !> Calculate the volume centroid of `mesh`.
1042 !> Assumes the mesh type, mesh%idir and mesh%s are already set.
1043 !===============================================================================
1044 FUNCTION volumecentroid (mesh) result (centre)
1045  class(mesh_t) :: mesh
1046  real(4) :: centre
1047  !.............................................................................
1048  call trace_begin('VolumeCentroid')
1049 
1050  associate(s => mesh%s, xmn => mesh%llc_nat)
1051  select type (mesh)
1052  type is (cartesian_mesh)
1053  ! In Cartesian coords., this is simply half the size of the patch.
1054  centre = xmn + 0.5 * s
1055  type is (cylindrical_mesh)
1056  select case (mesh%idir)
1057  case (1,3)
1058  centre = xmn + 0.5 * s
1059  case (2)
1060  ! radial coord.
1061  ! = (2/3) delta(R**3) / delta(R**2)
1062  centre = 2.0 * (3.0 * xmn**2 + 3.0 * xmn * s + s**2) &
1063  / (3.0 * (2.0 * xmn + s))
1064  end select
1065  type is (spherical_mesh)
1066  select case(mesh%idir)
1067  case (1)
1068  ! radial coord.
1069  ! = (3/4) delta(r**4) / delta(r**3)
1070  centre = 3.0 * (4.0 * xmn**3 + 6.0 * xmn**2 * s + 4.0 * xmn * s**2 + s**3) &
1071  / (4.0 * (3.0 * xmn**2 + 3.0 * xmn * s + s**2))
1072  case (2)
1073  ! theta coord.
1074  ! = delta(sin(theta) - theta * cos(theta)) / delta(-cos(theta))
1075  centre = (sin(xmn + s) - (xmn + s) * cos(xmn + s) - sin(xmn) + xmn * cos(xmn)) &
1076  / (cos(xmn) - cos(xmn + s))
1077  case (3)
1078  centre = mesh%llc_nat + 0.5 * s
1079  end select
1080  end select
1081  end associate
1082 
1083  call trace_end
1084 END FUNCTION volumecentroid
1085 
1086 !===============================================================================
1087 !> Round a single-precision floating point number to `digit` sig. figs.
1088 !===============================================================================
1089 FUNCTION round4(a,digit) result (b)
1090  real(4) :: a, b, rof, q1, q2, q3
1091  integer :: digit
1092  real(4) :: tiny = 1.0e-32
1093  integer :: minexp = -6
1094  integer :: i1,i2
1095 ! ----------------------------------------------------------------------
1096 
1097  q1 = abs(a)
1098  if (q1 <= tiny .or. digit == 0.) then
1099  b = a
1100  return
1101  end if
1102  q1 = log10(q1)
1103  q2 = sign(0.5, q1)
1104  i1 = digit - int(q1 + q2 + 0.5)
1105  i2 = min(-3, digit + minexp)
1106  q1 = 10.0**i1
1107  q2 = 10.0**i2
1108  q3 = sign(q2, a)
1109  b = anint(a * q1 + q3) / q1
1110 
1111 END FUNCTION round4
1112 
1113 !===============================================================================
1114 !> Round a double-precision floating point number to `digit` sig. figs.
1115 !===============================================================================
1116 FUNCTION round8(a,digit) result (b)
1117  real(8) :: a, b, rof, q1, q2,q3
1118  integer :: digit
1119  real(8) :: tiny = 1.0d-99
1120  integer :: minexp = -14
1121  integer :: i1,i2
1122 ! ----------------------------------------------------------------------
1123 
1124  q1 = abs(a)
1125  if (q1 <= tiny .or. digit == 0_8) then
1126  b = a
1127  return
1128  end if
1129  q1 = log10(q1)
1130  q2 = sign(0.5_8, q1)
1131  i1 = digit - int(q1 + q2 + 0.5_8)
1132  i2 = min(-3, digit + minexp)
1133  q1 = 10.0_8**i1
1134  q2 = 10.0_8**i2
1135  q3 = sign(q2, a)
1136  b = anint(a * q1 + q3) / q1
1137 
1138 END FUNCTION round8
1139 
1140 END MODULE mesh_mod
Template module for mesh.
Definition: mesh_mod.f90:4
Definition: io_mod.f90:4