Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_poisson_2d_mudpack_curvilinear.F90
Go to the documentation of this file.
1 !**************************************************************
2 !**************************************************************
3 ! Copyright INRIA
4 ! Authors :
5 ! CALVI project team
6 !
7 ! This code SeLaLib (for Semi-Lagrangian-Library)
8 ! is a parallel library for simulating the plasma turbulence
9 ! in a tokamak.
10 !
11 ! This software is governed by the CeCILL-B license
12 ! under French law and abiding by the rules of distribution
13 ! of free software. You can use, modify and redistribute
14 ! the software under the terms of the CeCILL-B license as
15 ! circulated by CEA, CNRS and INRIA at the following URL
16 ! "http://www.cecill.info".
17 !**************************************************************
21 !**************************************************************
22 
45 module sll_m_poisson_2d_mudpack_curvilinear
46 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
47 #include "sll_memory.h"
48 #include "sll_working_precision.h"
49 
50  use sll_m_boundary_condition_descriptors, only: &
51  sll_p_dirichlet, &
52  sll_p_periodic
53 
54  use sll_m_cubic_spline_interpolator_1d, only: &
55  sll_f_new_cubic_spline_interpolator_1d
56 
57  use sll_m_cubic_spline_interpolator_2d, only: &
58  sll_f_new_cubic_spline_interpolator_2d
59 
60  use sll_m_interpolators_1d_base, only: &
61  sll_c_interpolator_1d
62 
63  use sll_m_interpolators_2d_base, only: &
64  sll_c_interpolator_2d
65 
66  use sll_m_mudpack_curvilinear, only: &
67  sll_p_non_separable_with_cross_terms, &
68  sll_p_non_separable_without_cross_terms, &
69  sll_p_separable
70 
71  use sll_m_poisson_2d_base, only: &
72  sll_c_poisson_2d_base, &
73  sll_i_function_of_position
74 
75  implicit none
76 
77  public :: &
78  sll_f_new_poisson_2d_mudpack_curvilinear
79 
80  private
81 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
82 
84  type, extends(sll_c_poisson_2d_base) :: poisson_2d_mudpack_curvilinear
85 
87  sll_real64, dimension(:, :), pointer :: cxx_2d
89  sll_real64, dimension(:, :), pointer :: cxy_2d
91  sll_real64, dimension(:, :), pointer :: cyy_2d
93  sll_real64, dimension(:, :), pointer :: cx_2d
95  sll_real64, dimension(:, :), pointer :: cy_2d
97  sll_real64, dimension(:, :), pointer :: ce_2d
99  sll_real64, dimension(:), pointer :: cxx_1d
101  sll_real64, dimension(:), pointer :: cyy_1d
103  sll_real64, dimension(:), pointer :: cx_1d
105  sll_real64, dimension(:), pointer :: cy_1d
107  sll_real64, dimension(:), pointer :: cex_1d
109  sll_real64, dimension(:), pointer :: cey_1d
111  sll_real64 :: cxx
113  sll_real64 :: cyy
115  sll_real64 :: cx
117  sll_real64 :: cy
119  sll_real64 :: ce
121  sll_int32 :: mudpack_case
123  class(sll_c_interpolator_2d), pointer :: cxx_2d_interp
125  class(sll_c_interpolator_2d), pointer :: cxy_2d_interp
127  class(sll_c_interpolator_2d), pointer :: cyy_2d_interp
129  class(sll_c_interpolator_2d), pointer :: cx_2d_interp
131  class(sll_c_interpolator_2d), pointer :: cy_2d_interp
133  class(sll_c_interpolator_2d), pointer :: ce_2d_interp
135  class(sll_c_interpolator_1d), pointer :: cxx_1d_interp
137  class(sll_c_interpolator_1d), pointer :: cyy_1d_interp
139  class(sll_c_interpolator_1d), pointer :: cx_1d_interp
141  class(sll_c_interpolator_1d), pointer :: cy_1d_interp
143  class(sll_c_interpolator_1d), pointer :: cex_1d_interp
145  class(sll_c_interpolator_1d), pointer :: cey_1d_interp
146 
147  sll_real64, dimension(:), pointer :: work
148  sll_int32 :: mgopt(4)
149  sll_int32 :: iprm(16)
150  sll_real64 :: fprm(6)
151  sll_int32 :: iguess
152 
153  contains
154 
156  procedure, pass(poisson) :: initialize => initialize_poisson_2d_mudpack_curvilinear
158  procedure, pass(poisson) :: compute_phi_from_rho => compute_phi_from_rho_2d_mudpack
160  procedure, pass(poisson) :: compute_e_from_rho => compute_e_from_rho_2d_mudpack
161 
163  procedure :: &
164  l2norm_squared => l2norm_squarred_2d_mudpack
166  procedure :: &
167  compute_rhs_from_function => compute_rhs_from_function_2d_mudpack
169  procedure :: free => delete_2d_mudpack
170 
171  end type poisson_2d_mudpack_curvilinear
172 
174  class(poisson_2d_mudpack_curvilinear), pointer :: mudpack_wrapper => null()
175 
176 contains
177 
179  function sll_f_new_poisson_2d_mudpack_curvilinear( &
180  eta1_min, &
181  eta1_max, &
182  nc_eta1, &
183  eta2_min, &
184  eta2_max, &
185  nc_eta2, &
186  bc_eta1_left, &
187  bc_eta1_right, &
188  bc_eta2_left, &
189  bc_eta2_right, &
190  mudpack_case, &
191  cxx_2d, &
192  cxy_2d, &
193  cyy_2d, &
194  cx_2d, &
195  cy_2d, &
196  ce_2d, &
197  cxx_1d, &
198  cyy_1d, &
199  cx_1d, &
200  cy_1d, &
201  cex_1d, &
202  cey_1d, &
203  cxx, &
204  cxy, &
205  cyy, &
206  cx, &
207  cy, &
208  ce) &
209  result(poisson)
210 
211  type(poisson_2d_mudpack_curvilinear), pointer :: poisson
212  sll_real64, intent(in) :: eta1_min
213  sll_real64, intent(in) :: eta1_max
214  sll_int32, intent(in) :: nc_eta1
215  sll_real64, intent(in) :: eta2_min
216  sll_real64, intent(in) :: eta2_max
217  sll_int32, intent(in) :: nc_eta2
218  sll_int32, intent(in) :: bc_eta1_left
219  sll_int32, intent(in) :: bc_eta1_right
220  sll_int32, intent(in) :: bc_eta2_left
221  sll_int32, intent(in) :: bc_eta2_right
222  sll_int32, intent(in) :: mudpack_case
223  sll_real64, dimension(:, :), intent(in), optional :: cxx_2d
224  sll_real64, dimension(:, :), intent(in), optional :: cxy_2d
225  sll_real64, dimension(:, :), intent(in), optional :: cyy_2d
226  sll_real64, dimension(:, :), intent(in), optional :: cx_2d
227  sll_real64, dimension(:, :), intent(in), optional :: cy_2d
228  sll_real64, dimension(:, :), intent(in), optional :: ce_2d
229  sll_real64, dimension(:), intent(in), optional :: cxx_1d
230  sll_real64, dimension(:), intent(in), optional :: cyy_1d
231  sll_real64, dimension(:), intent(in), optional :: cx_1d
232  sll_real64, dimension(:), intent(in), optional :: cy_1d
233  sll_real64, dimension(:), intent(in), optional :: cex_1d
234  sll_real64, dimension(:), intent(in), optional :: cey_1d
235  sll_real64, intent(in), optional :: cxx
236  sll_real64, intent(in), optional :: cxy
237  sll_real64, intent(in), optional :: cyy
238  sll_real64, intent(in), optional :: cx
239  sll_real64, intent(in), optional :: cy
240  sll_real64, intent(in), optional :: ce
241 
242  sll_int32 :: ierr
243 
244  sll_allocate(poisson, ierr)
245  call initialize_poisson_2d_mudpack_curvilinear( &
246  poisson, &
247  eta1_min, &
248  eta1_max, &
249  nc_eta1, &
250  eta2_min, &
251  eta2_max, &
252  nc_eta2, &
253  bc_eta1_left, &
254  bc_eta1_right, &
255  bc_eta2_left, &
256  bc_eta2_right, &
257  mudpack_case, &
258  cxx_2d, &
259  cxy_2d, &
260  cyy_2d, &
261  cx_2d, &
262  cy_2d, &
263  ce_2d, &
264  cxx_1d, &
265  cyy_1d, &
266  cx_1d, &
267  cy_1d, &
268  cex_1d, &
269  cey_1d, &
270  cxx, &
271  cxy, &
272  cyy, &
273  cx, &
274  cy, &
275  ce)
276 
277  end function sll_f_new_poisson_2d_mudpack_curvilinear
278 
279  subroutine initialize_poisson_2d_mudpack_curvilinear( &
280  poisson, &
281  eta1_min, &
282  eta1_max, &
283  nc_eta1, &
284  eta2_min, &
285  eta2_max, &
286  nc_eta2, &
287  bc_eta1_left, &
288  bc_eta1_right, &
289  bc_eta2_left, &
290  bc_eta2_right, &
291  mudpack_case, &
292  cxx_2d, &
293  cxy_2d, &
294  cyy_2d, &
295  cx_2d, &
296  cy_2d, &
297  ce_2d, &
298  cxx_1d, &
299  cyy_1d, &
300  cx_1d, &
301  cy_1d, &
302  cex_1d, &
303  cey_1d, &
304  cxx, &
305  cxy, &
306  cyy, &
307  cx, &
308  cy, &
309  ce)
310  class(poisson_2d_mudpack_curvilinear), target :: poisson
311  sll_real64, intent(in) :: eta1_min
312  sll_real64, intent(in) :: eta1_max
313  sll_int32, intent(in) :: nc_eta1
314  sll_real64, intent(in) :: eta2_min
315  sll_real64, intent(in) :: eta2_max
316  sll_int32, intent(in) :: nc_eta2
317  sll_int32, intent(in) :: bc_eta1_left
318  sll_int32, intent(in) :: bc_eta1_right
319  sll_int32, intent(in) :: bc_eta2_left
320  sll_int32, intent(in) :: bc_eta2_right
321  sll_int32, intent(in) :: mudpack_case
322  sll_real64, dimension(:, :), intent(in), optional :: cxx_2d
323  sll_real64, dimension(:, :), intent(in), optional :: cxy_2d
324  sll_real64, dimension(:, :), intent(in), optional :: cyy_2d
325  sll_real64, dimension(:, :), intent(in), optional :: cx_2d
326  sll_real64, dimension(:, :), intent(in), optional :: cy_2d
327  sll_real64, dimension(:, :), intent(in), optional :: ce_2d
328  sll_real64, dimension(:), intent(in), optional :: cxx_1d
329  sll_real64, dimension(:), intent(in), optional :: cyy_1d
330  sll_real64, dimension(:), intent(in), optional :: cx_1d
331  sll_real64, dimension(:), intent(in), optional :: cy_1d
332  sll_real64, dimension(:), intent(in), optional :: cex_1d
333  sll_real64, dimension(:), intent(in), optional :: cey_1d
334  sll_real64, intent(in), optional :: cxx
335  sll_real64, intent(in), optional :: cxy
336  sll_real64, intent(in), optional :: cyy
337  sll_real64, intent(in), optional :: cx
338  sll_real64, intent(in), optional :: cy
339  sll_real64, intent(in), optional :: ce
340  sll_int32 :: ierr
341 
342  !!!! begin variables for mudpack
343  sll_int32, parameter :: iixp = 2, jjyq = 2
344  sll_int32 :: icall, iiex, jjey, llwork
345  sll_real64, pointer :: phi(:)
346  sll_real64, pointer :: rhs(:)
347  !put integer and floating point argument names in contiguous
348  !storeage for labelling in vectors iprm,fprm
349  sll_int32 :: iprm(16)
350  sll_real64 :: fprm(6)
351  !sll_int32 :: i
352  sll_int32 :: error
353  sll_int32 :: intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny
354  sll_int32 :: iguess, maxcy, method, nwork, lwrkqd, itero
355  common/itmud2sp/intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny, &
356  iguess, maxcy, method, nwork, lwrkqd, itero
357  sll_real64 :: xa, xb, yc, yd, tolmax, relmax
358  common/ftmud2sp/xa, xb, yc, yd, tolmax, relmax
359  equivalence(intl, iprm)
360  equivalence(xa, fprm)
361  !!!! end variables for mudpack
362 
363  nx = nc_eta1 + 1
364  ny = nc_eta2 + 1
365  allocate (phi(nx*ny))
366  allocate (rhs(nx*ny))
367  ! set minimum required work space
368  llwork = (7*(nx + 2)*(ny + 2) + 44*nx*ny)/3
369 
370  allocate (poisson%work(llwork))
371  icall = 0
372  iiex = ceiling(log((nx - 1.)/iixp)/log(2.)) + 1
373  jjey = ceiling(log((ny - 1.)/jjyq)/log(2.)) + 1
374 
375  !set input integer arguments
376  intl = 0
377 
378  !set boundary condition flags
379  nxa = bc_eta1_left
380  nxb = bc_eta1_right
381  nyc = bc_eta2_left
382  nyd = bc_eta2_right
383 
384  !set grid sizes from parameter statements
385  ixp = iixp
386  jyq = jjyq
387  iex = iiex
388  jey = jjey
389 
390  nx = ixp*(2**(iex - 1)) + 1
391  ny = jyq*(2**(jey - 1)) + 1
392 
393  if (nx /= nc_eta1 + 1 .or. ny /= nc_eta2 + 1) then
394  print *, "nx,nc_eta1+1=", nx, nc_eta1 + 1
395  print *, "ny,nc_eta2+1=", ny, nc_eta2 + 1
396  stop ' nx or ny different in sll_mudpack_cartesian '
397  end if
398 
399  !set multigrid arguments (w(2,1) cycling with fully weighted
400  !residual restriction and cubic prolongation)
401  poisson%mgopt(1) = 2
402  poisson%mgopt(2) = 2
403  poisson%mgopt(3) = 1
404  poisson%mgopt(4) = 3
405 
406  !set for three cycles to ensure second-order approximation is computed
407  maxcy = 3
408 
409  !set no initial guess forcing full multigrid cycling
410  poisson%iguess = 0
411  iguess = poisson%iguess
412 
413  !set work space length approximation from parameter statement
414  nwork = llwork
415 
416  !set point relaxation
417  method = 0
418 
419  !set end points of solution rectangle in (x,y) space
420  xa = eta1_min
421  xb = eta1_max
422  yc = eta2_min
423  yd = eta2_max
424 
425  !set for no error control flag
426  tolmax = 0.0_f64
427 
428 ! write(*,101) (iprm(i),i=1,15)
429 ! write(*,102) (poisson%mgopt(i),i=1,4)
430 ! write(*,103) xa,xb,yc,yd,tolmax
431 ! write(*,104) intl
432 
433 !call mud2sp(iprm,fprm,self%work,cofx,cofy,bndsp,rhs,phi,self%mgopt,error)
434 
435  poisson%mudpack_case = mudpack_case
436 
437  poisson%cxx_2d_interp => null()
438  poisson%cxy_2d_interp => null()
439  poisson%cyy_2d_interp => null()
440  poisson%cx_2d_interp => null()
441  poisson%cy_2d_interp => null()
442  poisson%ce_2d_interp => null()
443 
444  poisson%cxx_1d_interp => null()
445  poisson%cyy_1d_interp => null()
446  poisson%cx_1d_interp => null()
447  poisson%cy_1d_interp => null()
448  poisson%cex_1d_interp => null()
449  poisson%cey_1d_interp => null()
450 
451  select case (mudpack_case)
452  case (sll_p_separable)
453  if (present(cxx_2d) .or. present(cxy_2d) .or. present(cyy_2d) &
454  .or. present(cx_2d) .or. present(cy_2d) .or. present(ce_2d)) then
455  print *, '#2d arrays should not be here'
456  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
457  stop
458  end if
459 
460  if ((.not. (present(cxx_1d))) .and. (.not. (present(cxx)))) then
461  print *, '#1d/0d array should be here for cxx'
462  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
463  stop
464  end if
465  if (present(cxx_1d) .and. present(cxx)) then
466  print *, '#please choose between 1d or 0d array for cxx'
467  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
468  stop
469  end if
470  if (present(cxx_1d)) then
471  if (size(cxx_1d) < nc_eta1 + 1) then
472  print *, '#Bad size for cxx_1d', size(cxx_1d), nc_eta1 + 1
473  stop
474  end if
475  sll_allocate(poisson%cxx_1d(nc_eta1 + 1), ierr)
476  poisson%cxx_1d(1:nc_eta1 + 1) = cxx_1d(1:nc_eta1 + 1)
477  end if
478  if (present(cxx)) then
479  sll_allocate(poisson%cxx_1d(nc_eta1 + 1), ierr)
480  poisson%cxx_1d(1:nc_eta1 + 1) = cxx
481  end if
482 
483  if ((.not. (present(cyy_1d))) .and. (.not. (present(cyy)))) then
484  print *, '#1d/0d array should be here for cyy !'
485  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
486  stop
487  end if
488  if (present(cyy_1d) .and. present(cyy)) then
489  print *, '#please choose between 1d or 0d array for cyy'
490  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
491  stop
492  end if
493  if (present(cyy_1d)) then
494  if (size(cyy_1d) < nc_eta2 + 1) then
495  print *, '#Bad size for cyy_1d', size(cyy_1d), nc_eta2 + 1
496  stop
497  end if
498  sll_allocate(poisson%cyy_1d(nc_eta2 + 1), ierr)
499  poisson%cyy_1d(1:nc_eta2 + 1) = cyy_1d(1:nc_eta2 + 1)
500  end if
501  if (present(cyy)) then
502  sll_allocate(poisson%cyy_1d(nc_eta2 + 1), ierr)
503  poisson%cyy_1d(1:nc_eta2 + 1) = cyy
504  end if
505 
506  if ((.not. (present(cx_1d))) .and. (.not. (present(cx)))) then
507  sll_allocate(poisson%cx_1d(nc_eta1 + 1), ierr)
508  poisson%cx_1d(1:nc_eta1 + 1) = 0._f64
509  end if
510  if (present(cx_1d) .and. present(cx)) then
511  print *, '#please choose between 1d or 0d array for cx'
512  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
513  stop
514  end if
515  if (present(cx_1d)) then
516  if (size(cx_1d) < nc_eta1 + 1) then
517  print *, '#Bad size for cx_1d', size(cx_1d), nc_eta1 + 1
518  stop
519  end if
520  sll_allocate(poisson%cx_1d(nc_eta1 + 1), ierr)
521  poisson%cx_1d(1:nc_eta1 + 1) = cx_1d(1:nc_eta1 + 1)
522  end if
523  if (present(cx)) then
524  sll_allocate(poisson%cx_1d(nc_eta1 + 1), ierr)
525  poisson%cx_1d(1:nc_eta1 + 1) = cx
526  end if
527 
528  if ((.not. (present(cy_1d))) .and. (.not. (present(cy)))) then
529  sll_allocate(poisson%cy_1d(nc_eta2 + 1), ierr)
530  poisson%cy_1d(1:nc_eta2 + 1) = 0._f64
531  end if
532  if (present(cy_1d) .and. present(cy)) then
533  print *, '#please choose between 1d or 0d array for cy'
534  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
535  stop
536  end if
537  if (present(cy_1d)) then
538  if (size(cy_1d) < nc_eta2 + 1) then
539  print *, '#Bad size for cy_1d', size(cy_1d), nc_eta2 + 1
540  stop
541  end if
542  sll_allocate(poisson%cy_1d(nc_eta2 + 1), ierr)
543  poisson%cy_1d(1:nc_eta2 + 1) = cy_1d(1:nc_eta2 + 1)
544  end if
545  if (present(cy)) then
546  sll_allocate(poisson%cy_1d(nc_eta2 + 1), ierr)
547  poisson%cy_1d(1:nc_eta2 + 1) = cy
548  end if
549 
550  if ((.not. (present(cex_1d))) .and. (.not. (present(ce)))) then
551  sll_allocate(poisson%cex_1d(nc_eta1 + 1), ierr)
552  poisson%cex_1d = 0._f64
553  end if
554  if (present(cex_1d) .and. present(ce)) then
555  print *, '#please choose between 1d or 0d array for cex'
556  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
557  stop
558  end if
559  if (present(cex_1d)) then
560  if (size(cex_1d) < nc_eta1 + 1) then
561  print *, '#Bad size for cex_1d', size(cex_1d), nc_eta1 + 1
562  stop
563  end if
564  sll_allocate(poisson%cex_1d(nc_eta1 + 1), ierr)
565  poisson%cex_1d(1:nc_eta1 + 1) = cex_1d(1:nc_eta1 + 1)
566  end if
567  if (present(ce)) then
568  sll_allocate(poisson%cex_1d(nc_eta1 + 1), ierr)
569  poisson%cex_1d(1:nc_eta1 + 1) = 0.5_f64*ce
570  end if
571 
572  if ((.not. (present(cey_1d))) .and. (.not. (present(ce)))) then
573  sll_allocate(poisson%cey_1d(nc_eta2 + 1), ierr)
574  poisson%cey_1d(1:nc_eta2 + 1) = 0._f64
575  end if
576  if (present(cey_1d) .and. present(ce)) then
577  print *, '#please choose between 1d or 0d array for cey'
578  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
579  stop
580  end if
581  if (present(cey_1d)) then
582  if (size(cey_1d) < nc_eta2 + 1) then
583  print *, '#Bad size for cey_1d', size(cey_1d), nc_eta2 + 1
584  stop
585  end if
586  sll_allocate(poisson%cey_1d(nc_eta2 + 1), ierr)
587  poisson%cey_1d(1:nc_eta2 + 1) = cey_1d(1:nc_eta2 + 1)
588  end if
589  if (present(ce)) then
590  sll_allocate(poisson%cey_1d(nc_eta2 + 1), ierr)
591  poisson%cey_1d(1:nc_eta2 + 1) = 0.5_f64*ce
592  end if
593 
594  poisson%cxx_1d_interp => sll_f_new_cubic_spline_interpolator_1d( &
595  nx, &
596  eta1_min, &
597  eta1_max, &
598  sll_p_periodic)
599  call poisson%cxx_1d_interp%compute_interpolants(poisson%cxx_1d)
600 
601  poisson%cyy_1d_interp => sll_f_new_cubic_spline_interpolator_1d( &
602  ny, &
603  eta2_min, &
604  eta2_max, &
605  sll_p_periodic)
606  call poisson%cyy_1d_interp%compute_interpolants(poisson%cyy_1d)
607 
608  poisson%cx_1d_interp => sll_f_new_cubic_spline_interpolator_1d( &
609  nx, &
610  eta1_min, &
611  eta1_max, &
612  sll_p_periodic)
613  call poisson%cx_1d_interp%compute_interpolants(poisson%cx_1d)
614 
615  poisson%cy_1d_interp => sll_f_new_cubic_spline_interpolator_1d( &
616  ny, &
617  eta2_min, &
618  eta2_max, &
619  sll_p_periodic)
620  call poisson%cy_1d_interp%compute_interpolants(poisson%cy_1d)
621 
622  poisson%cex_1d_interp => sll_f_new_cubic_spline_interpolator_1d( &
623  nx, &
624  eta1_min, &
625  eta1_max, &
626  sll_p_periodic)
627  call poisson%cex_1d_interp%compute_interpolants(poisson%cex_1d)
628 
629  poisson%cey_1d_interp => sll_f_new_cubic_spline_interpolator_1d( &
630  ny, &
631  eta2_min, &
632  eta2_max, &
633  sll_p_periodic)
634  call poisson%cey_1d_interp%compute_interpolants(poisson%cey_1d)
635 
636  if (associated(mudpack_wrapper)) then
637  print *, '#Problem mudpack_wrapper is not null()'
638  stop
639  end if
640  mudpack_wrapper => poisson
641  call mud2sp(iprm, fprm, poisson%work, &
642  mudpack_cofx, &
643  mudpack_cofy, &
644  mudpack_bndsp, &
645  rhs, &
646  phi, &
647  poisson%mgopt, &
648  error)
649  mudpack_wrapper => null()
650 
651  case (sll_p_non_separable_without_cross_terms)
652 
653  if (present(cxx_1d) .or. present(cyy_1d) .or. present(cx_1d) &
654  .or. present(cy_1d) .or. present(cex_1d) .or. present(cey_1d)) then
655  print *, '#1d arrays should not be here'
656  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
657  stop
658  end if
659 
660  if ((.not. (present(cxx_2d))) .and. (.not. (present(cxx)))) then
661  print *, '#2d/0d array should be here for cxx'
662  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
663  stop
664  end if
665  if (present(cxx_2d) .and. present(cxx)) then
666  print *, '#please choose between 2d or 0d array for cxx'
667  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
668  stop
669  end if
670  if (present(cxx_2d)) then
671  if (size(cxx_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
672  print *, '#Bad size for cxx_2d', size(cxx_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
673  stop
674  end if
675  sll_allocate(poisson%cxx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
676  poisson%cxx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cxx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
677  end if
678  if (present(cxx)) then
679  sll_allocate(poisson%cxx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
680  poisson%cxx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cxx
681  end if
682 
683  if ((.not. (present(cyy_2d))) .and. (.not. (present(cyy)))) then
684  print *, '#2d/0d array should be here for cyy !'
685  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
686  stop
687  end if
688  if (present(cyy_2d) .and. present(cyy)) then
689  print *, '#please choose between 2d or 0d array for cyy'
690  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
691  stop
692  end if
693  if (present(cyy_2d)) then
694  if (size(cyy_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
695  print *, '#Bad size for cyy_2d', size(cyy_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
696  stop
697  end if
698  sll_allocate(poisson%cyy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
699  poisson%cyy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cyy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
700  end if
701  if (present(cyy)) then
702  sll_allocate(poisson%cyy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
703  poisson%cyy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cyy
704  end if
705 
706  if ((.not. (present(cx_2d))) .and. (.not. (present(cx)))) then
707  sll_allocate(poisson%cx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
708  poisson%cx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = 0._f64
709  end if
710  if (present(cx_2d) .and. present(cx)) then
711  print *, '#please choose between 2d or 0d array for cx'
712  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
713  stop
714  end if
715  if (present(cx_2d)) then
716  if (size(cx_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
717  print *, '#Bad size for cx_2d', size(cx_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
718  stop
719  end if
720  sll_allocate(poisson%cx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
721  poisson%cx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
722  end if
723  if (present(cx)) then
724  sll_allocate(poisson%cx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
725  poisson%cx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cx
726  end if
727 
728  if ((.not. (present(cy_2d))) .and. (.not. (present(cy)))) then
729  sll_allocate(poisson%cy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
730  poisson%cy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = 0._f64
731  end if
732  if (present(cy_2d) .and. present(cy)) then
733  print *, '#please choose between 2d or 0d array for cy'
734  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
735  stop
736  end if
737  if (present(cy_2d)) then
738  if (size(cy_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
739  print *, '#Bad size for cy_2d', size(cy_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
740  stop
741  end if
742  sll_allocate(poisson%cy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
743  poisson%cy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
744  end if
745  if (present(cy)) then
746  sll_allocate(poisson%cy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
747  poisson%cy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cy
748  end if
749 
750  if ((.not. (present(ce_2d))) .and. (.not. (present(ce)))) then
751  sll_allocate(poisson%ce_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
752  poisson%ce_2d = 0._f64
753  end if
754  if (present(ce_2d) .and. present(ce)) then
755  print *, '#please choose between 2d or 0d array for ce'
756  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
757  stop
758  end if
759  if (present(ce_2d)) then
760  if (size(ce_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
761  print *, '#Bad size for ce_2d', size(ce_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
762  stop
763  end if
764  sll_allocate(poisson%ce_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
765  poisson%ce_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = ce_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
766  end if
767  if (present(ce)) then
768  sll_allocate(poisson%ce_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
769  poisson%ce_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = ce
770  end if
771 
772  poisson%cxx_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
773  nx, &
774  ny, &
775  eta1_min, &
776  eta1_max, &
777  eta2_min, &
778  eta2_max, &
779  sll_p_periodic, &
780  sll_p_periodic)
781  call poisson%cxx_2d_interp%compute_interpolants(poisson%cxx_2d)
782 
783  poisson%cyy_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
784  nx, &
785  ny, &
786  eta1_min, &
787  eta1_max, &
788  eta2_min, &
789  eta2_max, &
790  sll_p_periodic, &
791  sll_p_periodic)
792  call poisson%cyy_2d_interp%compute_interpolants(poisson%cyy_2d)
793 
794  poisson%cx_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
795  nx, &
796  ny, &
797  eta1_min, &
798  eta1_max, &
799  eta2_min, &
800  eta2_max, &
801  sll_p_periodic, &
802  sll_p_periodic)
803  call poisson%cx_2d_interp%compute_interpolants(poisson%cx_2d)
804 
805  poisson%cy_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
806  nx, &
807  ny, &
808  eta1_min, &
809  eta1_max, &
810  eta2_min, &
811  eta2_max, &
812  sll_p_periodic, &
813  sll_p_periodic)
814  call poisson%cy_2d_interp%compute_interpolants(poisson%cy_2d)
815 
816  poisson%ce_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
817  nx, &
818  ny, &
819  eta1_min, &
820  eta1_max, &
821  eta2_min, &
822  eta2_max, &
823  sll_p_periodic, &
824  sll_p_periodic)
825  call poisson%ce_2d_interp%compute_interpolants(poisson%ce_2d)
826 
827  if (associated(mudpack_wrapper)) then
828  print *, '#Problem mudpack_wrapper is not null()'
829  stop
830  end if
831  mudpack_wrapper => poisson
832  call mud2(iprm, fprm, poisson%work, &
833  mudpack_cof, &
834  mudpack_bndsp, &
835  rhs, &
836  phi, &
837  poisson%mgopt, &
838  error)
839  mudpack_wrapper => null()
840 
841  case (sll_p_non_separable_with_cross_terms)
842 
843  if (present(cxx_1d) .or. present(cyy_1d) .or. present(cx_1d) &
844  .or. present(cy_1d) .or. present(cex_1d) .or. present(cey_1d)) then
845  print *, '#1d arrays should not be here'
846  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
847  stop
848  end if
849 
850  if ((.not. (present(cxx_2d))) .and. (.not. (present(cxx)))) then
851  print *, '#2d/0d array should be here for cxx'
852  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
853  stop
854  end if
855  if (present(cxx_2d) .and. present(cxx)) then
856  print *, '#please choose between 2d or 0d array for cxx'
857  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
858  stop
859  end if
860  if (present(cxx_2d)) then
861  if (size(cxx_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
862  print *, '#Bad size for cxx_2d', size(cxx_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
863  stop
864  end if
865  sll_allocate(poisson%cxx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
866  poisson%cxx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cxx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
867  end if
868  if (present(cxx)) then
869  sll_allocate(poisson%cxx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
870  poisson%cxx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cxx
871  end if
872 
873  if ((.not. (present(cxy_2d))) .and. (.not. (present(cxy)))) then
874  print *, '#2d array should be here for cxy'
875  print *, '# cxy=0. use another method'
876  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
877  stop
878  end if
879  if (present(cxy_2d) .and. present(cxy)) then
880  print *, '#please choose between 2d or 0d array for cxy'
881  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
882  stop
883  end if
884  if (present(cxy_2d)) then
885  if (size(cxy_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
886  print *, '#Bad size for cxy_2d', size(cxy_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
887  stop
888  end if
889  sll_allocate(poisson%cxy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
890  poisson%cxy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cxy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
891  end if
892  if (present(cxy)) then
893  sll_allocate(poisson%cxy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
894  poisson%cxy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cxy
895  end if
896  if ((.not. (present(cyy_2d))) .and. (.not. (present(cyy)))) then
897  print *, '#2d/0d array should be here for cyy !'
898  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
899  stop
900  end if
901  if (present(cyy_2d) .and. present(cyy)) then
902  print *, '#please choose between 2d or 0d array for cyy'
903  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
904  stop
905  end if
906  if (present(cyy_2d)) then
907  if (size(cyy_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
908  print *, '#Bad size for cyy_2d', size(cyy_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
909  stop
910  end if
911  sll_allocate(poisson%cyy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
912  poisson%cyy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cyy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
913  end if
914  if (present(cyy)) then
915  sll_allocate(poisson%cyy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
916  poisson%cyy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cyy
917  end if
918 
919  if ((.not. (present(cx_2d))) .and. (.not. (present(cx)))) then
920  sll_allocate(poisson%cx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
921  poisson%cx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = 0._f64
922  end if
923  if (present(cx_2d) .and. present(cx)) then
924  print *, '#please choose between 2d or 0d array for cx'
925  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
926  stop
927  end if
928  if (present(cx_2d)) then
929  if (size(cx_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
930  print *, '#Bad size for cx_2d', size(cx_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
931  stop
932  end if
933  sll_allocate(poisson%cx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
934  poisson%cx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
935  end if
936  if (present(cx)) then
937  sll_allocate(poisson%cx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
938  poisson%cx_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cx
939  end if
940 
941  if ((.not. (present(cy_2d))) .and. (.not. (present(cy)))) then
942  sll_allocate(poisson%cy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
943  poisson%cy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = 0._f64
944  end if
945  if (present(cy_2d) .and. present(cy)) then
946  print *, '#please choose between 2d or 0d array for cy'
947  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
948  stop
949  end if
950  if (present(cy_2d)) then
951  if (size(cy_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
952  print *, '#Bad size for cy_2d', size(cy_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
953  stop
954  end if
955  sll_allocate(poisson%cy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
956  poisson%cy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
957  end if
958  if (present(cy)) then
959  sll_allocate(poisson%cy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
960  poisson%cy_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = cy
961  end if
962 
963  if ((.not. (present(ce_2d))) .and. (.not. (present(ce)))) then
964  sll_allocate(poisson%ce_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
965  poisson%ce_2d = 0._f64
966  end if
967  if (present(ce_2d) .and. present(ce)) then
968  print *, '#please choose between 2d or 0d array for ce'
969  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
970  stop
971  end if
972  if (present(ce_2d)) then
973  if (size(ce_2d) < (nc_eta1 + 1)*(nc_eta2 + 1)) then
974  print *, '#Bad size for ce_2d', size(ce_2d), (nc_eta1 + 1)*(nc_eta2 + 1)
975  stop
976  end if
977  sll_allocate(poisson%ce_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
978  poisson%ce_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = ce_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1)
979  end if
980  if (present(ce)) then
981  sll_allocate(poisson%ce_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
982  poisson%ce_2d(1:nc_eta1 + 1, 1:nc_eta2 + 1) = ce
983  end if
984 
985  poisson%cxx_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
986  nx, &
987  ny, &
988  eta1_min, &
989  eta1_max, &
990  eta2_min, &
991  eta2_max, &
992  sll_p_periodic, &
993  sll_p_periodic)
994  call poisson%cxx_2d_interp%compute_interpolants(poisson%cxx_2d)
995 
996  poisson%cxy_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
997  nx, &
998  ny, &
999  eta1_min, &
1000  eta1_max, &
1001  eta2_min, &
1002  eta2_max, &
1003  sll_p_periodic, &
1004  sll_p_periodic)
1005  call poisson%cxy_2d_interp%compute_interpolants(poisson%cxy_2d)
1006 
1007  poisson%cyy_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
1008  nx, &
1009  ny, &
1010  eta1_min, &
1011  eta1_max, &
1012  eta2_min, &
1013  eta2_max, &
1014  sll_p_periodic, &
1015  sll_p_periodic)
1016  call poisson%cyy_2d_interp%compute_interpolants(poisson%cyy_2d)
1017 
1018  poisson%cx_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
1019  nx, &
1020  ny, &
1021  eta1_min, &
1022  eta1_max, &
1023  eta2_min, &
1024  eta2_max, &
1025  sll_p_periodic, &
1026  sll_p_periodic)
1027  call poisson%cx_2d_interp%compute_interpolants(poisson%cx_2d)
1028 
1029  poisson%cy_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
1030  nx, &
1031  ny, &
1032  eta1_min, &
1033  eta1_max, &
1034  eta2_min, &
1035  eta2_max, &
1036  sll_p_periodic, &
1037  sll_p_periodic)
1038  call poisson%cy_2d_interp%compute_interpolants(poisson%cy_2d)
1039 
1040  poisson%ce_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
1041  nx, &
1042  ny, &
1043  eta1_min, &
1044  eta1_max, &
1045  eta2_min, &
1046  eta2_max, &
1047  sll_p_periodic, &
1048  sll_p_periodic)
1049  call poisson%ce_2d_interp%compute_interpolants(poisson%ce_2d)
1050 
1051  if (associated(mudpack_wrapper)) then
1052  print *, '#Problem mudpack_wrapper is not null()'
1053  stop
1054  end if
1055  mudpack_wrapper => poisson
1056  call mud2cr(iprm, fprm, poisson%work, &
1057  mudpack_cofcr, &
1058  mudpack_bndsp, &
1059  rhs, &
1060  phi, &
1061  poisson%mgopt, &
1062  error)
1063  mudpack_wrapper => null()
1064 
1065  case default
1066  print *, '#bad mudpack_case', mudpack_case
1067  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
1068  stop
1069  end select
1070 
1071  end subroutine initialize_poisson_2d_mudpack_curvilinear
1072 
1073  ! solves \Delta phi = -rho in 2d
1074  subroutine compute_phi_from_rho_2d_mudpack(poisson, phi, rho)
1075  class(poisson_2d_mudpack_curvilinear) :: poisson
1076  sll_real64, dimension(:, :), intent(in) :: rho
1077  sll_real64, dimension(:, :), intent(out) :: phi
1078  !sll_real64 :: phi(:,:) !< Electric potential
1079  !sll_real64 :: rhs(:,:) !< Charge density
1080  !put integer and floating point argument names in contiguous
1081  !storeage for labelling in vectors iprm,fprm
1082  sll_int32 :: iprm(16)
1083  sll_real64 :: fprm(6)
1084  sll_int32 :: error
1085  sll_int32 :: i1
1086  sll_int32 :: i2
1087  sll_int32 :: intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny
1088  sll_int32 :: iguess, maxcy, method, nwork, lwrkqd, itero
1089  common/itmud2sp/intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny, &
1090  iguess, maxcy, method, nwork, lwrkqd, itero
1091  sll_real64 :: xa, xb, yc, yd, tolmax, relmax
1092  common/ftmud2sp/xa, xb, yc, yd, tolmax, relmax
1093 
1094  equivalence(intl, iprm)
1095  equivalence(xa, fprm)
1096 
1097  !set initial guess because solve should be called every time step in a
1098  !time dependent problem and the elliptic operator does not depend on time.
1099  iguess = poisson%iguess
1100 
1101  !attempt solution
1102  intl = 1
1103  !write(*,106) intl,method,iguess
1104 
1105  if (nxa == sll_p_dirichlet) then
1106  do i2 = 1, ny
1107  phi(1, i2) = 0._f64
1108  end do
1109  end if
1110  if (nxb == sll_p_dirichlet) then
1111  do i2 = 1, ny
1112  phi(nx, i2) = 0._f64
1113  end do
1114  end if
1115  if (nyc == sll_p_dirichlet) then
1116  do i1 = 1, nx
1117  phi(i1, 1) = 0._f64
1118  end do
1119  end if
1120  if (nyd == sll_p_dirichlet) then
1121  do i1 = 1, nx
1122  phi(i1, ny) = 0._f64
1123  end do
1124  end if
1125  select case (poisson%mudpack_case)
1126  case (sll_p_separable)
1127  if (associated(mudpack_wrapper)) then
1128  print *, '#Problem mudpack_wrapper is not null()'
1129  stop
1130  end if
1131  call associate_poisson(poisson)
1132 
1133  call mud2sp(iprm, &
1134  fprm, &
1135  poisson%work, &
1136  mudpack_cofx, &
1137  mudpack_cofy, &
1138  mudpack_bndsp, &
1139  -rho, &
1140  phi, &
1141  poisson%mgopt, &
1142  error)
1143  !write(*,107) error
1144  if (error > 0) stop 0
1145  ! attempt to improve approximation to fourth order
1146  call mud24sp(poisson%work, phi, error)
1147  !write (*,108) error
1148  if (error > 0) stop 0
1149 
1150  mudpack_wrapper => null()
1151 
1152  case (sll_p_non_separable_without_cross_terms)
1153  if (associated(mudpack_wrapper)) then
1154  print *, '#Problem mudpack_wrapper is not null()'
1155  stop
1156  end if
1157  call associate_poisson(poisson)
1158 
1159  call mud2(iprm, &
1160  fprm, &
1161  poisson%work, &
1162  mudpack_cof, &
1163  mudpack_bndsp, &
1164  -rho, &
1165  phi, &
1166  poisson%mgopt, &
1167  error)
1168  !write(*,107) error
1169  if (error > 0) stop 0
1170  ! attempt to improve approximation to fourth order
1171  call mud24(poisson%work, phi, error)
1172  !write (*,108) error
1173  if (error > 0) stop 0
1174 
1175  mudpack_wrapper => null()
1176 
1177  case (sll_p_non_separable_with_cross_terms)
1178  if (associated(mudpack_wrapper)) then
1179  print *, '#Problem mudpack_wrapper is not null()'
1180  stop
1181  end if
1182  call associate_poisson(poisson)
1183 
1184  call mud2cr(iprm, &
1185  fprm, &
1186  poisson%work, &
1187  mudpack_cofcr, &
1188  mudpack_bndsp, &
1189  -rho, &
1190  phi, &
1191  poisson%mgopt, &
1192  error)
1193  !write(*,107) error
1194  if (error > 0) stop 0
1195  ! attempt to improve approximation to fourth order
1196  call mud24cr(poisson%work, &
1197  mudpack_cofcr, &
1198  mudpack_bndsp, &
1199  phi, &
1200  error)
1201  !write (*,108) error
1202  if (error > 0) stop 0
1203 
1204  mudpack_wrapper => null()
1205 
1206  case default
1207  print *, '#bad mudpack_case', poisson%mudpack_case
1208  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear'
1209  stop
1210  end select
1211 
1212  end subroutine compute_phi_from_rho_2d_mudpack
1213 
1214  ! solves E = -\nabla Phi with -\Delta phi = rho in 2d
1215  subroutine compute_e_from_rho_2d_mudpack(poisson, E1, E2, rho)
1216  class(poisson_2d_mudpack_curvilinear) :: poisson
1217  sll_real64, dimension(:, :), intent(in) :: rho
1218  sll_real64, dimension(:, :), intent(out) :: e1
1219  sll_real64, dimension(:, :), intent(out) :: e2
1220 
1221  print *, '#compute_E_from_rho_2d_mudpack'
1222  print *, '#not implemented for the moment'
1223 
1224  e1 = 0._f64
1225  e2 = 0._f64
1226  print *, maxval(rho)
1227 
1228  if (.not. (associated(poisson%cxx_2d))) then
1229  print *, '#poisson%cxx_2d is not associated'
1230  end if
1231 
1232  stop
1233 
1234  !call solve( poisson%poiss, E1, E2, rho)
1235 
1236  end subroutine compute_e_from_rho_2d_mudpack
1237 
1238  function l2norm_squarred_2d_mudpack(poisson, coefs_dofs) result(r)
1239  class(poisson_2d_mudpack_curvilinear), intent(in) :: poisson
1240  sll_real64, intent(in) :: coefs_dofs(:, :)
1241  sll_real64 :: r
1242 
1243  print *, 'l2norm_squared not implemented for poisson_2d_mudpack.'
1244  r = 0.0_f64
1245 
1246  end function l2norm_squarred_2d_mudpack
1247 
1248  subroutine compute_rhs_from_function_2d_mudpack(poisson, func, coefs_dofs)
1249  class(poisson_2d_mudpack_curvilinear) :: poisson
1250  procedure(sll_i_function_of_position) :: func
1251  sll_real64, intent(out) :: coefs_dofs(:)
1252 
1253  print *, 'compute_rhs_from_function not implemented for poisson_2d_mudpack.'
1254 
1255  end subroutine compute_rhs_from_function_2d_mudpack
1256 
1257  subroutine delete_2d_mudpack(poisson)
1258  class(poisson_2d_mudpack_curvilinear) :: poisson
1259 
1260  end subroutine delete_2d_mudpack
1261 
1263  subroutine mudpack_cofx(x, cxx, cx, cex)
1264  real(8) :: x, cxx, cx, cex
1265  cxx = mudpack_wrapper%cxx_1d_interp%interpolate_from_interpolant_value(x)
1266  cx = mudpack_wrapper%cx_1d_interp%interpolate_from_interpolant_value(x)
1267  cex = mudpack_wrapper%cex_1d_interp%interpolate_from_interpolant_value(x)
1268  end subroutine mudpack_cofx
1269 
1271  subroutine mudpack_cofy(y, cyy, cy, cey)
1272  real(8) :: y, cyy, cy, cey
1273  cyy = mudpack_wrapper%cyy_1d_interp%interpolate_from_interpolant_value(y)
1274  cy = mudpack_wrapper%cy_1d_interp%interpolate_from_interpolant_value(y)
1275  cey = mudpack_wrapper%cey_1d_interp%interpolate_from_interpolant_value(y)
1276  end subroutine mudpack_cofy
1277 
1278  subroutine mudpack_cof(x, y, cxx, cyy, cx, cy, ce)
1279  real(8) :: x, cxx, cx
1280  real(8) :: y, cyy, cy, ce
1281  cxx = mudpack_wrapper%cxx_2d_interp%interpolate_from_interpolant_value(x, y)
1282  cyy = mudpack_wrapper%cyy_2d_interp%interpolate_from_interpolant_value(x, y)
1283  cx = mudpack_wrapper%cx_2d_interp%interpolate_from_interpolant_value(x, y)
1284  cy = mudpack_wrapper%cy_2d_interp%interpolate_from_interpolant_value(x, y)
1285  ce = mudpack_wrapper%ce_2d_interp%interpolate_from_interpolant_value(x, y)
1286  end subroutine mudpack_cof
1287 
1288  subroutine mudpack_cofcr(x, y, cxx, cxy, cyy, cx, cy, ce)
1289  real(8) :: x, cxx, cx, cxy
1290  real(8) :: y, cyy, cy, ce
1291  cxx = mudpack_wrapper%cxx_2d_interp%interpolate_from_interpolant_value(x, y)
1292  cxy = mudpack_wrapper%cxy_2d_interp%interpolate_from_interpolant_value(x, y)
1293  cyy = mudpack_wrapper%cyy_2d_interp%interpolate_from_interpolant_value(x, y)
1294  cx = mudpack_wrapper%cx_2d_interp%interpolate_from_interpolant_value(x, y)
1295  cy = mudpack_wrapper%cy_2d_interp%interpolate_from_interpolant_value(x, y)
1296  ce = mudpack_wrapper%ce_2d_interp%interpolate_from_interpolant_value(x, y)
1297  end subroutine mudpack_cofcr
1298 
1300  subroutine mudpack_bndsp(kbdy, xory, alfa, gbdy)
1301  integer :: kbdy
1302  real(8) :: xory, alfa, gbdy, x, y, pe, px, py
1303  real(8) :: xa, xb, yc, yd, tolmax, relmax
1304  common/ftmud2sp/xa, xb, yc, yd, tolmax, relmax
1305 
1306  pe = 0.0_8
1307  !subroutine not used in periodic case
1308  if (kbdy == 1) then ! x=xa boundary
1309  y = xory
1310  x = xa
1311  alfa = -1.0_8
1312  gbdy = px + alfa*pe
1313  return
1314  end if
1315 
1316  if (kbdy == 4) then ! y=yd boundary
1317  y = yd
1318  x = xory
1319  alfa = 1.0_8
1320  gbdy = py + alfa*pe
1321  return
1322  end if
1323  end subroutine mudpack_bndsp
1324 
1325  subroutine associate_poisson(poisson)
1326  type(poisson_2d_mudpack_curvilinear), target :: poisson
1327 
1328  mudpack_wrapper => poisson
1329 
1330  end subroutine associate_poisson
1331 
1332 end module sll_m_poisson_2d_mudpack_curvilinear
sll_m_poisson_2d_base::sll_i_function_of_position
nd real function
Definition: sll_m_poisson_2d_base.F90:46
sll_m_boundary_condition_descriptors
Describe different boundary conditions.
Definition: sll_m_boundary_condition_descriptors.F90:20
sll_m_cubic_spline_interpolator_1d
Interpolator 1d using cubic splines on regular mesh.
Definition: sll_m_cubic_spline_interpolator_1d.F90:26
sll_m_cubic_spline_interpolator_1d::sll_f_new_cubic_spline_interpolator_1d
type(sll_t_cubic_spline_interpolator_1d) function, pointer, public sll_f_new_cubic_spline_interpolator_1d(num_points, xmin, xmax, bc_type, slope_left, slope_right, fast_algorithm)
Definition: sll_m_cubic_spline_interpolator_1d.F90:260
sll_m_cubic_spline_interpolator_2d
Class for the cubic spline sll_c_interpolator_2d.
Definition: sll_m_cubic_spline_interpolator_2d.F90:23
sll_m_cubic_spline_interpolator_2d::sll_f_new_cubic_spline_interpolator_2d
type(sll_t_cubic_spline_interpolator_2d) function, pointer, public sll_f_new_cubic_spline_interpolator_2d(npts1, npts2, eta1_min, eta1_max, eta2_min, eta2_max, eta1_bc_type, eta2_bc_type, const_eta1_min_slope, const_eta1_max_slope, const_eta2_min_slope, const_eta2_max_slope, eta1_min_slopes, eta1_max_slopes, eta2_min_slopes, eta2_max_slopes)
Function that return a pointer to a cubic spline interpolator 2d object. The result can be the target...
Definition: sll_m_cubic_spline_interpolator_2d.F90:145
sll_m_interpolators_1d_base
Module for 1D interpolation and reconstruction.
Definition: sll_m_interpolators_1d_base.F90:8
sll_m_interpolators_2d_base
abstract data type for 2d interpolation
Definition: sll_m_interpolators_2d_base.F90:25
sll_m_poisson_2d_base
Module interface to solve Poisson equation in 2D.
Definition: sll_m_poisson_2d_base.F90:6
sll_m_poisson_2d_mudpack_curvilinear
Solves Poisson equation on 2d curvilinear mesh.
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:45
sll_m_poisson_2d_mudpack_curvilinear::sll_f_new_poisson_2d_mudpack_curvilinear
type(poisson_2d_mudpack_curvilinear) function, pointer, public sll_f_new_poisson_2d_mudpack_curvilinear(eta1_min, eta1_max, nc_eta1, eta2_min, eta2_max, nc_eta2, bc_eta1_left, bc_eta1_right, bc_eta2_left, bc_eta2_right, mudpack_case, cxx_2d, cxy_2d, cyy_2d, cx_2d, cy_2d, ce_2d, cxx_1d, cyy_1d, cx_1d, cy_1d, cex_1d, cey_1d, cxx, cxy, cyy, cx, cy, ce)
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:210
sll_m_poisson_2d_mudpack_curvilinear::delete_2d_mudpack
subroutine delete_2d_mudpack(poisson)
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1258
sll_m_poisson_2d_mudpack_curvilinear::mudpack_wrapper
class(poisson_2d_mudpack_curvilinear), pointer mudpack_wrapper
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:174
sll_m_poisson_2d_mudpack_curvilinear::compute_rhs_from_function_2d_mudpack
subroutine compute_rhs_from_function_2d_mudpack(poisson, func, coefs_dofs)
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1249
sll_m_poisson_2d_mudpack_curvilinear::mudpack_cof
subroutine mudpack_cof(x, y, cxx, cyy, cx, cy, ce)
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1279
sll_m_poisson_2d_mudpack_curvilinear::mudpack_cofy
subroutine mudpack_cofy(y, cyy, cy, cey)
input y dependent coefficients
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1272
sll_m_poisson_2d_mudpack_curvilinear::mudpack_cofcr
subroutine mudpack_cofcr(x, y, cxx, cxy, cyy, cx, cy, ce)
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1289
sll_m_poisson_2d_mudpack_curvilinear::mudpack_cofx
subroutine mudpack_cofx(x, cxx, cx, cex)
input x dependent coefficients
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1264
sll_m_poisson_2d_mudpack_curvilinear::associate_poisson
subroutine associate_poisson(poisson)
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1326
sll_m_poisson_2d_mudpack_curvilinear::mudpack_bndsp
subroutine mudpack_bndsp(kbdy, xory, alfa, gbdy)
input mixed derivative b.c. to mud2sp
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1301
sll_m_poisson_2d_mudpack_curvilinear::compute_e_from_rho_2d_mudpack
subroutine compute_e_from_rho_2d_mudpack(poisson, E1, E2, rho)
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1216
sll_m_poisson_2d_mudpack_curvilinear::initialize_poisson_2d_mudpack_curvilinear
subroutine initialize_poisson_2d_mudpack_curvilinear(poisson, eta1_min, eta1_max, nc_eta1, eta2_min, eta2_max, nc_eta2, bc_eta1_left, bc_eta1_right, bc_eta2_left, bc_eta2_right, mudpack_case, cxx_2d, cxy_2d, cyy_2d, cx_2d, cy_2d, ce_2d, cxx_1d, cyy_1d, cx_1d, cy_1d, cex_1d, cey_1d, cxx, cxy, cyy, cx, cy, ce)
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:310
sll_m_poisson_2d_mudpack_curvilinear::l2norm_squarred_2d_mudpack
real(kind=f64) function l2norm_squarred_2d_mudpack(poisson, coefs_dofs)
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1239
sll_m_poisson_2d_mudpack_curvilinear::compute_phi_from_rho_2d_mudpack
subroutine compute_phi_from_rho_2d_mudpack(poisson, phi, rho)
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:1075
sll_m_interpolators_1d_base::sll_c_interpolator_1d
Abstract class for 1D interpolation and reconstruction.
Definition: sll_m_interpolators_1d_base.F90:21
sll_m_interpolators_2d_base::sll_c_interpolator_2d
Base class/basic interface for 2D interpolators.
Definition: sll_m_interpolators_2d_base.F90:44
sll_m_poisson_2d_base::sll_c_poisson_2d_base
Definition: sll_m_poisson_2d_base.F90:20
sll_m_poisson_2d_mudpack_curvilinear::poisson_2d_mudpack_curvilinear
Derived type to solve Poisson equation on 2d curvilinear mesh.
Definition: sll_m_poisson_2d_mudpack_curvilinear.F90:84
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1