Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_advection_2d_tri_mesh.F90
Go to the documentation of this file.
1 
3 module sll_m_advection_2d_tri_mesh
4 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5 #include "sll_memory.h"
6 #include "sll_working_precision.h"
7 
8  use sll_m_triangular_meshes, only: &
9  sll_t_triangular_mesh_2d
10 
11  implicit none
12 
13  public :: &
14  sll_s_advection_2d, &
15  sll_f_new_advection_2d_tri_mesh, &
16  sll_t_advection_tri_mesh
17 
18  private
19 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
20 
22  type :: sll_t_advection_tri_mesh
23 
24  type(sll_t_triangular_mesh_2d), pointer :: mesh
25  sll_int32, dimension(:, :), pointer :: nvoiv
26  sll_int32, dimension(:), pointer :: nlpa
27  sll_int32, dimension(:), pointer :: nbpama
28  sll_int32, dimension(:), pointer :: iad1
29  sll_int32, dimension(:), pointer :: indice
30  sll_int32, dimension(:), pointer :: itabor
31  sll_int32, dimension(:), pointer :: itest
32  sll_int32, dimension(:), pointer :: nlmloc
33  sll_real64, dimension(:), pointer :: xbas
34  sll_real64, dimension(:, :), pointer :: xlm
35  sll_real64, dimension(:, :), pointer :: coef
36  sll_real64, dimension(:), pointer :: xp
37  sll_real64, dimension(:), pointer :: yp
38  sll_real64, dimension(:), pointer :: f_out
39  logical, dimension(:), pointer :: inzone
40  sll_int32, dimension(:), pointer :: numres
41 
42  end type sll_t_advection_tri_mesh
43 
44 ! type :: sll_degree_of_freedom
45 ! sll_real64, dimension(:), pointer :: func !> values of the distribution function
46 ! sll_real64, dimension(:), pointer :: deriv !> values of the derivatives of the distribution
47 ! sll_real64, dimension(:), pointer :: deriv2 !> values of the 2nd derivatives of the distribution
48 
49 ! sll_real64 :: epsilon
50 ! end type sll_degree_of_freedom
51 
52 contains
53 
60  function sll_f_new_advection_2d_tri_mesh(mesh) result(adv)
61 
62  type(sll_t_triangular_mesh_2d), intent(in), target :: mesh
63  type(sll_t_advection_tri_mesh), pointer :: adv
64 
65  sll_int32 :: ierr
66  sll_int32 :: i
67  sll_int32 :: j
68  sll_int32 :: is1, is2, is3
69 
70  sll_allocate(adv, ierr)
71 
72  adv%mesh => mesh
73 
74  ! --- Remplissage de "nvoiv" -----------------------------------
75  ! Identique a "nvois" sauf pour les aretes appartenant a
76  ! une frontiere. Le chiffre correspond ici a un code pour
77  ! le traitement des conditions aux limites sur les
78  ! particules, alors que dans "nvois" ce chiffre est
79  ! l'oppose du numero de reference de la frontiere concernee
80 
81  allocate (adv%nvoiv(3, mesh%num_triangles))
82  do i = 1, mesh%num_triangles
83  do j = 1, 3
84  if (mesh%nvois(j, i) > 0) then
85  adv%nvoiv(j, i) = mesh%nvois(j, i)
86  else
87  adv%nvoiv(j, i) = 0
88  end if
89  end do
90  end do
91 
92  allocate (adv%nlpa(mesh%num_nodes))
93  do i = 1, mesh%num_nodes
94  adv%nlpa(i) = mesh%npoel2(mesh%npoel1(i) + 1)
95  end do
96 
97  sll_allocate(adv%itest(mesh%num_nodes), ierr)
98  sll_allocate(adv%coef(4, mesh%num_nodes), ierr)
99  sll_allocate(adv%xlm(3, mesh%num_nodes), ierr)
100  sll_allocate(adv%nlmloc(mesh%num_nodes), ierr)
101  sll_allocate(adv%xp(mesh%num_nodes), ierr)
102  sll_allocate(adv%yp(mesh%num_nodes), ierr)
103  sll_allocate(adv%inzone(mesh%num_nodes), ierr)
104  allocate (adv%numres(mesh%num_nodes))
105  allocate (adv%iad1(mesh%num_triangles))
106  allocate (adv%indice(mesh%num_triangles))
107  allocate (adv%itabor(mesh%num_nodes))
108  allocate (adv%f_out(mesh%num_nodes))
109  allocate (adv%nbpama(mesh%num_triangles))
110  allocate (adv%xbas(mesh%num_nodes))
111 
112  do i = 1, mesh%num_triangles
113 
114  is1 = mesh%nodes(1, i)
115  is2 = mesh%nodes(2, i)
116  is3 = mesh%nodes(3, i)
117 
118  adv%xbas(is1) = adv%xbas(is1) + mesh%area(i)/3.
119  adv%xbas(is2) = adv%xbas(is2) + mesh%area(i)/3.
120  adv%xbas(is3) = adv%xbas(is3) + mesh%area(i)/3.
121  end do
122 
123  end function sll_f_new_advection_2d_tri_mesh
124 
138  subroutine compute_derivatives(f_val, f_der, f_der2, degree, epsilon)
139  sll_real64, dimension(:), intent(in) :: f_val
140  sll_real64, dimension(:), intent(out) :: f_der
141  sll_real64, dimension(:), intent(out) :: f_der2
142  sll_int32, intent(in) :: degree
143  sll_real64, intent(in) :: epsilon
144  sll_int32 :: i
145  sll_int32 :: i1
146  sll_int32 :: i2
147  sll_int32 :: i3
148 
149  i1 = 2
150  i2 = 3
151  i3 = 4
152  do i = 1, degree
153  f_der(i) = (f_val(i1) - f_val(1))/epsilon
154  f_der2(i) = (f_val(i2) - f_val(i1) - f_val(i3) + f_val(1))/(epsilon*epsilon)
155  i1 = i3
156  i2 = i2 + 2
157  i3 = modulo(i3 + 1, 2*degree) + 1
158  end do
159  end subroutine compute_derivatives
160 
174  subroutine compute_coordinates_dof(x1_adv, x2_adv, x1_coo, x2_coo, degree, epsilon, x1_dof, x2_dof)
175  sll_real64, intent(in) :: x1_adv
176  sll_real64, intent(in) :: x2_adv
177  sll_real64, dimension(:), intent(in) :: x1_coo
178  sll_real64, dimension(:), intent(in) :: x2_coo
179  sll_int32, intent(in) :: degree
180  sll_real64, intent(in) :: epsilon
181  sll_real64, dimension(:), intent(out) :: x1_dof
182  sll_real64, dimension(:), intent(out) :: x2_dof
183  sll_int32 :: i
184  sll_int32 :: i1
185  sll_int32 :: i2
186 
187  x1_dof(1) = x1_adv
188  x2_dof(1) = x2_adv
189  i1 = 2
190  i2 = 3
191 
192  do i = 1, degree
193  x1_dof(2*i) = x1_dof(1) + epsilon*(x1_coo(i1) - x1_coo(1))
194  x2_dof(2*i) = x2_dof(1) + epsilon*(x2_coo(i1) - x2_coo(1))
195 
196  x1_dof(2*i + 1) = x1_dof(1) + epsilon*(x1_coo(i1) - x1_coo(1)*2.0_f64 + x1_coo(i2))
197  x2_dof(2*i + 1) = x2_dof(1) + epsilon*(x2_coo(i1) - x2_coo(1)*2.0_f64 + x2_coo(i2))
198 
199  i1 = i2
200  i2 = modulo(i2 - 1, degree) + 2
201  end do
202  end subroutine compute_coordinates_dof
203 
212  function eval_at_lambda(lam, func_loc, der_loc, der2_loc) result(out)
213  sll_real64, dimension(3), intent(in) :: lam
214  sll_real64, dimension(3), intent(in) :: func_loc
215  sll_real64, dimension(6), intent(in) :: der_loc
216  sll_real64, dimension(3), intent(in) :: der2_loc
217  sll_real64 :: out
218 
223  out = lam(1)**2*(3._f64 - 2._f64*lam(1) + 6._f64*lam(2)*lam(3))*func_loc(1)
224  out = out + lam(2)**2*(3._f64 - 2._f64*lam(2) + 6._f64*lam(3)*lam(1))*func_loc(2)
225  out = out + lam(3)**2*(3._f64 - 2._f64*lam(3) + 6._f64*lam(1)*lam(2))*func_loc(3)
226 
227  !Sum over first derivatives
231  out = out + lam(1)**2*lam(2)*(1 + 2*lam(3))*der_loc(1) !df12
232  out = out + lam(1)**2*lam(3)*(1 + 2*lam(2))*der_loc(2) !df13
233  out = out + lam(2)**2*lam(3)*(1 + 2*lam(1))*der_loc(3) !df23
234  out = out + lam(2)**2*lam(1)*(1 + 2*lam(3))*der_loc(4) !df21
235  out = out + lam(3)**2*lam(1)*(1 + 2*lam(2))*der_loc(5) !df31
236  out = out + lam(3)**2*lam(2)*(1 + 2*lam(1))*der_loc(6) !df32
237 
242  out = out + lam(1)**2*lam(2)*lam(3)*der2_loc(1)
243  out = out + lam(2)**2*lam(3)*lam(1)*der2_loc(2)
244  out = out + lam(3)**2*lam(1)*lam(2)*der2_loc(3)
245 
246  end function eval_at_lambda
247 
248  subroutine interpolation_mitchell()
249 
250  end subroutine
306  subroutine sll_s_advection_2d(adv, f_in, ex, ey, dt)
307 
308  type(sll_t_advection_tri_mesh), intent(inout) :: adv
309  sll_real64, dimension(:), intent(inout) :: f_in
310  sll_real64, dimension(:), intent(in) :: ex
311  sll_real64, dimension(:), intent(in) :: ey
312  sll_real64, intent(in) :: dt
313 
314  sll_real64 :: eps
315  sll_int32 :: nbpert
316  sll_int32 :: nbpres
317  sll_int32 :: num
318  sll_int32 :: nrest
319  sll_int32 :: nfin
320  sll_int32 :: nbpr
321 
322  sll_int32 :: ip
323  sll_int32 :: jp
324 
325  sll_real64 :: pa1x, pa1y, pa2x, pa2y, pa3x, pa3y
326 
327  sll_real64 :: phi1, phi2, phi3
328  sll_real64 :: xm11, xm12, xm13
329  sll_int32 :: nprest
330  sll_int32 :: mpa, inum, ip1, ip2, ks, ind
331  sll_int32 :: it, is1, is2, is3
332  sll_real64 :: x1, x2, x3, y1, y2, y3, det
333 
334  eps = -adv%mesh%petitl**2
335 
336 !Set the position of the characteristic origin
337 !Cell is arbitrary set in one of the cell close to the ip vertex.
338  do ip = 1, adv%mesh%num_nodes
339  adv%numres(ip) = ip
340  adv%xp(ip) = adv%mesh%coord(1, ip) + ex(ip)*dt
341  adv%yp(ip) = adv%mesh%coord(2, ip) + ey(ip)*dt
342  adv%nlmloc(ip) = adv%nlpa(ip)
343  end do
344 
345  nbpres = adv%mesh%num_nodes
346  nbpert = 0
347  num = 0
348 
349  adv%inzone = .false.
350 
351  do while (nbpres > 0)
352 
353  !*** Premiere boucle dans le meme element
354 
355  nfin = 0
356  nrest = 0
357  num = num + 1
358 
359  if (num == 100) then
360 
361  write (*, "(//2x,'Arret dans POSITIONS:')", advance='no')
362  write (*, "('on n''arrive pas a positionner ')", advance='no')
363  write (*, "(i4)", advance='no') nbpres
364  write (*, "(' particules')", advance='no')
365  write (*, "(/2x,'Particules non positionnees :')", advance='no')
366  write (*, "(5x,'numero - position - vitesse - element'/)")
367  do ip = 1, min(50, nbpres)
368  write (*, "(i10,2x,4e12.4,i10)") adv%numres(ip), &
369  adv%xp(ip), adv%yp(ip), &
370  ex(ip), ey(ip), ip
371  end do
372  write (*, "(/5x,a)") &
373  'Il y a certainement un probleme avec le solveur de champ'
374  write (*, "(/5x,'(en general la CFL n''est pas verifiee)'/)")
375 
376  stop
377 
378  end if
379 
380  do ip = 1, nbpres
381 
382  jp = adv%numres(ip)
383  it = adv%nlmloc(ip)
384 
385  x1 = adv%mesh%coord(1, adv%mesh%nodes(1, it))
386  y1 = adv%mesh%coord(2, adv%mesh%nodes(1, it))
387  x2 = adv%mesh%coord(1, adv%mesh%nodes(2, it))
388  y2 = adv%mesh%coord(2, adv%mesh%nodes(2, it))
389  x3 = adv%mesh%coord(1, adv%mesh%nodes(3, it))
390  y3 = adv%mesh%coord(2, adv%mesh%nodes(3, it))
391 
392  pa1x = x1 - adv%xp(jp)
393  pa1y = y1 - adv%yp(jp)
394  pa2x = x2 - adv%xp(jp)
395  pa2y = y2 - adv%yp(jp)
396  pa3x = x3 - adv%xp(jp)
397  pa3y = y3 - adv%yp(jp)
398 
399  adv%coef(1, ip) = pa1x*pa2y - pa1y*pa2x
400  adv%coef(2, ip) = pa2x*pa3y - pa2y*pa3x
401  adv%coef(3, ip) = pa3x*pa1y - pa3y*pa1x
402 
403  if (adv%coef(1, ip) >= eps &
404  .and. adv%coef(2, ip) >= eps &
405  .and. adv%coef(3, ip) >= eps) then
406 
407  nfin = nfin + 1
408 
409  det = (x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1)
410 
411  adv%xlm(1, jp) = adv%coef(1, ip)/det
412  adv%xlm(2, jp) = adv%coef(2, ip)/det
413  adv%xlm(3, jp) = adv%coef(3, ip)/det
414 
415  adv%nlpa(jp) = adv%nlmloc(ip)
416  adv%itest(ip) = 0
417  adv%inzone(jp) = .true.
418 
419  end if
420 
421  end do
422 
423  !*** Deuxieme boucle pour celles qui sont sorties
424 
425  nbpr = nbpres - nfin
426 
427  if (nbpr .ne. 0) then
428 
429  do ip = 1, nbpres
430 
431  jp = adv%numres(ip)
432  it = adv%nlmloc(ip)
433 
434  if (adv%coef(1, ip) < eps &
435  .and. adv%coef(2, ip) >= eps &
436  .and. adv%coef(3, ip) >= eps) then
437 
438  adv%itest(ip) = 1 + 10*(1 - min(1, adv%nvoiv(1, it)))
439 
440  end if
441 
442  if (adv%coef(1, ip) >= eps &
443  .and. adv%coef(2, ip) < eps &
444  .and. adv%coef(3, ip) >= eps) then
445 
446  adv%itest(ip) = 2 + 10*(1 - min(1, adv%nvoiv(2, it)))
447 
448  end if
449 
450  if (adv%coef(1, ip) >= eps &
451  .and. adv%coef(2, ip) >= eps &
452  .and. adv%coef(3, ip) < eps) then
453 
454  adv%itest(ip) = 3 + 10*(1 - min(1, adv%nvoiv(3, it)))
455 
456  end if
457 
458  if (adv%coef(1, ip) < eps &
459  .and. adv%coef(2, ip) < eps) then
460 
461  pa2x = adv%mesh%coord(1, adv%mesh%nodes(2, it)) - adv%xp(jp)
462  pa2y = adv%mesh%coord(2, adv%mesh%nodes(2, it)) - adv%yp(jp)
463 
464  adv%coef(4, ip) = pa2x*ey(jp) - pa2y*ex(jp)
465 
466  adv%itest(ip) = 1 + max(0, nint(sign(1.0_f64, adv%coef(4, ip))))
467  adv%itest(ip) = adv%itest(ip) + 10*(1 - min(1, adv%nvoiv(adv%itest(ip), it)))
468 
469  end if
470 
471  if (adv%coef(2, ip) < eps &
472  .and. adv%coef(3, ip) < eps) then
473 
474  pa3x = adv%mesh%coord(1, adv%mesh%nodes(3, it)) - adv%xp(jp)
475  pa3y = adv%mesh%coord(2, adv%mesh%nodes(3, it)) - adv%yp(jp)
476 
477  adv%coef(4, ip) = pa3x*ey(jp) - pa3y*ex(jp)
478 
479  adv%itest(ip) = 2 + max(0, nint(sign(1.0_f64, adv%coef(4, ip))))
480  adv%itest(ip) = adv%itest(ip) + 10*(1 - min(1, adv%nvoiv(adv%itest(ip), it)))
481 
482  end if
483 
484  if (adv%coef(3, ip) < eps &
485  .and. adv%coef(1, ip) < eps) then
486 
487  pa1x = adv%mesh%coord(1, adv%mesh%nodes(1, it)) - adv%xp(jp)
488  pa1y = adv%mesh%coord(2, adv%mesh%nodes(1, it)) - adv%yp(jp)
489 
490  adv%coef(4, ip) = pa1x*ey(jp) - pa1y*ex(jp)
491 
492  adv%itest(ip) = 1 + mod(2 + max(0, nint(sign(1.0_f64, adv%coef(4, ip)))), 3)
493  adv%itest(ip) = adv%itest(ip) + 10*(1 - min(1, adv%nvoiv(adv%itest(ip), it)))
494 
495  end if
496 
497  end do
498 
499  !*** Particules absorbees par une frontiere -------------------
500 
501  do ip = 1, nbpres
502 
503  if (adv%itest(ip) > 10 .and. adv%itest(ip) < 14) then
504  nbpert = nbpert + 1
505  end if
506 
507  end do
508 
509  !*** Reorganisation des tableaux temporaires ------------------
510  ! Particules traversant un cote interne ou reflechie
511 
512  do ip = 1, nbpres
513 
514  if ((adv%itest(ip) > 0 .and. adv%itest(ip) < 4) .or. &
515  (adv%itest(ip) > 20 .and. adv%itest(ip) < 24)) then
516  nrest = nrest + 1
517  adv%numres(nrest) = adv%numres(ip)
518 
519  if (adv%itest(ip) > 20) then
520  adv%nlmloc(nrest) = adv%mesh%nvois(adv%itest(ip) - 20, adv%nlmloc(ip)) &
521  *max(0, sign(1, 10 - adv%itest(ip))) &
522  + adv%nlmloc(ip)*max(0, sign(1, adv%itest(ip) - 20))
523  else
524  adv%nlmloc(nrest) = adv%mesh%nvois(adv%itest(ip), adv%nlmloc(ip)) &
525  *max(0, sign(1, 10 - adv%itest(ip))) &
526  + adv%nlmloc(ip)*max(0, sign(1, adv%itest(ip) - 20))
527  end if
528 
529  end if
530 
531  end do
532 
533  end if
534 
535  nbpres = nrest
536 
537  end do
538 
539 !Recherche du nombre de particules de chaque maille -------
540 
541  adv%nbpama = 0
542  do ip = 1, adv%mesh%num_nodes
543  if (adv%inzone(ip)) then
544  mpa = adv%nlpa(ip)
545  adv%nbpama(mpa) = adv%nbpama(mpa) + 1
546  end if
547  end do
548 
549 !--- Determination de l'adresse de la premiere particule
550 ! de chaque maille dans le tableau ordonne
551 
552  adv%iad1 = 0
553  ks = 1
554  do it = 1, adv%mesh%num_triangles
555  if (adv%nbpama(it) .ne. 0) then
556  adv%iad1(it) = ks
557  ks = ks + adv%nbpama(it)
558  end if
559  end do
560 
561 !--- Construction du tableau ordonne des particules -----------
562 
563  adv%indice = 0
564  adv%itabor = 0
565 
566  do ip = 1, adv%mesh%num_nodes
567 
568  if (adv%inzone(ip)) then
569  mpa = adv%nlpa(ip)
570  ind = adv%iad1(mpa) + adv%indice(mpa)
571  adv%itabor(ind) = ip
572  adv%indice(mpa) = adv%indice(mpa) + 1
573  end if
574 
575  end do
576 
577  nprest = adv%mesh%num_nodes
578 
579  adv%f_out = 0.0_f64
580 
581  do it = 1, adv%mesh%num_triangles
582 
583  xm11 = 0.0_f64; xm12 = 0.0_f64; xm13 = 0.0_f64
584 
585  !nbpama(it) !Nombre de particules dans la maille numero it
586 
587  if (adv%nbpama(it) .ne. 0) then
588 
589  ip1 = adv%iad1(it)
590  ip2 = ip1 + adv%nbpama(it) - 1
591 
592  !Boucle sur les particules situees dans la maille courante
593 
594  do ip = ip1, ip2
595 
596  inum = adv%itabor(ip)
597 
598  phi1 = adv%xlm(2, inum)*f_in(inum)
599  phi2 = adv%xlm(3, inum)*f_in(inum)
600  phi3 = adv%xlm(1, inum)*f_in(inum)
601 
602  xm11 = xm11 + phi1
603  xm12 = xm12 + phi2
604  xm13 = xm13 + phi3
605 
606  end do
607 
608  is1 = adv%mesh%nodes(1, it)
609  is2 = adv%mesh%nodes(2, it)
610  is3 = adv%mesh%nodes(3, it)
611 
612  adv%f_out(is1) = adv%f_out(is1) + xm11
613  adv%f_out(is2) = adv%f_out(is2) + xm12
614  adv%f_out(is3) = adv%f_out(is3) + xm13
615 
616  if (adv%nbpama(it) == nprest) exit
617 
618  nprest = nprest - adv%nbpama(it)
619 
620  end if
621 
622  end do
623 
624  f_in = adv%f_out
625 
626  end subroutine sll_s_advection_2d
627 
628 end module sll_m_advection_2d_tri_mesh
sll_m_advection_2d_tri_mesh
advection on triangular mesh.
Definition: sll_m_advection_2d_tri_mesh.F90:3
sll_m_advection_2d_tri_mesh::compute_derivatives
subroutine compute_derivatives(f_val, f_der, f_der2, degree, epsilon)
Computes degrees of freedom on one point.
Definition: sll_m_advection_2d_tri_mesh.F90:139
sll_m_advection_2d_tri_mesh::sll_f_new_advection_2d_tri_mesh
type(sll_t_advection_tri_mesh) function, pointer, public sll_f_new_advection_2d_tri_mesh(mesh)
allocates the memory space for a new 2D advection on triangular mesh on the heap, initializes it with...
Definition: sll_m_advection_2d_tri_mesh.F90:61
sll_m_advection_2d_tri_mesh::sll_s_advection_2d
subroutine, public sll_s_advection_2d(adv, f_in, ex, ey, dt)
Compute characterisitic origin in triangular mesh.
Definition: sll_m_advection_2d_tri_mesh.F90:307
sll_m_advection_2d_tri_mesh::interpolation_mitchell
subroutine interpolation_mitchell()
Definition: sll_m_advection_2d_tri_mesh.F90:249
sll_m_advection_2d_tri_mesh::eval_at_lambda
real(kind=f64) function eval_at_lambda(lam, func_loc, der_loc, der2_loc)
evaluation at barycentric points of one cell
Definition: sll_m_advection_2d_tri_mesh.F90:213
sll_m_advection_2d_tri_mesh::compute_coordinates_dof
subroutine compute_coordinates_dof(x1_adv, x2_adv, x1_coo, x2_coo, degree, epsilon, x1_dof, x2_dof)
computes the coordinates of the poisition of the degrees of freedom (dof)
Definition: sll_m_advection_2d_tri_mesh.F90:175
sll_m_triangular_meshes
Definition: sll_m_triangular_meshes.F90:19
sll_m_advection_2d_tri_mesh::sll_t_advection_tri_mesh
2d advection on triangular mesh
Definition: sll_m_advection_2d_tri_mesh.F90:22
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:41 for Semi-Lagrangian Library by   doxygen 1.9.1