Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_poisson_2d_tri.F90
Go to the documentation of this file.
1 !File: Module Poisson
2 !
3 ! Poisson solver on unstructured mesh with triangles
4 !
5 module sll_m_poisson_2d_tri
6 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7 #include "sll_errors.h"
8 #include "sll_memory.h"
9 #include "sll_working_precision.h"
10 
11  use sll_m_choleski, only: &
12  sll_s_choles, &
13  sll_s_desrem
14 
15  use sll_m_triangular_meshes, only: &
16  sll_t_triangular_mesh_2d
17 
18  implicit none
19 
20  public :: &
21  sll_s_poisson_2d_triangular_init, &
22  sll_s_poisson_2d_triangular_free, &
23  sll_s_compute_e_from_phi, &
24  sll_s_compute_e_from_rho, &
25  sll_s_compute_phi_from_rho, &
26  sll_t_poisson_2d_triangular
27 
28  private
29 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
30 
31  interface sll_o_create
32  module procedure initialize_poisson_solver
33  module procedure initialize_poisson_solver_from_file
34  end interface sll_o_create
35 
36 ! Mesh parameters
37 !
38 ! nbs - number of nodes
39 ! nbt - number of triangles
40 ! ndiric - number of nodes Dirichlet
41 ! ndirb3 - number of nodes Dirichlet for Ampere
42 ! nbfrax - number of nodes axis/boundary
43 ! nmxfr - number of boundary referenced
44 ! nmxsd - number of subdomains referenced
45 ! nefro - number of elements with one node on the boundary
46 ! nelfr - number of elements on the boundaries
47 ! nelin - number of elements not on the boundaries
48 ! irefdir - boundary references Dirichlet non homogeneous
49 ! nnoeuf - number of nodes Dirichlet non homogeneous
50 !
51 ! Geometry :
52 !
53 ! xlml - limit below x of domain
54 ! xlmu - limit above x of domain
55 ! ylml - limit below y of domain
56 ! ylmu - limit above y of domain
57 ! petitl - small length reference
58 ! grandl - big length reference
59 ! imxref - undefined reference
60 !
61 ! Delaunay-Voronoi:
62 !
63 ! nbcoti - number of edges Delaunay inside
64 ! nbtcot - number of edges Delaunay total
65 ! nbcfli - number of edges inside with wrong CFL
66 ! ncotcu - number of edges fo each boundary
67 ! nnref - number of references Dirichlet non homogeneous boundary
68 !
69 ! Edges on boundaries:
70 !
71 ! nctfrt - number total de edges boundaries
72 ! nctfro - number de edges boundaries par reference
73 ! nctfrp - pointer des arrays de edges boundaries
74 !
75 ! Internal boundaries:
76 !
77 ! nnofnt - number nodes internal boundaries
78 ! ntrfnt - number total triangles
79 ! ntrfrn - number triangles par boundaries
80 ! ntrfrc - number cumsum triangles
81 ! nndfnt - nodes Dirichlet internal boundaries
82 ! ncdfnt - edges Dirichlet internal boundaries
83 !
84 ! Array elements types:
85 !
86 ! nmxcol - number of type
87 ! nclcol - number of elements in every type
88 !
89  sll_real64, parameter :: grandx = 1.d+20
90 
96  type :: sll_t_poisson_2d_triangular
97 
98  private
99  sll_real64, dimension(:), allocatable :: vnx
100  sll_real64, dimension(:), allocatable :: vny
101  logical, dimension(:), allocatable :: naux
102  sll_real64, dimension(:), allocatable :: gradx
103  sll_real64, dimension(:), allocatable :: grady
104  sll_real64, dimension(:), allocatable :: grgr
105  sll_int32, dimension(:), allocatable :: mors1
106  sll_int32, dimension(:), allocatable :: mors2
107  sll_int32, dimension(:), allocatable :: iprof
108  sll_int32, dimension(:), allocatable :: ifron
109  sll_real64, dimension(:), allocatable :: amass
110  sll_real64, dimension(:), allocatable :: vtantx
111  sll_real64, dimension(:), allocatable :: vtanty
112  sll_real64, dimension(:), allocatable :: sv1
113  sll_real64, dimension(:), allocatable :: sv2
114 
115  sll_int32 :: ndiric
116  sll_int32 :: ntypfr(5)
117  sll_real64 :: potfr(5)
118  sll_real64 :: eps0 = 1.0_f64 !8.8542e-12
119 
120  type(sll_t_triangular_mesh_2d), pointer :: mesh => null()
121 
122  end type sll_t_poisson_2d_triangular
123 
124 contains
125 
126  subroutine sll_s_poisson_2d_triangular_init(solver, mesh, ntypfr, potfr)
127 
128  type(sll_t_poisson_2d_triangular) :: solver
129  type(sll_t_triangular_mesh_2d), intent(in) :: mesh
130  sll_int32, dimension(:), intent(in) :: ntypfr
131  sll_real64, dimension(:), intent(in) :: potfr
132 
133  call initialize_poisson_solver(solver, mesh, ntypfr, potfr)
134 
135  end subroutine sll_s_poisson_2d_triangular_init
136 
138  subroutine sll_s_poisson_2d_triangular_free(self)
139  type(sll_t_poisson_2d_triangular) :: self
140 
141 #ifdef DEBUG
142  print *, 'delete poisson solver'
143 #endif
144 
145  end subroutine sll_s_poisson_2d_triangular_free
146 
153  subroutine sll_s_compute_e_from_rho(self, rho, phi, ex, ey)
154  type(sll_t_poisson_2d_triangular) :: self
155  sll_real64, intent(in) :: rho(:)
156  sll_real64, intent(out) :: phi(:)
157  sll_real64, intent(out) :: ex(:)
158  sll_real64, intent(out) :: ey(:)
159 
160  call poissn(self, rho, phi, ex, ey)
161  call poifrc(self, ex, ey)
162 
163  end subroutine sll_s_compute_e_from_rho
164 
165  subroutine sll_s_compute_phi_from_rho(self, rho, phi)
166  type(sll_t_poisson_2d_triangular) :: self
167  sll_real64, intent(in) :: rho(:)
168  sll_real64, intent(out) :: phi(:)
169 
170  call poissn(self, rho, phi)
171 
172  end subroutine sll_s_compute_phi_from_rho
173 
174  subroutine sll_s_compute_e_from_phi(self, phi, ex, ey)
175  type(sll_t_poisson_2d_triangular) :: self
176  sll_real64, intent(in) :: phi(:)
177  sll_real64, intent(out) :: ex(:)
178  sll_real64, intent(out) :: ey(:)
179 
180  call poliss(self, phi, ex, ey)
181  call poifrc(self, ex, ey)
182 
183  end subroutine sll_s_compute_e_from_phi
184 
185  subroutine read_data_solver(ntypfr, potfr)
186 
187  sll_int32, parameter :: nfrmx = 5
188  sll_int32, intent(out) :: ntypfr(nfrmx)
189  sll_real64, intent(out) :: potfr(nfrmx)
190 
191  sll_int32 :: ityp
192  sll_int32 :: ifr
193  sll_int32 :: i
194  character(len=72) :: argv
195  character(len=132) :: inpfil
196  logical :: lask
197  character(len=*), parameter :: self_sub_name = 'read_data_solver'
198 
199  namelist /nlcham/ ntypfr, potfr
200 
201  call get_command_argument(1, argv); write (*, '(1x, a)') argv
202 
203 !------------------------------------------------------------!
204 ! Reads in default parameters from input file (.inp) !
205 ! If LASK=t, ask user if file is to be read. !
206 !------------------------------------------------------------!
207 
208  lask = .true.
209 
210  inpfil = trim(argv)//'.inp'
211 
212  write (*, "(/10x,'Fichier d''entree : ',a)") inpfil
213 
214  open (10, file=inpfil, status='OLD', err=80)
215  write (*, 1050, advance='no') trim(inpfil)
216  lask = .false.
217 
218 80 continue
219 
220  if (lask) then
221  write (*, 1900) inpfil
222  write (*, 1800, advance='no')
223  read (*, "(a)") inpfil
224  write (*, 1700) trim(inpfil)
225  end if
226 
227  write (6, 900)
228 
229  ntypfr = 1
230  potfr = 0.0_f64
231 
232  open (10, file=inpfil)
233  write (*, *) "Lecture de la namelist nlcham"
234  read (10, nlcham, err=100)
235  goto 200
236 100 continue
237  write (*, *) "Error de lecture de la namelist nlcham"
238  write (*, *) "Valeurs par defaut"
239  write (*, nlcham)
240  stop
241 200 close (10)
242 
243  write (6, 921)
244 
245  do i = 1, nfrmx
246  write (6, 922) i, ntypfr(i), potfr(i)
247  end do
248 
249  write (6, 932)
250  do ifr = 1, nfrmx
251  ityp = ntypfr(ifr)
252  if (ityp == 1) then
253  write (6, 902) ifr
254  else if (ityp == 3) then
255  write (6, 904) ifr
256  else
257  write (6, 910) ifr, ityp
258  sll_error(self_sub_name, "Unspecified error.")
259  end if
260  end do
261 
262 900 format(//10x, 'Boundary conditions')
263 902 format(/10x, 'Reference :', i3, 5x, 'Dirichlet ')
264 904 format(/10x, 'Reference :', i3, 5x, 'Neumann')
265 910 format(/10x, 'Option not available'/ &
266  10x, 'Reference :', i3, 5x, 'Type :', i3/)
267 921 format(//5x, 'Frontiere', 5x, 'ntypfr', 7x, 'potfr')
268 922 format(5x, i3, 4x, i10, 2e12.3, 2i10, e12.3)
269 932 format(//10x, 'Boundary conditions'/)
270 1050 format(/' Read settings from file ', a, ' ? Y')
271 1700 format(/' New parameters write to file ', a,/)
272 1800 format(/' Settings may have been changed - New title :')
273 1900 format(/' Input file ', a, ' not found')
274 
275  end subroutine read_data_solver
276 
277  subroutine initialize_poisson_solver_from_file(self, mesh)
278 
279  type(sll_t_poisson_2d_triangular), intent(out) :: self
280  type(sll_t_triangular_mesh_2d), intent(in), target :: mesh
281  sll_int32 :: ntypfr(5)
282  sll_real64 :: potfr(5)
283 
284  call read_data_solver(ntypfr, potfr)
285 
286  call initialize_poisson_solver(self, mesh, ntypfr, potfr)
287 
288  end subroutine initialize_poisson_solver_from_file
289 
290 !Subroutine: init_solveur_poisson
291 ! Allocate arrays
292 !
293 ! amass - matrix mass diagonal
294 ! mors1 - matrix morse2
295 ! mors2 - elements matrices morse
296 ! prof - matrix profil
297 ! grgr - matrix grad-grad
298 ! gradx - matrix gradx
299 ! grady - matrix grady
300 ! fron - nodes Dirichlet
301 ! xaux - arrays auxiliar reals
302 ! iaux - arrays auxiliar integers
303 
304 ! Allocation arrays to stock matrix morse
305 ! Array number of last term of every line (mors1)
306  subroutine initialize_poisson_solver(self, mesh, ntypfr, potfr)
307 
308  type(sll_t_poisson_2d_triangular), intent(out) :: self
309  type(sll_t_triangular_mesh_2d), intent(in), target :: mesh
310  sll_int32, intent(in) :: ntypfr(:)
311  sll_real64, intent(in) :: potfr(:)
312 
313  sll_real64, allocatable :: tmp1(:)
314  character(len=*), parameter :: self_sub_name = 'read_data_solver'
315  character(len=128) :: err_msg
316 
317  sll_int32 :: nref, nn, ndir
318  sll_int32 :: i
319  sll_int32 :: ierr
320 
321  if (mesh%analyzed) then
322  self%mesh => mesh
323  else
324  err_msg = "Call sll_s_analyze_triangular_mesh before initialize_poisson_solver."
325  sll_error(self_sub_name, err_msg)
326  end if
327 
328  self%ntypfr = ntypfr
329  self%potfr = potfr
330 
331  allocate (self%sv1(mesh%nbtcot)); self%sv1 = 0.0_f64
332  allocate (self%sv2(mesh%nbtcot)); self%sv2 = 0.0_f64
333  allocate (self%vtantx(mesh%nbtcot)); self%vtantx = 0.0_f64
334  allocate (self%vtanty(mesh%nbtcot)); self%vtanty = 0.0_f64
335 
336  allocate (self%mors1(mesh%num_nodes + 1)); self%mors1 = 0
337 
338 ! Array witth number of every line term (mors2)
339 
340  allocate (self%mors2(12*mesh%num_nodes)); self%mors2 = 0
341 
342 ! mors1,mors2.
343 
344  call morse(mesh%npoel1, &
345  mesh%npoel2, &
346  mesh%nodes, &
347  mesh%num_triangles, &
348  mesh%num_nodes, &
349  self%mors1, &
350  self%mors2)
351 
352  allocate (self%iprof(mesh%num_nodes + 1)); self%iprof = 0
353 
354  call profil(mesh%nodes, &
355  mesh%num_nodes, &
356  mesh%npoel1, &
357  mesh%npoel2, &
358  self%iprof)
359 
360 !======================================================================
361 !--- 2.0 --- POISSON finite element -----------------
362 !======================================================================
363 
364 !matrix de masse diagonalized
365  allocate (self%amass(mesh%num_nodes)); self%amass = 0.0_f64
366 
367 !matrix "grad-grad" stocked profil form.
368  allocate (self%grgr(self%iprof(mesh%num_nodes + 1))); self%grgr = 0.0_f64
369 
370 !gradx grady
371  allocate (self%gradx(self%mors1(mesh%num_nodes + 1))); self%gradx = 0.0_f64
372  allocate (self%grady(self%mors1(mesh%num_nodes + 1))); self%grady = 0.0_f64
373 
374 !--- Array boundaries Dirichlet -------------------------
375 
376  ndir = 0
377  do nn = 1, mesh%num_nodes
378  nref = mesh%refs(nn)
379  if (nref > 0) then
380  if (self%ntypfr(nref) == 1) then
381  ndir = ndir + 1
382  end if
383  end if
384  end do
385 
386  allocate (self%ifron(ndir)); self%ifron = 0
387 
388  ndir = 0
389  do nn = 1, mesh%num_nodes
390  nref = mesh%refs(nn)
391  if (nref > 0) then
392  if (self%ntypfr(nref) == 1) then
393  ndir = ndir + 1
394  self%ifron(ndir) = nn
395  end if
396  end if
397  end do
398 
399  self%ndiric = ndir
400 
401 !--- Calcul des matrix ----------------------------------------------
402 
403  call poismc(self)
404 
405 !Calcul de la matrix B tel que B*Bt = A dans le cas Cholesky
406 
407  allocate (tmp1(self%iprof(mesh%num_nodes + 1))); tmp1 = 0.0_f64
408 
409 !write(*,"(//5x,a)")" *** Appel Choleski pour Poisson *** "
410  call sll_s_choles(self%iprof, self%grgr, tmp1)
411 
412  do i = 1, self%iprof(mesh%num_nodes + 1)
413  self%grgr(i) = tmp1(i)
414  end do
415 
416  deallocate (tmp1)
417 
418  sll_clear_allocate(self%vnx(1:self%mesh%num_nodes), ierr)
419  sll_clear_allocate(self%vny(1:self%mesh%num_nodes), ierr)
420  sll_allocate(self%naux(1:self%mesh%num_nodes), ierr)
421  self%naux = .false.
422 
423  end subroutine initialize_poisson_solver
424 
425 !======================================================================
426 
427 ! Function: morse
428 ! Calculation of the arrays containing the address of the last term of
429 ! each line of the "walrus" matrix and the arrays containing the
430 ! numbers of the terms of these matrixes. The diagonal term is
431 ! last place in each line.
432 !
433 ! npoel1 - location in npoel2 of the last element relative to each
434 ! node with npoel1(1)=0 and npoel1(i+1) relative to node i
435 ! npoel2 - arrays numbers of elements having a common vertex
436 ! ntri - triangles vertex numbers
437 ! nbt - number of mesh triangles
438 ! nbs - number of nodes
439 !
440 ! mors1 - arrays addresses of the last terms of each line
441 ! the matrix with the convention:
442 ! mors1(1)=0 and mors1(i+1) aii address
443 ! mors2 - arrays of the numbers of the "morse" matrix terms
444 
445  subroutine morse(npoel1, npoel2, ntri, nbt, nbs, mors1, mors2)
446 
447  sll_int32, intent(in) :: nbt
448  sll_int32, intent(in) :: nbs
449  sll_int32, dimension(:), intent(in) :: npoel1
450  sll_int32, dimension(:), intent(in) :: npoel2
451  sll_int32, dimension(3, nbt), intent(in) :: ntri
452  sll_int32, dimension(:), intent(out) :: mors1
453  sll_int32, dimension(:), intent(out) :: mors2
454  sll_int32, dimension(20) :: ilign
455 
456  sll_int32 :: l, itest1, itest2, js1, js2, is1, is2, is3, numel
457  sll_int32 :: iel, nlign, nel, is, im, k
458 
459  im = 0
460  k = 0
461 
462  mors1(1) = 0
463 
464  do is = 1, nbs !Boucle sur les nodes
465 
466  !elements 'is' node
467 
468  nel = npoel1(is + 1) - npoel1(is)
469  nlign = 0
470 
471  !Loop over this elements
472 
473  do iel = 1, nel
474 
475  k = k + 1
476 
477  numel = npoel2(k)
478  is1 = ntri(1, numel); is2 = ntri(2, numel); is3 = ntri(3, numel)
479  if (is1 == is) then
480  js1 = is2; js2 = is3
481  end if
482  if (is2 == is) then
483  js1 = is3; js2 = is1
484  end if
485  if (is3 == is) then
486  js1 = is1; js2 = is2
487  end if
488 
489  !We see if the 2 nodes other than is of the current element
490  !have already been taken into account in all the interacting nodes
491  !with is.
492 
493  itest1 = 0; itest2 = 0
494  if (nlign .NE. 0) then
495  do l = 1, nlign
496  if (js1 == ilign(l)) then
497  itest1 = 1
498  end if
499  if (js2 == ilign(l)) then
500  itest2 = 1
501  end if
502  end do
503  end if
504 
505  if (itest1 == 0) then
506  nlign = nlign + 1
507  ilign(nlign) = js1
508  end if
509  if (itest2 == 0) then
510  nlign = nlign + 1
511  ilign(nlign) = js2
512  end if
513 
514  end do
515 
516  !Definition of the address of the last term of the line
517  mors1(is + 1) = mors1(is) + nlign + 1
518 
519  !Filling the jaw arrays2 with the term numbers
520  !of the is line.
521 
522  if (nlign .NE. 0) then
523  do l = 1, nlign
524  im = im + 1
525  mors2(im) = ilign(l)
526  end do
527  end if
528  im = im + 1
529  mors2(im) = is
530  end do
531 
532  end subroutine morse
533 
534 ! ======================================================================
535 
536 !Function: poismc
537 !
538 ! Compute all matrices for cartesian Poisson
539 !
540 ! Input:
541 !
542 ! m - super-arrays
543 ! coor - coord nodes
544 ! refs - indice if node in ref boundary
545 ! ifron - arrays of nodes Dirichlet
546 ! mors1 - arrays of number terms per line matrix morse symetric
547 ! mors2 - arrays of number terms every line matrix morse symetric
548 ! ntri - numbers triangles nodes
549 ! area - area triangles
550 ! iprof - profile matrix grad-grad
551 !
552 ! noefnt - nodes interns Dirichlet
553 !
554 ! Output:
555 !
556 ! agrgr - matrix operator "grad-grad" "morse"
557 ! aliss - matrix fitting
558 ! amass - matrix mass diag
559 ! amclt - matrix mass diag absorbing bc
560 ! aroro - matrix operator "rot-rot" form "morse"
561 ! d1dx, - deriv functions basis every triangle
562 !..d3dy
563 ! gradx - matrix operator gradx
564 ! grady - matrix operator grady
565 ! rotx - matrix operator rotx
566 ! roty - matrix operator roty
567 ! area - every element
568 !
569 !Auteur: Puertolas - Version 1.0 Octobre 1992
570  subroutine poismc(self)
571 
572  type(sll_t_poisson_2d_triangular), intent(inout) :: self
573  sll_real64 :: amloc(3), aggloc(9), grxloc(9), gryloc(9)
574  sll_real64 :: dntx1, dntx2, dntx3, dnty1, dnty2, dnty3
575  sll_real64 :: x1t, x2t, x3t, y1t, y2t, y3t, coef
576  sll_int32 :: is1t, is2t, is3t, iel
577  sll_int32 :: is, j
578 
579 !Boucle sur les elements.
580 
581  do iel = 1, self%mesh%num_triangles
582 
583  !coefficients geometry
584 
585  is1t = self%mesh%nodes(1, iel)
586  is2t = self%mesh%nodes(2, iel)
587  is3t = self%mesh%nodes(3, iel)
588 
589  x1t = self%mesh%coord(1, is1t)
590  x2t = self%mesh%coord(1, is2t)
591  x3t = self%mesh%coord(1, is3t)
592 
593  y1t = self%mesh%coord(2, is1t)
594  y2t = self%mesh%coord(2, is2t)
595  y3t = self%mesh%coord(2, is3t)
596 
597  dntx1 = y2t - y3t
598  dntx2 = y3t - y1t
599  dntx3 = y1t - y2t
600 
601  dnty1 = x3t - x2t
602  dnty2 = x1t - x3t
603  dnty3 = x2t - x1t
604 
605  !Contribution matrix mass
606 
607  amloc(1) = self%mesh%area(iel)/3.
608  amloc(2) = self%mesh%area(iel)/3.
609  amloc(3) = self%mesh%area(iel)/3.
610 
611  !Assembling
612 
613  call asbld(amloc, is1t, is2t, is3t, self%amass)
614 
615  !Contribution matrix grad-grad
616 
617  coef = 1./(4.*self%mesh%area(iel))
618 
619  aggloc(1) = (dntx1**2 + dnty1**2)*coef
620  aggloc(2) = (dntx1*dntx2 + dnty1*dnty2)*coef
621  aggloc(3) = (dntx1*dntx3 + dnty1*dnty3)*coef
622  aggloc(4) = (dntx2**2 + dnty2**2)*coef
623  aggloc(5) = (dntx2*dntx3 + dnty2*dnty3)*coef
624  aggloc(6) = (dntx3**2 + dnty3**2)*coef
625 
626  !Contribution matrix gradx grady:
627 
628  grxloc(1) = -dntx1/6.; gryloc(1) = -dnty1/6.
629  grxloc(2) = -dntx2/6.; gryloc(2) = -dnty2/6.
630  grxloc(3) = -dntx3/6.; gryloc(3) = -dnty3/6.
631  grxloc(4) = -dntx1/6.; gryloc(4) = -dnty1/6.
632  grxloc(5) = -dntx2/6.; gryloc(5) = -dnty2/6.
633  grxloc(6) = -dntx3/6.; gryloc(6) = -dnty3/6.
634  grxloc(7) = -dntx1/6.; gryloc(7) = -dnty1/6.
635  grxloc(8) = -dntx2/6.; gryloc(8) = -dnty2/6.
636  grxloc(9) = -dntx3/6.; gryloc(9) = -dnty3/6.
637 
638  !Assembling
639 
640  call asblp(self%iprof, aggloc, is1t, is2t, is3t, self%grgr)
641 
642  call asblm2(grxloc, &
643  gryloc, &
644  self%mors1, &
645  self%mors2, &
646  is1t, &
647  is2t, &
648  is3t, &
649  self%gradx, &
650  self%grady)
651 
652  end do
653 
654 ! ======================================================================
655 ! ... Dirichlet
656 
657  do j = 1, self%ndiric
658  is = self%ifron(j)
659  self%grgr(self%iprof(is + 1)) = grandx
660  end do
661 
662  end subroutine poismc
663 
664 !Function: poissn
665 !potential and electric field
666 !Poisson en cartesien
667 !
668 ! ex - e1 nodes
669 ! ey - e2 nodes
670 ! rho - density charge nodes
671 ! phi - potential nodes
672 ! ifron - nodes boundary Dirichlet
673 ! noefnt - nodes internal boundaries Dirichlet
674 ! irffnt - reference of this nodes
675 ! grgr - matrix operator "grad-grad"
676 ! amass - matrix mass diag
677 ! gra1 - matrix operator "grad1"
678 ! gra2 - matrix operator "grad2"
679 ! mors1 - arrays matrix "morse"
680 ! mors2 - arrays matrix "morse"
681 ! iprof - matrix profil
682 !
683 !
684 ! grgrdd - product matrix grgr vector p (gradient conjugue)
685 ! dird - direction method gradient conjugue
686 ! res - residual method gradient conjugue
687 ! sdmb - second member system potential
688 ! sdmb12 - second member system electric field
689 ! precd - arrays local utilise gradient conjugue preconditionned
690 ! nbs - number nodes mesh
691 ! nmxfr - number boundaries ref
692 ! errpoi - precision
693 ! nitgc - number max iterations gradient conjugue
694 !
695  subroutine poissn(self, rho, phi, ex, ey)
696 
697  type(sll_t_poisson_2d_triangular) :: self
698  sll_real64, dimension(:), intent(in) :: rho
699  sll_real64, dimension(:), intent(out) :: phi
700  sll_real64, dimension(:), intent(out), optional :: ex
701  sll_real64, dimension(:), intent(out), optional :: ey
702 
703  sll_real64, dimension(:), allocatable :: sdmb
704  sll_real64, dimension(:), allocatable :: sdmb12
705 
706  sll_int32 :: i, is, nref, ierr
707 
708 ! ----------- CALCUL DU POTENTIEL -------------------------------------
709 !
710 !... RHS ...
711 
712  sll_allocate(sdmb(self%mesh%num_nodes), ierr)
713 
714 !!$OMP PARALLEL DEFAULT(SHARED), PRIVATE(is)
715 !!$OMP DO SCHEDULE(RUNTIME)
716  do is = 1, self%mesh%num_nodes
717  sdmb(is) = self%amass(is)*rho(is)/self%eps0
718  end do
719 !!$OMP END DO
720 !!$OMP END PARALLEL
721 
722 !... Dirichlet
723 
724  do is = 1, self%ndiric
725  nref = self%mesh%refs(self%ifron(is))
726  sdmb(self%ifron(is)) = self%potfr(nref)*grandx
727  end do
728 
729  call sll_s_desrem(self%iprof, self%grgr, sdmb, self%mesh%num_nodes, phi)
730 
731 !*** Ex Ey:
732 
733  if (present(ex) .and. present(ey)) then
734 
735  sll_allocate(sdmb12(self%mesh%num_nodes), ierr)
736  !*** RHS for Ex:
737 
738  call m1p(self%gradx, self%mors1, self%mors2, phi, self%mesh%num_nodes, sdmb12)
739 
740  do i = 1, self%mesh%num_nodes
741  ex(i) = sdmb12(i)/self%amass(i)
742  end do
743 
744  !*** RHS for Ey:
745 
746  call m1p(self%grady, self%mors1, self%mors2, phi, self%mesh%num_nodes, sdmb12)
747 
748  do i = 1, self%mesh%num_nodes
749  ey(i) = sdmb12(i)/self%amass(i)
750  end do
751 
752  end if
753 
754  end subroutine poissn
755 
756 !Function: poifrc
757 !Correction on boundaries
758 !
759 ! - E.tau = 0 boundaries Dirichlets
760 ! - E.nu = 0 boundaries Neumann
761 !
762 !Input:
763 !
764 ! ksofro - numbers edge nodes
765 ! krefro - numbers edge reference
766 ! vnofro - vector normal (inside)
767 ! vnx,vny - arrays normal vectors nodes
768 !
769 !Array local:
770 !
771 ! naux - arrays for nodes on boundary
772 !
773  subroutine poifrc(self, ex, ey)
774 
775  type(sll_t_poisson_2d_triangular) :: self
776  sll_real64 :: ex(:), ey(:)
777  sll_real64 :: pscal, xnor
778  sll_int32 :: is1, is2, ict
779  sll_int32 :: i
780 
781 !!$ ======================================================================
782 !!$ ... E.tau = 0 boundaries Dirichlet
783 
784 !!$ ... loop over edges boundaries for normals
785 !!$ ... nodes "Dirichlet"
786  self%vnx = 0.0_f64; self%vny = 0.0_f64; self%naux = .false.
787 
788  do ict = 1, self%mesh%nctfrt
789 
790  if (self%ntypfr(self%mesh%krefro(ict)) == 1) then
791 
792  is1 = self%mesh%ksofro(1, ict)
793  is2 = self%mesh%ksofro(2, ict)
794 
795  self%vnx(is1) = self%vnx(is1) + self%mesh%vnofro(1, ict)
796  self%vny(is1) = self%vny(is1) + self%mesh%vnofro(2, ict)
797  self%vnx(is2) = self%vnx(is2) + self%mesh%vnofro(1, ict)
798  self%vny(is2) = self%vny(is2) + self%mesh%vnofro(2, ict)
799 
800  self%naux(is1) = .true.
801  self%naux(is2) = .true.
802 
803  end if
804 
805  end do
806 
807 !... E.tau=0
808 
809  do i = 1, self%mesh%num_nodes
810 
811  if (self%naux(i)) then
812 
813  xnor = sqrt(self%vnx(i)**2 + self%vny(i)**2)
814  if (xnor > self%mesh%petitl) then
815  self%vnx(i) = self%vnx(i)/xnor
816  self%vny(i) = self%vny(i)/xnor
817 
818  pscal = self%vnx(i)*ex(i) + self%vny(i)*ey(i)
819  ex(i) = self%vnx(i)*pscal
820  ey(i) = self%vny(i)*pscal
821  end if
822 
823  end if
824 
825  end do
826 
827 !======================================================================
828 ! ... E.nu = 0 boundaries Neumann
829 !
830  self%vnx = 0.0_f64; self%vny = 0.0_f64; self%naux = .false.
831 
832 ! ... Loop over edges boundaries nodes "Neumann"
833 
834  do ict = 1, self%mesh%nctfrt
835 
836  if (self%ntypfr(self%mesh%krefro(ict)) == 3) then
837 
838  is1 = self%mesh%ksofro(1, ict)
839  is2 = self%mesh%ksofro(2, ict)
840 
841  self%vnx(is1) = self%vnx(is1) + self%mesh%vnofro(1, ict)
842  self%vny(is1) = self%vny(is1) + self%mesh%vnofro(2, ict)
843  self%vnx(is2) = self%vnx(is2) + self%mesh%vnofro(1, ict)
844  self%vny(is2) = self%vny(is2) + self%mesh%vnofro(2, ict)
845 
846  self%naux(is1) = .true.
847  self%naux(is2) = .true.
848 
849  end if
850 
851  end do
852 
853 !... E.nu=0
854 
855  do i = 1, self%mesh%num_nodes
856 
857  if (self%naux(i)) then
858 
859  xnor = sqrt(self%vnx(i)**2 + self%vny(i)**2)
860  if (xnor > self%mesh%petitl) then
861  self%vnx(i) = self%vnx(i)/xnor
862  self%vny(i) = self%vny(i)/xnor
863  pscal = self%vnx(i)*ex(i) + self%vny(i)*ey(i)
864  ex(i) = ex(i) - self%vnx(i)*pscal
865  ey(i) = ey(i) - self%vny(i)*pscal
866  end if
867 
868  end if
869 
870  end do
871 
872  end subroutine poifrc
873 
874 !Function: profil
875 !numbers elements on diagonal matrix profil associate
876 !
877 ! Input:
878 ! npoel1 - index in npoel2 of last element of every node
879 ! with npoel1(1)=0 and npoel1(i+1) link to node i
880 ! npoel2 - arrays numbers elements with common node
881 ! ntri - numbers nodes triangles
882 ! nbt - number triangles
883 ! noe - number nodes
884 !
885 ! Ouput:
886 ! iprof - arrays numbers terms diagonal with
887 ! - iprof(i+1) term diag line i
888 ! - iprof(1)=0
889 !
890 !Auteur:
891 ! J. Segre - Juillet 89
892  subroutine profil(ntri, nbs, npoel1, npoel2, iprof)
893 
894  sll_int32, dimension(:) :: iprof
895  sll_int32, dimension(:), intent(in) :: npoel1, npoel2
896  sll_int32, dimension(:, :), intent(in) :: ntri
897  sll_int32 :: in, k, is1, is2, is3, numel, ind, iel, nel
898  sll_int32, intent(in) :: nbs
899 
900 !************************************************************************
901 !*
902 !* Compute length of line i matrix profil :
903 !* for every node i we look for nodes j linked;
904 !* if node j not found, length line j is set to
905 !* j-i+1
906 !
907 !* loop over nodes
908 
909  k = 0
910  do in = 1, nbs
911 
912  !* number elements with in node
913  nel = npoel1(in + 1) - npoel1(in)
914 
915  !* loop over these elements
916 
917  do iel = 1, nel
918 
919  k = k + 1
920  !* number element
921  numel = npoel2(k)
922  is1 = ntri(1, numel)
923  ind = is1 + 1
924  if (iprof(ind) .eq. 0) then
925  iprof(ind) = is1 - in + 1
926  end if
927  is2 = ntri(2, numel)
928  ind = is2 + 1
929  if (iprof(ind) .eq. 0) then
930  iprof(ind) = is2 - in + 1
931  end if
932  is3 = ntri(3, numel)
933  ind = is3 + 1
934  if (iprof(ind) .eq. 0) then
935  iprof(ind) = is3 - in + 1
936  end if
937 
938  end do
939 
940  end do
941 
942 !* compute position terms diags matrix
943 !* profil (sum numbers elements of lines
944 !* before and current line ).
945 
946  do ind = 3, nbs + 1
947  iprof(ind) = iprof(ind - 1) + iprof(ind)
948  end do
949 
950  end subroutine profil
951 
952 !Function: asbld
953 ! Assemble matrix element in matrix global if diagonal
954 !
955 !Parametres d'entree:
956 !
957 ! aele - Matrix diag
958 ! (3 terms per element triangular)
959 ! i1,i2,i3 - numbers nodes element
960 !
961 !Parametre resultat:
962 !
963 ! xmass - matrix global diagonal
964 !
965 !Auteur:
966 ! J. Segre - Version 1.0 Juillet 1989
967  subroutine asbld(aele, i1, i2, i3, xmass)
968 
969  sll_real64, dimension(:) :: aele, xmass
970  sll_int32 :: i1, i2, i3
971 
972  xmass(i1) = xmass(i1) + aele(1)
973  xmass(i2) = xmass(i2) + aele(2)
974  xmass(i3) = xmass(i3) + aele(3)
975 
976  end subroutine
977 
978 !Function: asblm2
979 !
980 ! assemble 3 matrix element
981 ! 3 matrix global stocked
982 ! form "morse" non symetric
983 !
984 ! Input:
985 ! aele1,aele2 - Matrices elements
986 ! mors1 - Array number of terms per
987 ! line matrix morse
988 ! mors2 - Array numbers terms
989 ! of every line matrix morse
990 ! i1,i2,i3 - Num nodes element
991 !
992 ! Output:
993 ! a1,a2 - Matrix global stocked
994 ! format "morse"
995 !
996 !Auteur:
997 !J. Segre - Version 1.0 Aout 1989
998  subroutine asblm2(aele1, aele2, mors1, mors2, i1, i2, i3, a1, a2)
999 
1000  sll_int32, intent(in) :: i1, i2, i3
1001  sll_int32, dimension(:), intent(in) :: mors1, mors2
1002  sll_real64, dimension(:), intent(in) :: aele1, aele2
1003  sll_real64, dimension(:), intent(out) :: a1, a2
1004  sll_int32 :: j, j1, j2, ind1, ind2, ind3
1005 
1006 ! --- 1.1 --- terms diag ---------------------------
1007 
1008  ind1 = mors1(i1 + 1)
1009  a1(ind1) = a1(ind1) + aele1(1)
1010  a2(ind1) = a2(ind1) + aele2(1)
1011 
1012  ind2 = mors1(i2 + 1)
1013  a1(ind2) = a1(ind2) + aele1(5)
1014  a2(ind2) = a2(ind2) + aele2(5)
1015 
1016  ind3 = mors1(i3 + 1)
1017  a1(ind3) = a1(ind3) + aele1(9)
1018  a2(ind3) = a2(ind3) + aele2(9)
1019 
1020 ! --- 1.2 --- other terms ------------------------------
1021 
1022  j2 = ind1 - 1
1023  j1 = mors1(i1) + 1
1024 
1025  if (j2 >= j1) then
1026  do j = j1, j2
1027  if (i2 == mors2(j)) then
1028  a1(j) = a1(j) + aele1(2)
1029  a2(j) = a2(j) + aele2(2)
1030  end if
1031  if (i3 == mors2(j)) then
1032  a1(j) = a1(j) + aele1(3)
1033  a2(j) = a2(j) + aele2(3)
1034  end if
1035  end do
1036  end if
1037 
1038  j2 = ind2 - 1
1039  j1 = mors1(i2) + 1
1040 
1041  if (j2 >= j1) then
1042  do j = j1, j2
1043  if (i1 == mors2(j)) then
1044  a1(j) = a1(j) + aele1(4)
1045  a2(j) = a2(j) + aele2(4)
1046  end if
1047  if (i3 == mors2(j)) then
1048  a1(j) = a1(j) + aele1(6)
1049  a2(j) = a2(j) + aele2(6)
1050  end if
1051  end do
1052  end if
1053 
1054  j2 = ind3 - 1
1055  j1 = mors1(i3) + 1
1056 
1057  if (j2 >= j1) then
1058  do j = j1, j2
1059  if (i1 == mors2(j)) then
1060  a1(j) = a1(j) + aele1(7)
1061  a2(j) = a2(j) + aele2(7)
1062  end if
1063  if (i2 == mors2(j)) then
1064  a1(j) = a1(j) + aele1(8)
1065  a2(j) = a2(j) + aele2(8)
1066  end if
1067  end do
1068  end if
1069 
1070  end subroutine
1071 
1072 !----------------------------------------------------------------------
1073 
1074 !Function: asblp
1075 !
1076 ! assemble matrix element
1077 ! in matrix global symetric
1078 ! format "profil"
1079 !
1080 !Input:
1081 !
1082 ! aele - Matrix element (6 - 10 terms
1083 ! element triangle or quad)
1084 ! iprof - Description matrix
1085 ! i1,i2,i3 - Num nodes element
1086 !
1087 !Output:
1088 !
1089 ! xmass - matrix global stockee format "profil"
1090 !
1091 !Auteur:
1092 ! J. Segre - Version 1.0 Aout 1989
1093  subroutine asblp(iprof, aele, i1, i2, i3, xmass)
1094 
1095  sll_int32, dimension(:) :: iprof
1096  sll_real64, dimension(:) :: aele(*), xmass(*)
1097  sll_int32 :: i1, i2, i3, idiag1, idiag2, idiag3, ind
1098 
1099 !--- 1.1 --- terms diag ---------------------------
1100 
1101  idiag1 = iprof(i1 + 1)
1102  idiag2 = iprof(i2 + 1)
1103  idiag3 = iprof(i3 + 1)
1104 
1105  xmass(idiag1) = xmass(idiag1) + aele(1)
1106  xmass(idiag2) = xmass(idiag2) + aele(4)
1107  xmass(idiag3) = xmass(idiag3) + aele(6)
1108 
1109 !--- 1.2 --- others terms ------------------------------
1110 ! (save only aij if i>j)
1111 
1112  if (i1 < i2) then
1113  ind = idiag2 + i1 - i2
1114  else
1115  ind = idiag1 + i2 - i1
1116  end if
1117 
1118  xmass(ind) = xmass(ind) + aele(2)
1119 
1120  if (i1 < i3) then
1121  ind = idiag3 + i1 - i3
1122  else
1123  ind = idiag1 + i3 - i1
1124  end if
1125 
1126  xmass(ind) = xmass(ind) + aele(3)
1127 
1128  if (i2 < i3) then
1129  ind = idiag3 + i2 - i3
1130  else
1131  ind = idiag2 + i3 - i2
1132  end if
1133 
1134  xmass(ind) = xmass(ind) + aele(5)
1135 
1136 !----------------------------------------------------------------------
1137 
1138  end subroutine asblp
1139 
1140 !Function: m1p
1141 !
1142 ! yvect = xmors.xvect
1143 ! xvect,yvect vectors
1144 !
1145 !Input:
1146 ! xmors - matrix "morse"
1147 ! xvect - vector
1148 ! mors1,mors2 - arrays matrix "morse"
1149 ! nlign - number lines matrices
1150 !
1151 !Output:
1152 ! yvect - vector result
1153 !
1154 !Auteur:
1155 ! J. Segre - Version 1.0 Decembre 1989
1156  subroutine m1p(xmors, mors1, mors2, xvect, nlign, yvect)
1157 
1158  sll_real64 :: xmors(*), xvect(*), yvect(*)
1159  sll_int32 :: mors1(*), mors2(*)
1160  sll_int32 :: il, nlign, noeui
1161 
1162  do il = 1, nlign
1163 
1164  noeui = mors1(il + 1) - mors1(il)
1165 
1166  select case (noeui)
1167  case (6)
1168  yvect(il) = xmors(mors1(il + 1))*xvect(mors2(mors1(il + 1))) + &
1169  xmors(mors1(il + 1) - 1)*xvect(mors2(mors1(il + 1) - 1)) + &
1170  xmors(mors1(il + 1) - 2)*xvect(mors2(mors1(il + 1) - 2)) + &
1171  xmors(mors1(il + 1) - 3)*xvect(mors2(mors1(il + 1) - 3)) + &
1172  xmors(mors1(il + 1) - 4)*xvect(mors2(mors1(il + 1) - 4)) + &
1173  xmors(mors1(il + 1) - 5)*xvect(mors2(mors1(il + 1) - 5))
1174 
1175  case (5)
1176  yvect(il) = xmors(mors1(il + 1))*xvect(mors2(mors1(il + 1))) + &
1177  xmors(mors1(il + 1) - 1)*xvect(mors2(mors1(il + 1) - 1)) + &
1178  xmors(mors1(il + 1) - 2)*xvect(mors2(mors1(il + 1) - 2)) + &
1179  xmors(mors1(il + 1) - 3)*xvect(mors2(mors1(il + 1) - 3)) + &
1180  xmors(mors1(il + 1) - 4)*xvect(mors2(mors1(il + 1) - 4))
1181 
1182  case (7)
1183  yvect(il) = xmors(mors1(il + 1))*xvect(mors2(mors1(il + 1))) + &
1184  xmors(mors1(il + 1) - 1)*xvect(mors2(mors1(il + 1) - 1)) + &
1185  xmors(mors1(il + 1) - 2)*xvect(mors2(mors1(il + 1) - 2)) + &
1186  xmors(mors1(il + 1) - 3)*xvect(mors2(mors1(il + 1) - 3)) + &
1187  xmors(mors1(il + 1) - 4)*xvect(mors2(mors1(il + 1) - 4)) + &
1188  xmors(mors1(il + 1) - 5)*xvect(mors2(mors1(il + 1) - 5)) + &
1189  xmors(mors1(il + 1) - 6)*xvect(mors2(mors1(il + 1) - 6))
1190 
1191  case (4)
1192  yvect(il) = xmors(mors1(il + 1))*xvect(mors2(mors1(il + 1))) + &
1193  xmors(mors1(il + 1) - 1)*xvect(mors2(mors1(il + 1) - 1)) + &
1194  xmors(mors1(il + 1) - 2)*xvect(mors2(mors1(il + 1) - 2)) + &
1195  xmors(mors1(il + 1) - 3)*xvect(mors2(mors1(il + 1) - 3))
1196 
1197  case (8)
1198  yvect(il) = xmors(mors1(il + 1))*xvect(mors2(mors1(il + 1))) + &
1199  xmors(mors1(il + 1) - 1)*xvect(mors2(mors1(il + 1) - 1)) + &
1200  xmors(mors1(il + 1) - 2)*xvect(mors2(mors1(il + 1) - 2)) + &
1201  xmors(mors1(il + 1) - 3)*xvect(mors2(mors1(il + 1) - 3)) + &
1202  xmors(mors1(il + 1) - 4)*xvect(mors2(mors1(il + 1) - 4)) + &
1203  xmors(mors1(il + 1) - 5)*xvect(mors2(mors1(il + 1) - 5)) + &
1204  xmors(mors1(il + 1) - 6)*xvect(mors2(mors1(il + 1) - 6)) + &
1205  xmors(mors1(il + 1) - 7)*xvect(mors2(mors1(il + 1) - 7))
1206 
1207  case (3)
1208  yvect(il) = xmors(mors1(il + 1))*xvect(mors2(mors1(il + 1))) + &
1209  xmors(mors1(il + 1) - 1)*xvect(mors2(mors1(il + 1) - 1)) + &
1210  xmors(mors1(il + 1) - 2)*xvect(mors2(mors1(il + 1) - 2))
1211 
1212  case (9)
1213  yvect(il) = xmors(mors1(il + 1))*xvect(mors2(mors1(il + 1))) + &
1214  xmors(mors1(il + 1) - 1)*xvect(mors2(mors1(il + 1) - 1)) + &
1215  xmors(mors1(il + 1) - 2)*xvect(mors2(mors1(il + 1) - 2)) + &
1216  xmors(mors1(il + 1) - 3)*xvect(mors2(mors1(il + 1) - 3)) + &
1217  xmors(mors1(il + 1) - 4)*xvect(mors2(mors1(il + 1) - 4)) + &
1218  xmors(mors1(il + 1) - 5)*xvect(mors2(mors1(il + 1) - 5)) + &
1219  xmors(mors1(il + 1) - 6)*xvect(mors2(mors1(il + 1) - 6)) + &
1220  xmors(mors1(il + 1) - 7)*xvect(mors2(mors1(il + 1) - 7)) + &
1221  xmors(mors1(il + 1) - 8)*xvect(mors2(mors1(il + 1) - 8))
1222 
1223  case (10)
1224  yvect(il) = xmors(mors1(il + 1))*xvect(mors2(mors1(il + 1))) + &
1225  xmors(mors1(il + 1) - 1)*xvect(mors2(mors1(il + 1) - 1)) + &
1226  xmors(mors1(il + 1) - 2)*xvect(mors2(mors1(il + 1) - 2)) + &
1227  xmors(mors1(il + 1) - 3)*xvect(mors2(mors1(il + 1) - 3)) + &
1228  xmors(mors1(il + 1) - 4)*xvect(mors2(mors1(il + 1) - 4)) + &
1229  xmors(mors1(il + 1) - 5)*xvect(mors2(mors1(il + 1) - 5)) + &
1230  xmors(mors1(il + 1) - 6)*xvect(mors2(mors1(il + 1) - 6)) + &
1231  xmors(mors1(il + 1) - 7)*xvect(mors2(mors1(il + 1) - 7)) + &
1232  xmors(mors1(il + 1) - 8)*xvect(mors2(mors1(il + 1) - 8)) + &
1233  xmors(mors1(il + 1) - 9)*xvect(mors2(mors1(il + 1) - 9))
1234 
1235  case (11)
1236  yvect(il) = xmors(mors1(il + 1))*xvect(mors2(mors1(il + 1))) + &
1237  xmors(mors1(il + 1) - 1)*xvect(mors2(mors1(il + 1) - 1)) + &
1238  xmors(mors1(il + 1) - 2)*xvect(mors2(mors1(il + 1) - 2)) + &
1239  xmors(mors1(il + 1) - 3)*xvect(mors2(mors1(il + 1) - 3)) + &
1240  xmors(mors1(il + 1) - 4)*xvect(mors2(mors1(il + 1) - 4)) + &
1241  xmors(mors1(il + 1) - 5)*xvect(mors2(mors1(il + 1) - 5)) + &
1242  xmors(mors1(il + 1) - 6)*xvect(mors2(mors1(il + 1) - 6)) + &
1243  xmors(mors1(il + 1) - 7)*xvect(mors2(mors1(il + 1) - 7)) + &
1244  xmors(mors1(il + 1) - 8)*xvect(mors2(mors1(il + 1) - 8)) + &
1245  xmors(mors1(il + 1) - 9)*xvect(mors2(mors1(il + 1) - 9)) + &
1246  xmors(mors1(il + 1) - 10)*xvect(mors2(mors1(il + 1) - 10))
1247 
1248  case (12)
1249  yvect(il) = xmors(mors1(il + 1))*xvect(mors2(mors1(il + 1))) + &
1250  xmors(mors1(il + 1) - 1)*xvect(mors2(mors1(il + 1) - 1)) + &
1251  xmors(mors1(il + 1) - 2)*xvect(mors2(mors1(il + 1) - 2)) + &
1252  xmors(mors1(il + 1) - 3)*xvect(mors2(mors1(il + 1) - 3)) + &
1253  xmors(mors1(il + 1) - 4)*xvect(mors2(mors1(il + 1) - 4)) + &
1254  xmors(mors1(il + 1) - 5)*xvect(mors2(mors1(il + 1) - 5)) + &
1255  xmors(mors1(il + 1) - 6)*xvect(mors2(mors1(il + 1) - 6)) + &
1256  xmors(mors1(il + 1) - 7)*xvect(mors2(mors1(il + 1) - 7)) + &
1257  xmors(mors1(il + 1) - 8)*xvect(mors2(mors1(il + 1) - 8)) + &
1258  xmors(mors1(il + 1) - 9)*xvect(mors2(mors1(il + 1) - 9)) + &
1259  xmors(mors1(il + 1) - 10)*xvect(mors2(mors1(il + 1) - 10)) + &
1260  xmors(mors1(il + 1) - 11)*xvect(mors2(mors1(il + 1) - 11))
1261 
1262  end select
1263 
1264  end do
1265 
1266  end subroutine m1p
1267 
1268 ! Compute Ex et Ey from potential
1269  subroutine poliss(self, phi, ex, ey)
1270 
1271  type(sll_t_poisson_2d_triangular), intent(inout) :: self
1272  sll_real64, dimension(:), intent(in) :: phi
1273  sll_real64, dimension(:), intent(out) :: ex
1274  sll_real64, dimension(:), intent(out) :: ey
1275 
1276  sll_int32 :: is, ic, nbc, iac
1277  LOGICAL :: lerr
1278 
1279  lerr = .false.
1280 
1281  do ic = 1, self%mesh%nbtcot
1282  self%vtanty(ic) = (phi(self%mesh%nuvac(1, ic)) - phi(self%mesh%nuvac(2, ic)))/self%mesh%xlcod(ic)
1283  end do
1284 
1285  do ic = 1, self%mesh%nbtcot
1286  self%vtantx(ic) = self%vtanty(ic)*self%mesh%vtaux(ic)
1287  end do
1288 
1289  do ic = 1, self%mesh%nbtcot
1290  self%vtanty(ic) = self%vtanty(ic)*self%mesh%vtauy(ic)
1291  end do
1292 
1293  do is = 1, self%mesh%num_nodes
1294 
1295  iac = self%mesh%nbcov(is) + 1
1296  nbc = self%mesh%nbcov(is + 1) - self%mesh%nbcov(is)
1297 
1298  if (nbc == 6) then
1299 
1300  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1301  + self%vtantx(self%mesh%nugcv(iac + 1)) &
1302  + self%vtantx(self%mesh%nugcv(iac + 2)) &
1303  + self%vtantx(self%mesh%nugcv(iac + 3)) &
1304  + self%vtantx(self%mesh%nugcv(iac + 4)) &
1305  + self%vtantx(self%mesh%nugcv(iac + 5))
1306 
1307  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1308  + self%vtanty(self%mesh%nugcv(iac + 1)) &
1309  + self%vtanty(self%mesh%nugcv(iac + 2)) &
1310  + self%vtanty(self%mesh%nugcv(iac + 3)) &
1311  + self%vtanty(self%mesh%nugcv(iac + 4)) &
1312  + self%vtanty(self%mesh%nugcv(iac + 5))
1313 
1314  else if (nbc == 5) then
1315 
1316  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1317  + self%vtantx(self%mesh%nugcv(iac + 1)) &
1318  + self%vtantx(self%mesh%nugcv(iac + 2)) &
1319  + self%vtantx(self%mesh%nugcv(iac + 3)) &
1320  + self%vtantx(self%mesh%nugcv(iac + 4))
1321 
1322  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1323  + self%vtanty(self%mesh%nugcv(iac + 1)) &
1324  + self%vtanty(self%mesh%nugcv(iac + 2)) &
1325  + self%vtanty(self%mesh%nugcv(iac + 3)) &
1326  + self%vtanty(self%mesh%nugcv(iac + 4))
1327 
1328  else if (nbc == 7) then
1329 
1330  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1331  + self%vtantx(self%mesh%nugcv(iac + 1)) &
1332  + self%vtantx(self%mesh%nugcv(iac + 2)) &
1333  + self%vtantx(self%mesh%nugcv(iac + 3)) &
1334  + self%vtantx(self%mesh%nugcv(iac + 4)) &
1335  + self%vtantx(self%mesh%nugcv(iac + 5)) &
1336  + self%vtantx(self%mesh%nugcv(iac + 6))
1337 
1338  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1339  + self%vtanty(self%mesh%nugcv(iac + 1)) &
1340  + self%vtanty(self%mesh%nugcv(iac + 2)) &
1341  + self%vtanty(self%mesh%nugcv(iac + 3)) &
1342  + self%vtanty(self%mesh%nugcv(iac + 4)) &
1343  + self%vtanty(self%mesh%nugcv(iac + 5)) &
1344  + self%vtanty(self%mesh%nugcv(iac + 6))
1345 
1346  else if (nbc == 4) then
1347 
1348  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1349  + self%vtantx(self%mesh%nugcv(iac + 1)) &
1350  + self%vtantx(self%mesh%nugcv(iac + 2)) &
1351  + self%vtantx(self%mesh%nugcv(iac + 3))
1352 
1353  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1354  + self%vtanty(self%mesh%nugcv(iac + 1)) &
1355  + self%vtanty(self%mesh%nugcv(iac + 2)) &
1356  + self%vtanty(self%mesh%nugcv(iac + 3))
1357 
1358  else if (nbc == 8) then
1359 
1360  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1361  + self%vtantx(self%mesh%nugcv(iac + 1)) &
1362  + self%vtantx(self%mesh%nugcv(iac + 2)) &
1363  + self%vtantx(self%mesh%nugcv(iac + 3)) &
1364  + self%vtantx(self%mesh%nugcv(iac + 4)) &
1365  + self%vtantx(self%mesh%nugcv(iac + 5)) &
1366  + self%vtantx(self%mesh%nugcv(iac + 6)) &
1367  + self%vtantx(self%mesh%nugcv(iac + 7))
1368 
1369  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1370  + self%vtanty(self%mesh%nugcv(iac + 1)) &
1371  + self%vtanty(self%mesh%nugcv(iac + 2)) &
1372  + self%vtanty(self%mesh%nugcv(iac + 3)) &
1373  + self%vtanty(self%mesh%nugcv(iac + 4)) &
1374  + self%vtanty(self%mesh%nugcv(iac + 5)) &
1375  + self%vtanty(self%mesh%nugcv(iac + 6)) &
1376  + self%vtanty(self%mesh%nugcv(iac + 7))
1377 
1378  else if (nbc == 3) then
1379 
1380  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1381  + self%vtantx(self%mesh%nugcv(iac + 1)) &
1382  + self%vtantx(self%mesh%nugcv(iac + 2))
1383 
1384  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1385  + self%vtanty(self%mesh%nugcv(iac + 1)) &
1386  + self%vtanty(self%mesh%nugcv(iac + 2))
1387 
1388  else if (nbc == 9) then
1389 
1390  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1391  + self%vtantx(self%mesh%nugcv(iac + 1)) &
1392  + self%vtantx(self%mesh%nugcv(iac + 2)) &
1393  + self%vtantx(self%mesh%nugcv(iac + 3)) &
1394  + self%vtantx(self%mesh%nugcv(iac + 4)) &
1395  + self%vtantx(self%mesh%nugcv(iac + 5)) &
1396  + self%vtantx(self%mesh%nugcv(iac + 6)) &
1397  + self%vtantx(self%mesh%nugcv(iac + 7)) &
1398  + self%vtantx(self%mesh%nugcv(iac + 8))
1399 
1400  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1401  + self%vtanty(self%mesh%nugcv(iac + 1)) &
1402  + self%vtanty(self%mesh%nugcv(iac + 2)) &
1403  + self%vtanty(self%mesh%nugcv(iac + 3)) &
1404  + self%vtanty(self%mesh%nugcv(iac + 4)) &
1405  + self%vtanty(self%mesh%nugcv(iac + 5)) &
1406  + self%vtanty(self%mesh%nugcv(iac + 6)) &
1407  + self%vtanty(self%mesh%nugcv(iac + 7)) &
1408  + self%vtanty(self%mesh%nugcv(iac + 8))
1409 
1410  else if (nbc == 2) then
1411 
1412  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1413  + self%vtantx(self%mesh%nugcv(iac + 1))
1414 
1415  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1416  + self%vtanty(self%mesh%nugcv(iac + 1))
1417 
1418  else if (nbc == 10) then
1419 
1420  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1421  + self%vtantx(self%mesh%nugcv(iac + 1)) &
1422  + self%vtantx(self%mesh%nugcv(iac + 2)) &
1423  + self%vtantx(self%mesh%nugcv(iac + 3)) &
1424  + self%vtantx(self%mesh%nugcv(iac + 4)) &
1425  + self%vtantx(self%mesh%nugcv(iac + 5)) &
1426  + self%vtantx(self%mesh%nugcv(iac + 6)) &
1427  + self%vtantx(self%mesh%nugcv(iac + 7)) &
1428  + self%vtantx(self%mesh%nugcv(iac + 8)) &
1429  + self%vtantx(self%mesh%nugcv(iac + 9))
1430 
1431  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1432  + self%vtanty(self%mesh%nugcv(iac + 1)) &
1433  + self%vtanty(self%mesh%nugcv(iac + 2)) &
1434  + self%vtanty(self%mesh%nugcv(iac + 3)) &
1435  + self%vtanty(self%mesh%nugcv(iac + 4)) &
1436  + self%vtanty(self%mesh%nugcv(iac + 5)) &
1437  + self%vtanty(self%mesh%nugcv(iac + 6)) &
1438  + self%vtanty(self%mesh%nugcv(iac + 7)) &
1439  + self%vtanty(self%mesh%nugcv(iac + 8)) &
1440  + self%vtanty(self%mesh%nugcv(iac + 9))
1441 
1442  else if (nbc == 11) then
1443 
1444  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1445  + self%vtantx(self%mesh%nugcv(iac + 1)) &
1446  + self%vtantx(self%mesh%nugcv(iac + 2)) &
1447  + self%vtantx(self%mesh%nugcv(iac + 3)) &
1448  + self%vtantx(self%mesh%nugcv(iac + 4)) &
1449  + self%vtantx(self%mesh%nugcv(iac + 5)) &
1450  + self%vtantx(self%mesh%nugcv(iac + 6)) &
1451  + self%vtantx(self%mesh%nugcv(iac + 7)) &
1452  + self%vtantx(self%mesh%nugcv(iac + 8)) &
1453  + self%vtantx(self%mesh%nugcv(iac + 9)) &
1454  + self%vtantx(self%mesh%nugcv(iac + 10))
1455 
1456  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1457  + self%vtanty(self%mesh%nugcv(iac + 1)) &
1458  + self%vtanty(self%mesh%nugcv(iac + 2)) &
1459  + self%vtanty(self%mesh%nugcv(iac + 3)) &
1460  + self%vtanty(self%mesh%nugcv(iac + 4)) &
1461  + self%vtanty(self%mesh%nugcv(iac + 5)) &
1462  + self%vtanty(self%mesh%nugcv(iac + 6)) &
1463  + self%vtanty(self%mesh%nugcv(iac + 7)) &
1464  + self%vtanty(self%mesh%nugcv(iac + 8)) &
1465  + self%vtanty(self%mesh%nugcv(iac + 9)) &
1466  + self%vtanty(self%mesh%nugcv(iac + 10))
1467 
1468  else if (nbc == 12) then
1469 
1470  self%sv1(is) = self%vtantx(self%mesh%nugcv(iac)) &
1471  + self%vtantx(self%mesh%nugcv(iac + 1)) &
1472  + self%vtantx(self%mesh%nugcv(iac + 2)) &
1473  + self%vtantx(self%mesh%nugcv(iac + 3)) &
1474  + self%vtantx(self%mesh%nugcv(iac + 4)) &
1475  + self%vtantx(self%mesh%nugcv(iac + 5)) &
1476  + self%vtantx(self%mesh%nugcv(iac + 6)) &
1477  + self%vtantx(self%mesh%nugcv(iac + 7)) &
1478  + self%vtantx(self%mesh%nugcv(iac + 8)) &
1479  + self%vtantx(self%mesh%nugcv(iac + 9)) &
1480  + self%vtantx(self%mesh%nugcv(iac + 10)) &
1481  + self%vtantx(self%mesh%nugcv(iac + 11))
1482 
1483  self%sv2(is) = self%vtanty(self%mesh%nugcv(iac)) &
1484  + self%vtanty(self%mesh%nugcv(iac + 1)) &
1485  + self%vtanty(self%mesh%nugcv(iac + 2)) &
1486  + self%vtanty(self%mesh%nugcv(iac + 3)) &
1487  + self%vtanty(self%mesh%nugcv(iac + 4)) &
1488  + self%vtanty(self%mesh%nugcv(iac + 5)) &
1489  + self%vtanty(self%mesh%nugcv(iac + 6)) &
1490  + self%vtanty(self%mesh%nugcv(iac + 7)) &
1491  + self%vtanty(self%mesh%nugcv(iac + 8)) &
1492  + self%vtanty(self%mesh%nugcv(iac + 9)) &
1493  + self%vtanty(self%mesh%nugcv(iac + 10)) &
1494  + self%vtanty(self%mesh%nugcv(iac + 11))
1495  else
1496 
1497  lerr = .true.
1498 
1499  end if
1500 
1501  end do
1502 
1503  if (lerr) then
1504  write (6, 900)
1505  stop
1506  end if
1507 
1508 ! --- 3.0 --- Solve systems 2*2 --------------------
1509 
1510  do is = 1, self%mesh%num_nodes
1511  ex(is) = self%mesh%xmal2(is)*self%sv1(is) - self%mesh%xmal3(is)*self%sv2(is)
1512  ey(is) = self%mesh%xmal1(is)*self%sv2(is) - self%mesh%xmal3(is)*self%sv1(is)
1513  end do
1514 
1515 ! --- 9.0 --- Formats --------------------------------------------------
1516 
1517 900 format(//10x, 'Error dans POLISS' &
1518  /10x, 'On a trouve more de 12 edges')
1519 
1520  end subroutine poliss
1521 
1522 end module sll_m_poisson_2d_tri
1523 
sll_m_poisson_2d_tri::sll_o_create
Definition: sll_m_poisson_2d_tri.F90:31
sll_m_choleski
Definition: sll_m_choleski.F90:1
sll_m_choleski::sll_s_choles
subroutine, public sll_s_choles(mudl, ae, as)
Definition: sll_m_choleski.F90:34
sll_m_choleski::sll_s_desrem
subroutine, public sll_s_desrem(mudl, a, be, ntdl, bs)
Definition: sll_m_choleski.F90:114
sll_m_poisson_2d_tri
Definition: sll_m_poisson_2d_tri.F90:5
sll_m_poisson_2d_tri::asblp
subroutine asblp(iprof, aele, i1, i2, i3, xmass)
Definition: sll_m_poisson_2d_tri.F90:1094
sll_m_poisson_2d_tri::read_data_solver
subroutine read_data_solver(ntypfr, potfr)
Definition: sll_m_poisson_2d_tri.F90:186
sll_m_poisson_2d_tri::morse
subroutine morse(npoel1, npoel2, ntri, nbt, nbs, mors1, mors2)
Definition: sll_m_poisson_2d_tri.F90:446
sll_m_poisson_2d_tri::initialize_poisson_solver
subroutine initialize_poisson_solver(self, mesh, ntypfr, potfr)
Definition: sll_m_poisson_2d_tri.F90:307
sll_m_poisson_2d_tri::poifrc
subroutine poifrc(self, ex, ey)
Definition: sll_m_poisson_2d_tri.F90:774
sll_m_poisson_2d_tri::m1p
subroutine m1p(xmors, mors1, mors2, xvect, nlign, yvect)
Definition: sll_m_poisson_2d_tri.F90:1157
sll_m_poisson_2d_tri::poismc
subroutine poismc(self)
Definition: sll_m_poisson_2d_tri.F90:571
sll_m_poisson_2d_tri::poissn
subroutine poissn(self, rho, phi, ex, ey)
Definition: sll_m_poisson_2d_tri.F90:696
sll_m_poisson_2d_tri::initialize_poisson_solver_from_file
subroutine initialize_poisson_solver_from_file(self, mesh)
Definition: sll_m_poisson_2d_tri.F90:278
sll_m_poisson_2d_tri::sll_s_poisson_2d_triangular_free
subroutine, public sll_s_poisson_2d_triangular_free(self)
Delete the solver derived type.
Definition: sll_m_poisson_2d_tri.F90:139
sll_m_poisson_2d_tri::asblm2
subroutine asblm2(aele1, aele2, mors1, mors2, i1, i2, i3, a1, a2)
Definition: sll_m_poisson_2d_tri.F90:999
sll_m_poisson_2d_tri::sll_s_poisson_2d_triangular_init
subroutine, public sll_s_poisson_2d_triangular_init(solver, mesh, ntypfr, potfr)
Definition: sll_m_poisson_2d_tri.F90:127
sll_m_poisson_2d_tri::sll_s_compute_e_from_phi
subroutine, public sll_s_compute_e_from_phi(self, phi, ex, ey)
Definition: sll_m_poisson_2d_tri.F90:175
sll_m_poisson_2d_tri::profil
subroutine profil(ntri, nbs, npoel1, npoel2, iprof)
Definition: sll_m_poisson_2d_tri.F90:893
sll_m_poisson_2d_tri::sll_s_compute_e_from_rho
subroutine, public sll_s_compute_e_from_rho(self, rho, phi, ex, ey)
Compute electric field from charge density.
Definition: sll_m_poisson_2d_tri.F90:154
sll_m_poisson_2d_tri::grandx
real(kind=f64), parameter grandx
Definition: sll_m_poisson_2d_tri.F90:89
sll_m_poisson_2d_tri::asbld
subroutine asbld(aele, i1, i2, i3, xmass)
Definition: sll_m_poisson_2d_tri.F90:968
sll_m_poisson_2d_tri::poliss
subroutine poliss(self, phi, ex, ey)
Definition: sll_m_poisson_2d_tri.F90:1270
sll_m_poisson_2d_tri::sll_s_compute_phi_from_rho
subroutine, public sll_s_compute_phi_from_rho(self, rho, phi)
Definition: sll_m_poisson_2d_tri.F90:166
sll_m_triangular_meshes
Definition: sll_m_triangular_meshes.F90:19
sll_m_poisson_2d_tri::sll_t_poisson_2d_triangular
Derived type for Poisson solver on unstructured mesh with triangles.
Definition: sll_m_poisson_2d_tri.F90:96
sll_m_triangular_meshes::sll_t_triangular_mesh_2d
2d hexagonal mesh
Definition: sll_m_triangular_meshes.F90:58
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1