Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_mudpack.F90
Go to the documentation of this file.
1 #ifndef DOXYGEN_SHOULD_SKIP_THIS
2 
10 module sll_m_mudpack
11 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12 #include "sll_assert.h"
13 #include "sll_working_precision.h"
14 
15 ! use F77_mudpack, only: &
16 ! mud24cr, &
17 ! mud24sp, &
18 ! mud2cr, &
19 ! mud2sp
20 
21  implicit none
22 
23  public :: &
24  sll_s_delete_mudpack_cartesian, &
25  sll_s_mudpack_cartesian_init, &
26  sll_s_mudpack_polar_init, &
27  sll_o_create, &
28  sll_o_delete, &
29  sll_t_mudpack_solver, &
30  sll_o_solve, &
31  sll_s_solve_mudpack_cartesian, &
32  sll_s_solve_mudpack_polar
33 
34  private
35 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
36 
38  type :: sll_t_mudpack_solver
39 
40  sll_real64, dimension(:), allocatable :: work
41  sll_int32 :: mgopt(4)
42  sll_int32 :: iprm(16)
43  sll_real64 :: fprm(6)
44  sll_int32 :: iguess
45  sll_int32, pointer :: iwork(:, :)
46 
47  end type sll_t_mudpack_solver
48 
49  integer, parameter :: CARTESIAN_2D = 2
50  integer, parameter :: CARTESIAN_3D = 3
51  integer, parameter :: POLAR = 11
52  integer, parameter :: CYLINDRICAL = 12
53 
54  interface sll_o_create
55  module procedure sll_s_mudpack_cartesian_init
56  end interface sll_o_create
57 
58  interface sll_o_solve
59  module procedure sll_s_solve_mudpack_cartesian
60  end interface sll_o_solve
61 
62  interface sll_o_delete
63  module procedure sll_s_delete_mudpack_cartesian
64  end interface sll_o_delete
65 
66 contains
67 
69  subroutine sll_s_mudpack_cartesian_init(self, &
70  eta1_min, eta1_max, nc_eta1, &
71  eta2_min, eta2_max, nc_eta2, &
72  bc_eta1_left, bc_eta1_right, &
73  bc_eta2_left, bc_eta2_right)
74  implicit none
75 
76 ! set grid size params
77  type(sll_t_mudpack_solver):: self
78  sll_real64, intent(in) :: eta1_min
79  sll_real64, intent(in) :: eta1_max
80  sll_real64, intent(in) :: eta2_min
81  sll_real64, intent(in) :: eta2_max
82  sll_int32, intent(in) :: nc_eta1
83  sll_int32, intent(in) :: nc_eta2
84  sll_int32, optional :: bc_eta1_left
85  sll_int32, optional :: bc_eta1_right
86  sll_int32, optional :: bc_eta2_left
87  sll_int32, optional :: bc_eta2_right
88  sll_int32, parameter :: iixp = 4, jjyq = 4
89  sll_int32 :: iiex, jjey, llwork
90 
91  sll_real64 :: phi(nc_eta1 + 1, nc_eta2 + 1)
92  sll_real64 :: rhs(nc_eta1 + 1, nc_eta2 + 1)
93 
94 !put integer and floating point argument names in contiguous
95 !storeage for labelling in vectors iprm,fprm
96  sll_int32 :: iprm(16)
97  sll_real64 :: fprm(6)
98  sll_int32 :: error
99  sll_int32 :: intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny
100  sll_int32 :: iguess, maxcy, method, nwork, lwrkqd, itero
101  common/itmud2sp/intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny, &
102  iguess, maxcy, method, nwork, lwrkqd, itero
103  sll_real64 :: xa, xb, yc, yd, tolmax, relmax
104  common/ftmud2sp/xa, xb, yc, yd, tolmax, relmax
105 !sll_real64,dimension(:),allocatable :: cxx_array
106 
107 #ifdef DEBUG
108  sll_int32 :: i
109 #endif
110 
111  equivalence(intl, iprm)
112  equivalence(xa, fprm)
113 
114  nx = nc_eta1 + 1
115  ny = nc_eta2 + 1
116 
117 ! set minimum required work space
118  llwork = (7*(nx + 2)*(ny + 2) + 44*nx*ny)/3
119 
120  allocate (self%work(llwork))
121  iiex = ceiling(log((nx - 1.)/iixp)/log(2.)) + 1
122  jjey = ceiling(log((ny - 1.)/jjyq)/log(2.)) + 1
123 
124 !set input integer arguments
125  intl = 0
126 
127 !set boundary condition flags
128  nxa = merge(bc_eta1_left, 0, present(bc_eta1_left))
129  nxb = merge(bc_eta1_right, 0, present(bc_eta1_right))
130  nyc = merge(bc_eta2_left, 0, present(bc_eta2_left))
131  nyd = merge(bc_eta2_right, 0, present(bc_eta2_right))
132 
133 !set grid sizes from parameter statements
134  ixp = iixp
135  jyq = jjyq
136  iex = iiex
137  jey = jjey
138 
139  nx = ixp*(2**(iex - 1)) + 1
140  ny = jyq*(2**(jey - 1)) + 1
141 
142  if (nx /= nc_eta1 + 1 .or. ny /= nc_eta2 + 1) then
143  print *, "nx,nc_eta1+1=", nx, nc_eta1 + 1
144  print *, "ny,nc_eta2+1=", ny, nc_eta2 + 1
145  stop ' nx or ny different in sll_mudpack_cartesian '
146  end if
147 
148 !set multigrid arguments (w(2,1) cycling with fully weighted
149 !residual restriction and cubic prolongation)
150  self%mgopt(1) = 2
151  self%mgopt(2) = 2
152  self%mgopt(3) = 1
153  self%mgopt(4) = 3
154 
155 !set for three cycles to ensure second-order approximation is computed
156  maxcy = 5
157 
158 !set no initial guess forcing full multigrid cycling
159  self%iguess = 0
160  iguess = self%iguess
161 
162 !set work space length approximation from parameter statement
163  nwork = llwork
164 
165 !set point relaxation
166  method = 0
167 
168 !set end points of solution rectangle in (x,y) space
169  xa = eta1_min
170  xb = eta1_max
171  yc = eta2_min
172  yd = eta2_max
173 
174 !set for no error control flag
175  tolmax = 0.0_f64
176 
177 #ifdef DEBUG
178  write (*, 101) (iprm(i), i=1, 15)
179  write (*, 102) (self%mgopt(i), i=1, 4)
180  write (*, 103) xa, xb, yc, yd, tolmax
181  write (*, 104) intl
182 #endif
183 
184  call mud2sp(iprm, fprm, self%work, cofx, cofy, bndsp, rhs, phi, self%mgopt, error)
185 
186 #ifdef DEBUG
187 !print error parameter and minimum work space requirement
188  write (*, 200) error, iprm(16)
189  if (error > 0) stop 0
190 
191 101 format(/'# integer input arguments ', &
192  /'#intl = ', i2, ' nxa = ', i2, ' nxb = ', i2, ' nyc = ', i2, ' nyd = ', i2, &
193  /'# ixp = ', i2, ' jyq = ', i2, ' iex = ', i2, ' jey = ', i2 &
194  /'# nx = ', i3, ' ny = ', i3, ' iguess = ', i2, ' maxcy = ', i2, &
195  /'# method = ', i2, ' work space estimate = ', i7)
196 
197 102 format(/'# multigrid option arguments ', &
198  /'# kcycle = ', i2, &
199  /'# iprer = ', i2, &
200  /'# ipost = ', i2, &
201  /'# intpol = ', i2)
202 
203 103 format(/'# floating point input parameters ', &
204  /'# xa = ', f6.3, ' xb = ', f6.3, ' yc = ', f6.3, ' yd = ', f6.3, &
205  /'# tolerance (error control) = ', e10.3)
206 
207 104 format(/'# discretization call to mud2sp', ' intl = ', i2)
208 
209 200 format('# error = ', i2, ' minimum work space = ', i7)
210 
211 #endif
212 
213  end subroutine sll_s_mudpack_cartesian_init
214 
216  subroutine sll_s_solve_mudpack_cartesian(self, phi, rhs, ex, ey, nrj)
217 
218  type(sll_t_mudpack_solver) :: self
219  sll_real64 :: phi(:, :)
220  sll_real64 :: rhs(:, :)
221  sll_real64, optional :: ex(:, :)
222  sll_real64, optional :: ey(:, :)
223  sll_real64, optional :: nrj
224 
225 !put integer and floating point argument names in contiguous
226 !storeage for labelling in vectors iprm,fprm
227  sll_int32 :: iprm(16)
228  sll_real64 :: fprm(6)
229  sll_int32 :: error
230  sll_int32 :: i, j
231  sll_real64 :: dx, dy
232 
233  sll_int32 :: intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny
234  sll_int32 :: iguess, maxcy, method, nwork, lwrkqd, itero
235  common/itmud2sp/intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny, &
236  iguess, maxcy, method, nwork, lwrkqd, itero
237  sll_real64 :: xa, xb, yc, yd, tolmax, relmax
238  common/ftmud2sp/xa, xb, yc, yd, tolmax, relmax
239 
240  equivalence(intl, iprm)
241  equivalence(xa, fprm)
242 
243 !set initial guess because solve should be called every time step in a
244 !time dependent problem and the elliptic operator does not depend on time.
245  if (self%iguess == 0) then
246  iguess = 0
247  else
248  self%iguess = 1
249  iguess = 1
250  end if
251 
252 !attempt solution
253  intl = 1
254 #ifdef DEBUG
255  write (*, 106) intl, method, iguess
256 #endif
257 
258  if (nxa == 0 .and. nyc == 0) &
259  rhs = rhs - sum(rhs)/real(nx*ny, f64)
260 
261  call mud2sp(iprm, fprm, self%work, cofx, cofy, bndsp, rhs, phi, self%mgopt, error)
262 
263 #ifdef DEBUG
264  write (*, 107) error
265  if (error > 0) stop 0
266 #endif
267 
268  if (nxa == 0 .and. nyc == 0) &
269  phi = phi - sum(phi)/real(nx*ny, f64)
270 
271  iguess = 1
272 ! attempt to improve approximation to fourth order
273  call mud24sp(self%work, phi, error)
274 
275 #ifdef DEBUG
276  write (*, 108) error
277  if (error > 0) stop 0
278 
279 106 format(/'#approximation call to mud2sp', &
280  /'# intl = ', i2, ' method = ', i2, ' iguess = ', i2)
281 107 format('#error = ', i2)
282 108 format(/'# mud24sp test ', ' error = ', i2)
283 
284 #endif
285 
286  if (present(ex) .and. present(ey)) then
287  dx = (xb - xa)/(nx - 1)
288  dy = (yd - yc)/(ny - 1)
289  do i = 2, nx - 1
290  ex(i, :) = (phi(i + 1, :) - phi(i - 1, :))/(2*dx)
291  end do
292  do j = 2, ny - 1
293  ey(:, j) = (phi(:, j + 1) - phi(:, j - 1))/(2*dy)
294  end do
295 
296  if (nxa == 0) then
297  ex(nx, :) = (phi(2, :) - phi(nx - 1, :))/(2*dx)
298  ex(1, :) = ex(nx, :)
299  end if
300  if (nyc == 0) then
301  ey(:, ny) = (phi(:, 2) - phi(:, ny - 1))/(2*dy)
302  ey(:, 1) = ey(:, ny)
303  end if
304 
305  if (present(nrj)) then
306  nrj = sum(ex*ex + ey*ey)*dx*dy
307  end if
308 
309  end if
310 
311  end subroutine sll_s_solve_mudpack_cartesian
312 
314  subroutine sll_s_delete_mudpack_cartesian(self)
315  type(sll_t_mudpack_solver) :: self
316 
317  deallocate (self%work)
318 
319  end subroutine sll_s_delete_mudpack_cartesian
320 
323  subroutine sll_s_mudpack_polar_init(self, &
324  r_min, r_max, nr, &
325  theta_min, theta_max, nth, &
326  bc_r_min, bc_r_max, &
327  bc_theta_min, bc_theta_max)
328  implicit none
329 
330  type(sll_t_mudpack_solver) :: self
331  sll_real64, intent(in) :: r_min
332  sll_real64, intent(in) :: r_max
333  sll_real64, intent(in) :: theta_min
334  sll_real64, intent(in) :: theta_max
335  sll_int32, intent(in) :: nr
336  sll_int32, intent(in) :: nth
337  sll_int32 :: icall
338 !sll_int32 :: iiex,jjey
339  sll_int32 :: llwork
340  sll_int32 :: bc_r_min
341  sll_int32 :: bc_r_max
342  sll_int32 :: bc_theta_min
343  sll_int32 :: bc_theta_max
344 
345  sll_real64 :: phi(nr, nth)
346  sll_real64 :: rhs(nr, nth)
347 
348 ! put sll_int32 and floating point argument names in contiguous
349 ! storeage for labelling in vectors iprm,fprm
350 
351  sll_int32 :: intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny, &
352  iguess, maxcy, method, nwork, lwrkqd, itero
353  common/itmud2cr/intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny, &
354  iguess, maxcy, method, nwork, lwrkqd, itero
355  sll_real64 :: xa, xb, yc, yd, tolmax, relmax
356  common/ftmud2cr/xa, xb, yc, yd, tolmax, relmax
357  sll_int32 :: i
358 !sll_int32 :: j
359  sll_int32 :: ierror
360  sll_int32 :: iprm(16)
361  sll_real64 :: fprm(6)
362 
363  equivalence(intl, iprm)
364  equivalence(xa, fprm)
365 
366  nx = nr
367  ny = nth
368 
369 ! set minimum required work space
370  llwork = (7*(nx + 2)*(ny + 2) + 44*nx*ny)/3
371 
372  allocate (self%work(llwork))
373  icall = 0
374 
375 ! set input sll_int32 arguments
376  intl = 0
377 
378 ! set boundary condition flags
379  nxa = bc_r_min
380  nxb = bc_r_max
381  nyc = bc_theta_min
382  nyd = bc_theta_max
383 
384 ! set grid sizes from parameter statements
385  ixp = 2
386  jyq = 2
387  iex = ceiling(log((nx - 1.)/ixp)/log(2.)) + 1
388  jey = ceiling(log((ny - 1.)/jyq)/log(2.)) + 1
389 
390  nx = ixp*(2**(iex - 1)) + 1
391  ny = jyq*(2**(jey - 1)) + 1
392 
393  if (nx /= nr) then
394  print *, "nx,nr=", nx, nr
395  stop ' nx different de nr dans mg_polar_poisson'
396  end if
397  if (ny /= nth) then
398  print *, "ny,nth=", ny, nth
399  stop ' ny different de nth dans mg_polar_poisson'
400  end if
401 
402 ! set multigrid arguments (w(2,1) cycling with fully weighted
403 ! residual restriction and cubic prolongation)
404  self%mgopt(1) = 2
405  self%mgopt(2) = 2
406  self%mgopt(3) = 1
407  self%mgopt(4) = 3
408 
409 ! set three cycles to ensure second-order approx
410  maxcy = 3
411 
412 ! set no initial guess forcing full multigrid cycling
413  iguess = 0
414 
415 ! set work space length approximation from parameter statement
416  nwork = llwork
417 
418 ! set point relaxation
419  method = 0
420 
421 ! set mesh increments
422  xa = r_min
423  xb = r_max
424  yc = theta_min
425  yd = theta_max
426 
427 ! set for no error control flag
428  tolmax = 0.0_f64
429 
430  write (*, 100)
431  write (*, 101) (iprm(i), i=1, 15)
432  write (*, 102) (self%mgopt(i), i=1, 4)
433  write (*, 103) xa, xb, yc, yd, tolmax
434  write (*, 104) intl
435 
436  call mud2cr(iprm, fprm, self%work, coef_polar, bndcr, rhs, phi, self%mgopt, ierror)
437 
438  write (*, 200) ierror, iprm(16)
439  if (ierror > 0) stop 0
440 
441 100 format(//' multigrid poisson solver in polar coordinates ')
442 101 format(/' integer input arguments ', &
443  /'intl = ', i2, ' nxa = ', i2, ' nxb = ', i2, ' nyc = ', i2, ' nyd = ', i2, &
444  /' ixp = ', i2, ' jyq = ', i2, ' iex = ', i2, ' jey = ', i2 &
445  /' nx = ', i3, ' ny = ', i3, ' iguess = ', i2, ' maxcy = ', i2, &
446  /' method = ', i2, ' work space estimate = ', i7)
447 102 format(/' multigrid option arguments ', &
448  /' kcycle = ', i2, &
449  /' iprer = ', i2, &
450  /' ipost = ', i2 &
451  /' intpol = ', i2)
452 103 format(/' floating point input parameters ', &
453  /' xa = ', f6.3, ' xb = ', f6.3, ' yc = ', f6.3, ' yd = ', f6.3, &
454  /' tolerance (error control) = ', e10.3)
455 104 format(/' discretization call to mud2cr', ' intl = ', i2)
456 200 format(' ierror = ', i2, ' minimum work space = ', i7)
457 
458  return
459  end subroutine sll_s_mudpack_polar_init
460 
462  subroutine sll_s_solve_mudpack_polar(self, phi, rhs)
463  implicit none
464 
465 ! set grid size params
466  type(sll_t_mudpack_solver) :: self
467  sll_int32 :: icall
468  sll_int32, parameter :: iixp = 2, jjyq = 2
469 
470  sll_real64, intent(inout) :: phi(:, :)
471  sll_real64, intent(inout) :: rhs(:, :)
472 
473 ! put sll_int32 and floating point argument names in contiguous
474 ! storeage for labelling in vectors iprm,fprm
475 
476  sll_int32 :: intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny
477  sll_int32 :: iguess, maxcy, method, nwork, lwrkqd, itero
478  sll_real64 :: xa, xb, yc, yd, tolmax, relmax
479  sll_int32 :: ierror
480  sll_int32 :: iprm(16)
481  sll_real64 :: fprm(6)
482 
483  common/itmud2cr/intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny, &
484  iguess, maxcy, method, nwork, lwrkqd, itero
485  common/ftmud2cr/xa, xb, yc, yd, tolmax, relmax
486 
487  equivalence(intl, iprm)
488  equivalence(xa, fprm)
489 
490  icall = 1
491  intl = 1
492  write (*, 106) intl, method, iguess
493 ! attempt solution
494  call mud2cr(iprm, fprm, self%work, coef_polar, bndcr, rhs, phi, self%mgopt, ierror)
495  sll_assert(ierror == 0)
496 ! attempt fourth order approximation
497  call mud24cr(self%work, coef_polar, bndcr, phi, ierror)
498  sll_assert(ierror == 0)
499 
500 106 format(/' approximation call to mud2cr', &
501  /' intl = ', i2, ' method = ', i2, ' iguess = ', i2)
502 
503  end subroutine sll_s_solve_mudpack_polar
504 
507  subroutine coef_polar(x, y, cxx, cxy, cyy, cx, cy, ce)
508  real(8) :: x, y, cxx, cxy, cyy, cx, cy, ce
509  cxx = 1.0_8 + 0.0_8*y
510  cxy = 0.0_8
511  cyy = 1.0_8/(x*x)
512  cx = 1.0_8/x
513  cy = 0.0_8
514  ce = 0.0_8
515  end subroutine coef_polar
516 
519  subroutine bndcr(kbdy, xory, alfa, beta, gama, gbdy)
520  integer :: kbdy
521  real(8) :: xory, alfa, beta, gama, gbdy
522 
523  if (kbdy .eq. 2) then
524 
525  ! x=xb boundary.
526  ! b.c. has the form alfyd(x)*px+betyd(x)*py+gamyd(x)*pe = gbdyd(x)
527  ! where x = yorx. alfa,beta,gama,gbdy corresponding to alfyd(x),
528  ! betyd(x),gamyd(x),gbdyd(y) must be output.
529 
530  alfa = 1.0_8 + 0.0_8*xory
531  beta = 0.0_8
532  gama = 0.0_8
533  gbdy = 0.0_8
534 
535  end if
536 
537  end subroutine bndcr
538 
539 !the form of the pde solved is:
540 !
541 !
542 ! cxx(x)*pxx + cx(x)*px + cex(x)*p(x,y) +
543 !
544 ! cyy(y)*pyy + cy(y)*py + cey(y)*p(x,y) = r(x,y)
545 !
546 ! pxx,pyy,px,py are second and first partial derivatives of the
547 ! unknown real solution function p(x,y) with respect to the
548 ! independent variables x,y. cxx,cx,cex,cyy,cy,cey are the known
549 ! real coefficients of the elliptic pde and r(x,y) is the known
550 ! real right hand side of the equation. cxx and cyy should be
551 ! positive for all x,y in the solution region. If some of the
552 ! coefficients depend on both x and y then the PDE is nonseparable.
553 
555  subroutine cofx(x, cxx, cx, cex)
556  real(8) :: x, cxx, cx, cex
557  cxx = 1.0_8 + 0.0_8*x
558  cx = 0.0_8
559  cex = 0.0_8
560  end subroutine cofx
561 
563  subroutine cofy(y, cyy, cy, cey)
564  real(8) :: y, cyy, cy, cey
565  cyy = 1.0_8 + 0.0_8*y
566  cy = 0.0_8
567  cey = 0.0_8
568  end subroutine cofy
569 
571  subroutine bndsp(kbdy, xory, alfa, gbdy)
572  integer :: kbdy
573  real(8) :: xory, alfa, gbdy, x, y, pe, px, py
574  real(8) :: xa, xb, yc, yd, tolmax, relmax
575  common/ftmud2sp/xa, xb, yc, yd, tolmax, relmax
576 
577 !subroutine not used in periodic case
578  if (kbdy == 1) then ! x=xa boundary
579  y = xory
580  x = xa
581  alfa = -1.0_8
582  gbdy = px + alfa*pe
583  return
584  end if
585 
586  if (kbdy == 4) then ! y=yd boundary
587  y = yd
588  x = xory
589  alfa = 1.0_8
590  gbdy = py + alfa*pe
591  return
592  end if
593  end subroutine bndsp
594 
595 end module sll_m_mudpack
596 
597 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1