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_old.F90
Go to the documentation of this file.
1 #ifndef DOXYGEN_SHOULD_SKIP_THIS
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 
23 !solves \sum_{i,j=1}^2 A_{i,j}\partial_{i,j} phi
24 ! +\sum_{i=1}^2B_i\partial_i phi
25 ! +C \phi = rho
26 !in polar coordinates
27 !self leads when A_{1,2}=A_{2,1}=0 and B_2 = 0
28 ! A_11\partial_{1,1}\hat{phi}+B_1\partial_{1}\hat{phi}+(C+A_{2,2}k^2)\hat{phi} = \hat{rho}
29 
31 module sll_m_poisson_2d_mudpack_curvilinear_old
32 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33 #include "sll_memory.h"
34 #include "sll_working_precision.h"
35 
36 ! use F77_mudpack, only: &
37 ! mud2, &
38 ! mud24, &
39 ! mud24cr, &
40 ! mud2cr
41 
42  use sll_m_boundary_condition_descriptors, only: &
43  sll_p_dirichlet
44 
45  use sll_m_cartesian_meshes, only: &
46  sll_t_cartesian_mesh_2d
47 
48  use sll_m_coordinate_transformation_2d_base, only: &
49  sll_c_coordinate_transformation_2d_base
50 
51  use sll_m_cubic_spline_interpolator_2d, only: &
52  sll_f_new_cubic_spline_interpolator_2d
53 
54  use sll_m_interpolators_2d_base, only: &
55  sll_c_interpolator_2d
56 
57  use sll_m_mudpack_curvilinear, only: &
58  sll_p_non_separable_with_cross_terms, &
59  sll_p_non_separable_without_cross_terms
60 
61  use sll_m_poisson_2d_base, only: &
62  sll_c_poisson_2d_base, &
63  sll_i_function_of_position
64 
65  implicit none
66 
67  public :: &
68  sll_f_new_poisson_2d_mudpack_curvilinear_old
69 
70  private
71 
72 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
73 
74  type, extends(sll_c_poisson_2d_base) :: poisson_2d_mudpack_curvilinear_old
75 
76  !type(sll_plan_poisson_polar), pointer :: poiss
77  sll_real64, dimension(:, :), pointer :: cxx_2d
78  sll_real64, dimension(:, :), pointer :: cxy_2d
79  sll_real64, dimension(:, :), pointer :: cyy_2d
80  sll_real64, dimension(:, :), pointer :: cx_2d
81  sll_real64, dimension(:, :), pointer :: cy_2d
82  sll_real64, dimension(:, :), pointer :: ce_2d
83  sll_real64, dimension(:, :), pointer :: rho
84  class(sll_c_interpolator_2d), pointer :: cxx_2d_interp
85  class(sll_c_interpolator_2d), pointer :: cyy_2d_interp
86  class(sll_c_interpolator_2d), pointer :: cxy_2d_interp
87  class(sll_c_interpolator_2d), pointer :: cx_2d_interp
88  class(sll_c_interpolator_2d), pointer :: cy_2d_interp
89  class(sll_c_interpolator_2d), pointer :: ce_2d_interp
90  class(sll_c_interpolator_2d), pointer :: a11_interp
91  class(sll_c_interpolator_2d), pointer :: a22_interp
92  class(sll_c_interpolator_2d), pointer :: a12_interp
93  class(sll_c_interpolator_2d), pointer :: a21_interp
94  class(sll_c_coordinate_transformation_2d_base), pointer :: transformation
95  sll_int32 :: mudpack_curvilinear_case
96  sll_real64, dimension(:), pointer :: work
97  sll_int32 :: mgopt(4)
98  sll_int32 :: iprm(16)
99  sll_real64 :: fprm(6)
100  sll_int32 :: iguess
101 
102  contains
103  procedure, pass(poisson) :: initialize => &
104  initialize_poisson_2d_mudpack_curvilinear_old
105  procedure, pass(poisson) :: compute_phi_from_rho => &
106  compute_phi_from_rho_2d_mudpack_curvilinear
107  procedure, pass(poisson) :: compute_E_from_rho => &
108  compute_e_from_rho_2d_mudpack_curvilinear
109 ! procedure, pass(poisson) :: compute_E_from_phi => &
110 ! compute_E_from_phi_2d_polar
111 
113  procedure :: &
114  l2norm_squared => l2norm_squarred_2d_mudpack_curvilinear
116  procedure :: &
117  compute_rhs_from_function => compute_rhs_from_function_2d_mudpack_curvilinear
119  procedure :: free => delete_2d_mudpack_curvilinear_solver
120 
121  end type poisson_2d_mudpack_curvilinear_old
122 
123  class(poisson_2d_mudpack_curvilinear_old), pointer :: mudpack_curvilinear_wrapper => null()
124 
125 contains
126  function sll_f_new_poisson_2d_mudpack_curvilinear_old( &
127  transf, &
128  eta1_min, &
129  eta1_max, &
130  nc_eta1, &
131  eta2_min, &
132  eta2_max, &
133  nc_eta2, &
134  bc_eta1_left, &
135  bc_eta1_right, &
136  bc_eta2_left, &
137  bc_eta2_right, &
138  bc_interp2d_eta1, &
139  bc_interp2d_eta2, &
140  b11, &
141  b12, &
142  b21, &
143  b22, &
144  c, &
145  mudpack_curvilinear_case) &
146  result(poisson)
147 
148  type(poisson_2d_mudpack_curvilinear_old), pointer :: poisson
149  sll_real64, intent(in) :: eta1_min
150  sll_real64, intent(in) :: eta1_max
151  sll_int32, intent(in) :: nc_eta1
152  sll_real64, intent(in) :: eta2_min
153  sll_real64, intent(in) :: eta2_max
154  sll_int32, intent(in) :: nc_eta2
155  sll_int32, intent(in) :: bc_eta1_left
156  sll_int32, intent(in) :: bc_eta1_right
157  sll_int32, intent(in) :: bc_eta2_left
158  sll_int32, intent(in) :: bc_eta2_right
159  sll_int32, intent(in) :: bc_interp2d_eta1
160  sll_int32, intent(in) :: bc_interp2d_eta2
161  sll_int32, intent(in), optional :: mudpack_curvilinear_case
162  sll_real64, dimension(:, :), intent(in) :: b11
163  sll_real64, dimension(:, :), intent(in) :: b12
164  sll_real64, dimension(:, :), intent(in) :: b21
165  sll_real64, dimension(:, :), intent(in) :: b22
166  sll_real64, dimension(:, :), intent(in) :: c
167  class(sll_c_coordinate_transformation_2d_base), pointer, intent(in) :: transf
168 
169  sll_int32 :: ierr
170 
171  sll_allocate(poisson, ierr)
172 
173  call initialize_poisson_2d_mudpack_curvilinear_old( &
174  poisson, &
175  transf, &
176  eta1_min, &
177  eta1_max, &
178  nc_eta1, &
179  eta2_min, &
180  eta2_max, &
181  nc_eta2, &
182  bc_eta1_left, &
183  bc_eta1_right, &
184  bc_eta2_left, &
185  bc_eta2_right, &
186  bc_interp2d_eta1, &
187  bc_interp2d_eta2, &
188  b11, &
189  b12, &
190  b21, &
191  b22, &
192  c, &
193  mudpack_curvilinear_case)
194 
195  end function sll_f_new_poisson_2d_mudpack_curvilinear_old
196 
197  subroutine initialize_poisson_2d_mudpack_curvilinear_old( &
198  poisson, &
199  transf, &
200  eta1_min, &
201  eta1_max, &
202  nc_eta1, &
203  eta2_min, &
204  eta2_max, &
205  nc_eta2, &
206  bc_eta1_left, &
207  bc_eta1_right, &
208  bc_eta2_left, &
209  bc_eta2_right, &
210  bc_interp2d_eta1, &
211  bc_interp2d_eta2, &
212  b11, &
213  b12, &
214  b21, &
215  b22, &
216  c, &
217  mudpack_curvilinear_case)
218  class(poisson_2d_mudpack_curvilinear_old), target :: poisson
219  sll_real64, intent(in) :: eta1_min
220  sll_real64, intent(in) :: eta1_max
221  sll_int32, intent(in) :: nc_eta1
222  sll_real64, intent(in) :: eta2_min
223  sll_real64, intent(in) :: eta2_max
224  sll_int32, intent(in) :: nc_eta2
225  sll_int32, intent(in) :: bc_eta1_left
226  sll_int32, intent(in) :: bc_eta1_right
227  sll_int32, intent(in) :: bc_eta2_left
228  sll_int32, intent(in) :: bc_eta2_right
229  sll_int32, intent(in) :: bc_interp2d_eta1
230  sll_int32, intent(in) :: bc_interp2d_eta2
231  sll_int32, intent(in), optional :: mudpack_curvilinear_case
232  sll_real64, dimension(:, :), intent(in) :: b11
233  sll_real64, dimension(:, :), intent(in) :: b12
234  sll_real64, dimension(:, :), intent(in) :: b21
235  sll_real64, dimension(:, :), intent(in) :: b22
236  sll_real64, dimension(:, :), intent(in) :: c
237  class(sll_c_coordinate_transformation_2d_base), pointer :: transf
238  sll_real64, dimension(:, :), allocatable :: a12_array
239  sll_real64, dimension(:, :), allocatable :: a21_array
240  sll_int32 :: ierr
241  !!!! begin variables for mudpack_curvilinear
242  sll_int32, parameter :: iixp = 2, jjyq = 2
243  sll_int32 :: icall, iiex, jjey, llwork
244  sll_real64, pointer :: phi(:)
245  sll_real64, pointer :: rhs(:)
246  !put integer and floating point argument names in contiguous
247  !storeage for labelling in vectors iprm,fprm
248  sll_int32 :: iprm(16)
249  sll_real64 :: fprm(6)
250  !sll_int32 :: i
251  sll_int32 :: error
252  sll_int32 :: intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny
253  sll_int32 :: iguess, maxcy, method, nwork, lwrkqd, itero
254  common/itmud2sp/intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny, &
255  iguess, maxcy, method, nwork, lwrkqd, itero
256  sll_real64 :: xa, xb, yc, yd, tolmax, relmax
257  common/ftmud2sp/xa, xb, yc, yd, tolmax, relmax
258  equivalence(intl, iprm)
259  equivalence(xa, fprm)
260 
261  !!!! end variables for mudpack_curvilinear
262  sll_real64 :: delta1, delta2
263 
264  nx = nc_eta1 + 1
265  ny = nc_eta2 + 1
266 
267  allocate (phi(nx*ny))
268  allocate (rhs(nx*ny))
269 
270  delta1 = (eta1_max - eta1_min)/real(nc_eta1, f64)
271  delta2 = (eta2_max - eta2_min)/real(nc_eta2, f64)
272  ! set minimum required work space
273  llwork = (7*(nx + 2)*(ny + 2) + 44*nx*ny)/3
274 
275  allocate (poisson%work(llwork))
276  icall = 0
277  iiex = ceiling(log((nx - 1.)/iixp)/log(2.)) + 1
278  jjey = ceiling(log((ny - 1.)/jjyq)/log(2.)) + 1
279 
280  !set input integer arguments
281  intl = 0
282 
283  !set boundary condition flags
284  nxa = bc_eta1_left
285  nxb = bc_eta1_right
286  nyc = bc_eta2_left
287  nyd = bc_eta2_right
288 
289  !set grid sizes from parameter statements
290  ixp = iixp
291  jyq = jjyq
292  iex = iiex
293  jey = jjey
294 
295  nx = ixp*(2**(iex - 1)) + 1
296  ny = jyq*(2**(jey - 1)) + 1
297 
298  if (nx /= nc_eta1 + 1 .or. ny /= nc_eta2 + 1) then
299  print *, "nx,nc_eta1+1=", nx, nc_eta1 + 1
300  print *, "ny,nc_eta2+1=", ny, nc_eta2 + 1
301  stop ' nx or ny different in sll_mudpack_curvilinear_cartesian '
302  end if
303 
304  !set multigrid arguments (w(2,1) cycling with fully weighted
305  !residual restriction and cubic prolongation)
306  poisson%mgopt(1) = 2
307  poisson%mgopt(2) = 2
308  poisson%mgopt(3) = 1
309  poisson%mgopt(4) = 3
310 
311  !set for three cycles to ensure second-order approximation is computed
312  maxcy = 3
313 
314  !set no initial guess forcing full multigrid cycling
315  poisson%iguess = 0
316  iguess = poisson%iguess
317 
318  !set work space length approximation from parameter statement
319  nwork = llwork
320 
321  !set point relaxation
322  method = 0
323 
324  !set end points of solution rectangle in (x,y) space
325  xa = eta1_min
326  xb = eta1_max
327  yc = eta2_min
328  yd = eta2_max
329 
330  !set for no error control flag
331  tolmax = 0.0_f64
332 
333 ! write(*,101) (iprm(i),i=1,15)
334 ! write(*,102) (poisson%mgopt(i),i=1,4)
335 ! write(*,103) xa,xb,yc,yd,tolmax
336 ! write(*,104) intl
337 
338 !call mud2sp(iprm,fprm,self%work,cofx,cofy,bndsp,rhs,phi,self%mgopt,error)
339  if (present(mudpack_curvilinear_case)) then
340  poisson%mudpack_curvilinear_case = mudpack_curvilinear_case
341  else
342  poisson%mudpack_curvilinear_case = sll_p_non_separable_with_cross_terms
343  end if
344  poisson%transformation => transf
345  poisson%cxx_2d_interp => null()
346  poisson%cyy_2d_interp => null()
347  poisson%cx_2d_interp => null()
348  poisson%cy_2d_interp => null()
349  poisson%ce_2d_interp => null()
350  poisson%a12_interp => null()
351  poisson%a21_interp => null()
352  sll_allocate(poisson%rho(nc_eta1 + 1, nc_eta2 + 1), ierr)
353 
354  sll_allocate(poisson%cxx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
355  !SLL_ALLOCATE(poisson%cxy_2d(nc_eta1+1,nc_eta2+1),ierr)
356  sll_allocate(poisson%cyy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
357  sll_allocate(poisson%cx_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
358  sll_allocate(poisson%cy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
359  sll_allocate(poisson%ce_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
360 
361  poisson%a12_interp => sll_f_new_cubic_spline_interpolator_2d( &
362  nx, &
363  ny, &
364  eta1_min, &
365  eta1_max, &
366  eta2_min, &
367  eta2_max, &
368  bc_interp2d_eta1, &
369  bc_interp2d_eta2)
370  poisson%a21_interp => sll_f_new_cubic_spline_interpolator_2d( &
371  nx, &
372  ny, &
373  eta1_min, &
374  eta1_max, &
375  eta2_min, &
376  eta2_max, &
377  bc_interp2d_eta1, &
378  bc_interp2d_eta2)
379  sll_allocate(a12_array(nx, ny), ierr)
380  sll_allocate(a21_array(nx, ny), ierr)
381  call a12_a21_array(b11, b12, b21, b22, transf, eta1_min, &
382  eta2_min, delta1, delta2, nx, ny, a12_array, a21_array)
383  call poisson%a12_interp%compute_interpolants(a12_array)
384  call poisson%a21_interp%compute_interpolants(a21_array)
385 
386  poisson%cxx_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
387  nx, &
388  ny, &
389  eta1_min, &
390  eta1_max, &
391  eta2_min, &
392  eta2_max, &
393  bc_interp2d_eta1, &
394  bc_interp2d_eta2)
395  poisson%cyy_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
396  nx, &
397  ny, &
398  eta1_min, &
399  eta1_max, &
400  eta2_min, &
401  eta2_max, &
402  bc_interp2d_eta1, &
403  bc_interp2d_eta2)
404  call coefxxyy_array(b11, b12, b21, b22, transf, eta1_min, eta2_min, &
405  delta1, delta2, nx, ny, poisson%cxx_2d, poisson%cyy_2d)
406  call poisson%cxx_2d_interp%compute_interpolants(poisson%cxx_2d)
407  call poisson%cyy_2d_interp%compute_interpolants(poisson%cyy_2d)
408 
409  poisson%cx_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
410  nx, &
411  ny, &
412  eta1_min, &
413  eta1_max, &
414  eta2_min, &
415  eta2_max, &
416  bc_interp2d_eta1, &
417  bc_interp2d_eta2)
418  call coefx_array(eta1_min, eta2_min, delta1, delta2, nx, ny, &
419  poisson%cxx_2d_interp, poisson%a21_interp, poisson%cx_2d)
420  call poisson%cx_2d_interp%compute_interpolants(poisson%cx_2d)
421 
422  poisson%cy_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
423  nx, &
424  ny, &
425  eta1_min, &
426  eta1_max, &
427  eta2_min, &
428  eta2_max, &
429  bc_interp2d_eta1, &
430  bc_interp2d_eta2)
431  call coefy_array(eta1_min, eta2_min, delta1, delta2, nx, ny, &
432  poisson%cyy_2d_interp, poisson%a12_interp, poisson%cy_2d)
433  call poisson%cy_2d_interp%compute_interpolants(poisson%cy_2d)
434 
435  poisson%ce_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
436  nx, &
437  ny, &
438  eta1_min, &
439  eta1_max, &
440  eta2_min, &
441  eta2_max, &
442  bc_interp2d_eta1, &
443  bc_interp2d_eta2)
444  poisson%ce_2d = -c
445  call poisson%ce_2d_interp%compute_interpolants(poisson%ce_2d)
446 
447  !******sll_p_non_separable_with_cross_terms)
448  select case (poisson%mudpack_curvilinear_case)
449 
450  case (sll_p_non_separable_without_cross_terms)
451  if (associated(mudpack_curvilinear_wrapper)) then
452  print *, '#Problem mudpack_curvilinear_wrapper is not null()'
453  stop
454  end if
455  call associate_poisson(poisson)
456  call mud2(iprm, fprm, poisson%work, &
457  mudpack_curvilinear_cof, &
458  mudpack_curvilinear_bndcr, &
459  rhs, &
460  phi, &
461  poisson%mgopt, &
462  error)
463  mudpack_curvilinear_wrapper => null()
464 
465  case (sll_p_non_separable_with_cross_terms)
466  sll_allocate(poisson%cxy_2d(nc_eta1 + 1, nc_eta2 + 1), ierr)
467  poisson%cxy_2d_interp => sll_f_new_cubic_spline_interpolator_2d( &
468  nx, &
469  ny, &
470  eta1_min, &
471  eta1_max, &
472  eta2_min, &
473  eta2_max, &
474  bc_interp2d_eta1, &
475  bc_interp2d_eta2)
476  call coefxy_array(b11, b12, b21, b22, transf, eta1_min, eta2_min, &
477  delta1, delta2, nx, ny, poisson%cxy_2d)
478  call poisson%cxy_2d_interp%compute_interpolants(poisson%cxy_2d)
479  if (associated(mudpack_curvilinear_wrapper)) then
480  print *, '#Problem mudpack_curvilinear_wrapper is not null()'
481  stop
482  end if
483  call associate_poisson(poisson)
484  call mud2cr(iprm, fprm, poisson%work, &
485  mudpack_curvilinear_cofcr, &
486  mudpack_curvilinear_bndcr, &
487  rhs, &
488  phi, &
489  poisson%mgopt, &
490  error)
491  mudpack_curvilinear_wrapper => null()
492  case default
493  print *, '#bad mudpack_curvilinear_case', poisson%mudpack_curvilinear_case
494  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear_old'
495  stop
496  end select
497 
498  end subroutine initialize_poisson_2d_mudpack_curvilinear_old
499 
500  ! solves -\Delta phi = rho in 2d
501  subroutine compute_phi_from_rho_2d_mudpack_curvilinear(poisson, phi, rho)
502  class(poisson_2d_mudpack_curvilinear_old) :: poisson
503  sll_real64, dimension(:, :), intent(in) :: rho
504  sll_real64, dimension(:, :), intent(out) :: phi
505  sll_int32 :: nc_eta1
506  sll_int32 :: nc_eta2
507  sll_real64 :: eta1_min
508  sll_real64 :: eta2_min
509  sll_real64 :: delta_eta1
510  sll_real64 :: delta_eta2
511  sll_real64 :: eta1
512  sll_real64 :: eta2
513  sll_int32 :: i1
514  sll_int32 :: i2
515  !sll_real64 :: phi(:,:) !< Electric potential
516  !sll_real64 :: rhs(:,:) !< Charge density
517  !put integer and floating point argument names in contiguous
518  !storeage for labelling in vectors iprm,fprm
519  sll_int32 :: iprm(16)
520  sll_real64 :: fprm(6)
521  sll_int32 :: error
522  sll_int32 :: intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny
523  sll_int32 :: iguess, maxcy, method, nwork, lwrkqd, itero
524 
525  common/itmud2sp/intl, nxa, nxb, nyc, nyd, ixp, jyq, iex, jey, nx, ny, &
526  iguess, maxcy, method, nwork, lwrkqd, itero
527  sll_real64 :: xa, xb, yc, yd, tolmax, relmax
528  common/ftmud2sp/xa, xb, yc, yd, tolmax, relmax
529 
530  equivalence(intl, iprm)
531  equivalence(xa, fprm)
532 
533  !set initial guess because solve should be called every time step in a
534  !time dependent problem and the elliptic operator does not depend on time.
535  iguess = poisson%iguess
536 
537  !attempt solution
538  intl = 1
539  !write(*,106) intl,method,iguess
540 
541  associate(mesh => poisson%transformation%mesh)
542 
543  nc_eta1 = mesh%num_cells1
544  nc_eta2 = mesh%num_cells2
545  eta1_min = mesh%eta1_min
546  eta2_min = mesh%eta2_min
547  delta_eta1 = mesh%delta_eta1
548  delta_eta2 = mesh%delta_eta2
549 
550  end associate
551 
552  poisson%rho(1:nc_eta1 + 1, 1:nc_eta2 + 1) = 0._f64
553  do i2 = 1, nc_eta2 + 1
554  eta2 = eta2_min + real(i2 - 1, f64)*delta_eta2
555  do i1 = 1, nc_eta1 + 1
556  eta1 = eta1_min + real(i1 - 1, f64)*delta_eta1
557  poisson%rho(i1, i2) = -rho(i1, i2)*poisson%transformation%jacobian(eta1, eta2)
558  end do
559  end do
560 
561  if (nxa == sll_p_dirichlet) then
562  do i2 = 1, nc_eta2 + 1
563  phi(1, i2) = 0._f64
564  end do
565  end if
566  if (nxb == sll_p_dirichlet) then
567  do i2 = 1, nc_eta2 + 1
568  phi(nc_eta1 + 1, i2) = 0._f64
569  end do
570  end if
571  if (nyc == sll_p_dirichlet) then
572  do i1 = 1, nc_eta1 + 1
573  phi(i1, 1) = 0._f64
574  end do
575  end if
576  if (nyd == sll_p_dirichlet) then
577  do i1 = 1, nc_eta1 + 1
578  phi(i1, nc_eta2 + 1) = 0._f64
579  end do
580  end if
581 
582  select case (poisson%mudpack_curvilinear_case)
583 
584  case (sll_p_non_separable_without_cross_terms)
585  if (associated(mudpack_curvilinear_wrapper)) then
586  print *, '#Problem mudpack_curvilinear_wrapper is not null()'
587  stop
588  end if
589  call associate_poisson(poisson)
590  call mud2(iprm, &
591  fprm, &
592  poisson%work, &
593  mudpack_curvilinear_cof, &
594  mudpack_curvilinear_bndcr, &
595  poisson%rho, &
596  phi, &
597  poisson%mgopt, &
598  error)
599  !write(*,107) error
600  if (error > 0) stop 0
601  ! attempt to improve approximation to fourth order
602  ! seems not to work for the moment
603  call mud24(poisson%work, phi, error)
604  !write (*,108) error
605  if (error > 0) stop 0
606  mudpack_curvilinear_wrapper => null()
607 
608  case (sll_p_non_separable_with_cross_terms)
609  if (associated(mudpack_curvilinear_wrapper)) then
610  print *, '#Problem mudpack_curvilinear_wrapper is not null()'
611  stop
612  end if
613  call associate_poisson(poisson)
614  call mud2cr(iprm, &
615  fprm, &
616  poisson%work, &
617  mudpack_curvilinear_cofcr, &
618  mudpack_curvilinear_bndcr, &
619  poisson%rho, &
620  phi, &
621  poisson%mgopt, &
622  error)
623  !write(*,107) error
624  if (error > 0) stop 0
625  ! attempt to improve approximation to fourth order
626  ! seems not to work for the moment
627  call mud24cr(poisson%work, &
628  mudpack_curvilinear_cofcr, &
629  mudpack_curvilinear_bndcr, &
630  phi, &
631  error)
632  !write (*,108) error
633  if (error > 0) stop 0
634  mudpack_curvilinear_wrapper => null()
635 
636  case default
637  print *, '#bad mudpack_curvilinear_case', poisson%mudpack_curvilinear_case
638  print *, '#in subroutine initialize_poisson_2d_mudpack_curvilinear_old'
639  stop
640  end select
641 
642  end subroutine compute_phi_from_rho_2d_mudpack_curvilinear
643 
644  ! solves E = -\nabla Phi with -\Delta phi = rho in 2d
645  subroutine compute_e_from_rho_2d_mudpack_curvilinear(poisson, E1, E2, rho)
646  class(poisson_2d_mudpack_curvilinear_old) :: poisson
647  sll_real64, dimension(:, :), intent(in) :: rho
648  sll_real64, dimension(:, :), intent(out) :: e1
649  sll_real64, dimension(:, :), intent(out) :: e2
650 
651  print *, '#compute_E_from_rho_2d_mudpack_curvilinear'
652  print *, '#not implemented for the moment'
653  e1 = 0._f64
654  e2 = 0._f64
655  print *, maxval(rho)
656 
657  if (.not. (associated(poisson%cxx_2d))) then
658  print *, '#poisson%cxx_2d is not associated'
659  end if
660 
661  stop
662 
663  !call solve( poisson%poiss, E1, E2, rho)
664 
665  end subroutine compute_e_from_rho_2d_mudpack_curvilinear
666 
667  function l2norm_squarred_2d_mudpack_curvilinear(poisson, coefs_dofs) result(r)
668  class(poisson_2d_mudpack_curvilinear_old), intent(in) :: poisson
669  sll_real64, intent(in) :: coefs_dofs(:, :)
670  sll_real64 :: r
671 
672  print *, 'l2norm_squared not implemented for poisson_2d_mudpack_curvilinear_solver.'
673  r = 0.0_f64
674 
675  end function l2norm_squarred_2d_mudpack_curvilinear
676 
677  subroutine compute_rhs_from_function_2d_mudpack_curvilinear(poisson, func, coefs_dofs)
678  class(poisson_2d_mudpack_curvilinear_old) :: poisson
679  procedure(sll_i_function_of_position) :: func
680  sll_real64, intent(out) :: coefs_dofs(:)
681 
682  print *, 'compute_rhs_from_function not implemented for poisson_2d_mudpack_curvilinear_solver.'
683 
684  end subroutine compute_rhs_from_function_2d_mudpack_curvilinear
685 
686  subroutine delete_2d_mudpack_curvilinear_solver(poisson)
687  class(poisson_2d_mudpack_curvilinear_old) :: poisson
688 
689  end subroutine delete_2d_mudpack_curvilinear_solver
690 
691  subroutine coefxxyy_array(b11, b12, b21, b22, transf, eta1_min, eta2_min, &
692  delta1, delta2, nx, ny, cxx_array, cyy_array)
693  implicit none
694  sll_real64 :: eta1, eta1_min, eta2_min
695  sll_real64 :: eta2, delta1, delta2
696  sll_int32 :: i, j, nx, ny
697  sll_real64, dimension(:, :):: cxx_array, cyy_array
698  sll_real64, dimension(1:2, 1:2) :: jac_m
699  class(sll_c_coordinate_transformation_2d_base), pointer :: transf
700  sll_real64, dimension(:, :) :: b11
701  sll_real64, dimension(:, :) :: b12
702  sll_real64, dimension(:, :) :: b21
703  sll_real64, dimension(:, :) :: b22
704 
705  do j = 1, ny
706  eta2 = eta2_min + real(j - 1, f64)*delta2
707  do i = 1, nx
708  eta1 = eta1_min + real(i - 1, f64)*delta1
709  jac_m = transf%jacobian_matrix(eta1, eta2)
710  cxx_array(i, j) = (b11(i, j)*(jac_m(1, 2)*jac_m(1, 2) + jac_m(2, 2)*jac_m(2, 2)) - &
711  & b12(i, j)*(jac_m(2, 1)*jac_m(2, 2) + jac_m(1, 1)*jac_m(1, 2))) &
712  /transf%jacobian(eta1, eta2)
713  cyy_array(i, j) = (b22(i, j)*(jac_m(2, 1)*jac_m(2, 1) + jac_m(1, 1)*jac_m(1, 1)) - &
714  & b21(i, j)*(jac_m(2, 1)*jac_m(2, 2) + jac_m(1, 1)*jac_m(1, 2))) &
715  /transf%jacobian(eta1, eta2)
716  end do
717  end do
718  end subroutine coefxxyy_array
719 
720  subroutine coefxy_array(b11, b12, b21, b22, transf, eta1_min, eta2_min, &
721  delta1, delta2, nx, ny, cxy_array)
722  implicit none
723  sll_real64 :: eta1, eta1_min, eta2_min
724  sll_real64 :: eta2, delta1, delta2
725  sll_real64 :: a12, a21
726  sll_int32 :: i, j, nx, ny
727  sll_real64, dimension(:, :):: cxy_array
728  sll_real64, dimension(1:2, 1:2) :: jac_m
729  class(sll_c_coordinate_transformation_2d_base), pointer :: transf
730  sll_real64, dimension(:, :) :: b11
731  sll_real64, dimension(:, :) :: b12
732  sll_real64, dimension(:, :) :: b21
733  sll_real64, dimension(:, :) :: b22
734 
735  do j = 1, ny
736  eta2 = eta2_min + real(j - 1, f64)*delta2
737  do i = 1, nx
738  eta1 = eta1_min + real(i - 1, f64)*delta1
739  jac_m = transf%jacobian_matrix(eta1, eta2)
740  a12 = b12(i, j)*(jac_m(2, 1)*jac_m(2, 1) + jac_m(1, 1)*jac_m(1, 1)) - &
741  & b11(i, j)*(jac_m(2, 1)*jac_m(2, 2) + jac_m(1, 1)*jac_m(1, 2))
742 
743  a21 = b21(i, j)*(jac_m(1, 2)*jac_m(1, 2) + jac_m(2, 2)*jac_m(2, 2)) - &
744  & b22(i, j)*(jac_m(2, 1)*jac_m(2, 2) + jac_m(1, 1)*jac_m(1, 2))
745  cxy_array(i, j) = (a12 + a21)/transf%jacobian(eta1, eta2)
746  !write(100,*) eta1,eta2, cxy_array(i,j), transf%jacobian(eta1,eta2)
747  end do
748  end do
749  end subroutine coefxy_array
750 
751  subroutine a12_a21_array(b11, b12, b21, b22, transf, eta1_min, eta2_min, delta1, delta2, nx, ny, a12_array, a21_array)
752  implicit none
753  sll_real64 :: eta1, eta1_min, eta2_min
754  sll_real64 :: eta2, delta1, delta2
755  sll_real64 :: a12, a21
756  sll_int32 :: i, j, nx, ny
757  sll_real64, dimension(:, :):: a12_array
758  sll_real64, dimension(:, :):: a21_array
759  class(sll_c_coordinate_transformation_2d_base), pointer :: transf
760  sll_real64, dimension(:, :) :: b11
761  sll_real64, dimension(:, :) :: b12
762  sll_real64, dimension(:, :) :: b21
763  sll_real64, dimension(:, :) :: b22
764  sll_real64, dimension(1:2, 1:2) :: jac_m
765 
766  do j = 1, ny
767  eta2 = eta2_min + real(j - 1, f64)*delta2
768  do i = 1, nx
769  eta1 = eta1_min + real(i - 1, f64)*delta1
770  jac_m = transf%jacobian_matrix(eta1, eta2)
771  a12 = b12(i, j)*(jac_m(2, 1)*jac_m(2, 1) + jac_m(1, 1)*jac_m(1, 1)) - &
772  & b11(i, j)*(jac_m(2, 1)*jac_m(2, 2) + jac_m(1, 1)*jac_m(1, 2))
773  a12_array(i, j) = a12/transf%jacobian(eta1, eta2)
774  a21 = b21(i, j)*(jac_m(1, 2)*jac_m(1, 2) + jac_m(2, 2)*jac_m(2, 2)) - &
775  & b22(i, j)*(jac_m(2, 1)*jac_m(2, 2) + jac_m(1, 1)*jac_m(1, 2))
776  a21_array(i, j) = a21/transf%jacobian(eta1, eta2)
777  end do
778  end do
779  end subroutine a12_a21_array
780 
781  subroutine coefx_array(eta1_min, eta2_min, &
782  delta1, delta2, nx, ny, cxx_2d_interp, a21_interp, cx_array)
783  implicit none
784  sll_real64 :: eta1, eta1_min, eta2_min
785  sll_real64 :: eta2, delta1, delta2
786  sll_int32 :: i, j, nx, ny
787  sll_real64, dimension(:, :):: cx_array
788  class(sll_c_interpolator_2d), pointer :: cxx_2d_interp
789  class(sll_c_interpolator_2d), pointer :: a21_interp
790 
791  do j = 1, ny
792  eta2 = eta2_min + real(j - 1, f64)*delta2
793  do i = 1, nx
794  eta1 = eta1_min + real(i - 1, f64)*delta1
795  cx_array(i, j) = cxx_2d_interp%interpolate_from_interpolant_derivative_eta1(eta1, eta2) + &
796  a21_interp%interpolate_from_interpolant_derivative_eta2(eta1, eta2)
797  end do
798  end do
799  end subroutine coefx_array
800 
801  subroutine coefy_array(eta1_min, eta2_min, &
802  delta1, delta2, nx, ny, cyy_2d_interp, a12_interp, cy_array)
803  implicit none
804  sll_real64 :: eta1, eta1_min, eta2_min
805  sll_real64 :: eta2, delta1, delta2
806  sll_int32 :: i, j, nx, ny
807  sll_real64, dimension(:, :):: cy_array
808  class(sll_c_interpolator_2d), pointer :: cyy_2d_interp
809  class(sll_c_interpolator_2d), pointer :: a12_interp
810 
811  do j = 1, ny
812  eta2 = eta2_min + real(j - 1, f64)*delta2
813  do i = 1, nx
814  eta1 = eta1_min + real(i - 1, f64)*delta1
815  cy_array(i, j) = cyy_2d_interp%interpolate_from_interpolant_derivative_eta2(eta1, eta2) + &
816  a12_interp%interpolate_from_interpolant_derivative_eta1(eta1, eta2)
817  end do
818  end do
819  end subroutine coefy_array
820 
821  subroutine mudpack_curvilinear_cof(x, y, cxx, cyy, cx, cy, ce)
822  real(8) :: x, cxx, cx
823  real(8) :: y, cyy, cy, ce
824  cxx = mudpack_curvilinear_wrapper%cxx_2d_interp%interpolate_from_interpolant_value(x, y)
825  cyy = mudpack_curvilinear_wrapper%cyy_2d_interp%interpolate_from_interpolant_value(x, y)
826  cx = mudpack_curvilinear_wrapper%cx_2d_interp%interpolate_from_interpolant_value(x, y)
827  cy = mudpack_curvilinear_wrapper%cy_2d_interp%interpolate_from_interpolant_value(x, y)
828  ce = mudpack_curvilinear_wrapper%ce_2d_interp%interpolate_from_interpolant_value(x, y)
829  return
830  end subroutine
831 
832  subroutine mudpack_curvilinear_cofcr(x, y, cxx, cxy, cyy, cx, cy, ce)
833  real(8) :: x, cxx, cx, cxy
834  real(8) :: y, cyy, cy, ce
835  cxx = mudpack_curvilinear_wrapper%cxx_2d_interp%interpolate_from_interpolant_value(x, y)
836  cxy = mudpack_curvilinear_wrapper%cxy_2d_interp%interpolate_from_interpolant_value(x, y)
837  cyy = mudpack_curvilinear_wrapper%cyy_2d_interp%interpolate_from_interpolant_value(x, y)
838  cx = mudpack_curvilinear_wrapper%cx_2d_interp%interpolate_from_interpolant_value(x, y)
839  cy = mudpack_curvilinear_wrapper%cy_2d_interp%interpolate_from_interpolant_value(x, y)
840  ce = mudpack_curvilinear_wrapper%ce_2d_interp%interpolate_from_interpolant_value(x, y)
841 
842  return
843  end subroutine
844 
846  subroutine mudpack_curvilinear_bndcr(kbdy, xory, alfa, gbdy)
847  integer :: kbdy
848  real(8) :: xory, alfa, gbdy, x, y, pe, px, py
849  real(8) :: xa, xb, yc, yd, tolmax, relmax
850  common/ftmud2sp/xa, xb, yc, yd, tolmax, relmax
851 
852 !subroutine not used in periodic case
853  if (kbdy == 1) then ! x=xa boundary
854  y = xory
855  x = xa
856  alfa = -1._f64
857  gbdy = px + alfa*pe
858  return
859  end if
860 
861  if (kbdy == 4) then ! y=yd boundary
862  y = yd
863  x = xory
864  alfa = 1._f64
865  gbdy = py + alfa*pe
866  return
867  end if
868  end subroutine
869 
870  subroutine associate_poisson(poisson)
871  type(poisson_2d_mudpack_curvilinear_old), target :: poisson
872 
873  mudpack_curvilinear_wrapper => poisson
874 
875  end subroutine associate_poisson
876 
877 end module sll_m_poisson_2d_mudpack_curvilinear_old
878 
879 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
880 
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_cartesian_meshes
Cartesian mesh basic types.
Definition: sll_m_cartesian_meshes.F90:20
sll_m_coordinate_transformation_2d_base
Abstract class for coordinate transformations.
Definition: sll_m_coordinate_transformation_2d_base.F90:21
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_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_cartesian_meshes::sll_t_cartesian_mesh_2d
2D cartesian mesh
Definition: sll_m_cartesian_meshes.F90:81
sll_m_coordinate_transformation_2d_base::sll_c_coordinate_transformation_2d_base
Abstract class for coordinate transformation.
Definition: sll_m_coordinate_transformation_2d_base.F90:56
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
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1