Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_hex_poisson.F90
Go to the documentation of this file.
1 module sll_m_hex_poisson
2 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 #include "sll_memory.h"
4 #include "sll_working_precision.h"
5 
6  use sll_m_hexagonal_meshes, only: &
7  sll_f_new_hex_mesh_2d, &
8  sll_t_hex_mesh_2d
9 
10  implicit none
11 
12  public :: &
13  sll_s_compute_hex_fields, &
14  sll_s_hex_matrix_poisson, &
15  sll_s_hex_second_terme_poisson
16 
17  private
18 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
19 
20  !***********************************************************
21  ! self module is used to solve the Poison equation on a hexagonal mesh
22  ! by creating the poisson matrix for the hex mesh : matrix_poisson
23  ! and the second term of the linear equation : second_term
24  ! therefore we solve with a linear solver : matrix_poisson * X = second_term
25  ! we also compute the fields and it derivatives from the potential we get from
26  ! the Poisson equation
27  ! Precision is of order 4 ( for both poisson solver and fields computation )
28  ! for any question : prouveur@math.univ-lyon1.fr
29  !***********************************************************
30 
31 contains
32 
33  subroutine sll_s_hex_matrix_poisson(matrix_poisson, mesh, type)
34  type(sll_t_hex_mesh_2d), pointer :: mesh
35  sll_real64, dimension(:, :), intent(out):: matrix_poisson
36 
37  sll_int32, intent(in):: type ! unused parameter atm
38 
39  sll_int32 :: num_cells ! number of hex cells
40  sll_int32 :: global ! global index
41 
42  ! index on the matrix
43  sll_int32 :: index_tab, index_tabij
44  sll_int32 :: index_tabi_1j, index_tabij_1
45  sll_int32 :: index_tabij1, index_tabi1j
46  sll_int32 :: index_tabi1j1, index_tabi_1j_1
47  sll_int32 :: k1, k2, n
48 
49  num_cells = mesh%num_cells
50 
51  matrix_poisson = 0._f64
52 
53  ! indexation of the matrix : we fill row by row with
54  ! j + (i-1)*n(j) the index given by index_tab
55  ! therefore we get (i-1,j-1), in hex coordinate, then (i-1,j)
56  ! (i,j-1) ; (i,j) ; (i,j+1)
57  ! (i+1,j) ; (i,j+1)
58 
59  if (type == 1) then
60  n = mesh%num_pts_tot
61  elseif (type == 2) then
62  n = mesh%num_triangles
63  elseif (type == 3) then
64  n = mesh%num_edges
65  end if
66 
67  do global = 1, n
68 
69  k1 = mesh%hex_coord(1, global)
70  k2 = mesh%hex_coord(2, global)
71 
72  call mesh%index_hex_to_global(k1, k2, index_tab)
73 
74  call mesh%index_hex_to_global(k1 - 1, k2 - 1, index_tabi_1j_1)
75  call mesh%index_hex_to_global(k1 - 1, k2, index_tabi_1j)
76  call mesh%index_hex_to_global(k1, k2 - 1, index_tabij_1)
77  call mesh%index_hex_to_global(k1, k2 + 1, index_tabij1)
78  call mesh%index_hex_to_global(k1 + 1, k2, index_tabi1j)
79  call mesh%index_hex_to_global(k1 + 1, k2 + 1, index_tabi1j1)
80 
81  index_tabij = index_tab
82 
83  index_tabi_1j_1 = index_tabi_1j_1 - index_tab + 2*num_cells + 2
84  index_tabi_1j = index_tabi_1j - index_tab + 2*num_cells + 2
85  index_tabij_1 = index_tabij_1 - index_tab + 2*num_cells + 2
86  index_tabij = index_tabij - index_tab + 2*num_cells + 2
87  index_tabij1 = index_tabij1 - index_tab + 2*num_cells + 2
88  index_tabi1j = index_tabi1j - index_tab + 2*num_cells + 2
89  index_tabi1j1 = index_tabi1j1 - index_tab + 2*num_cells + 2
90 
91  if (abs(k1) == num_cells .and. abs(k2) == num_cells .or. &
92  abs(k1) == num_cells .and. k2 == 0 .or. &
93  abs(k2) == num_cells .and. k1 == 0 &
94  ) then ! corners
95 
96  if (k1 == num_cells .and. k2 == 0) then! corner top right
97  matrix_poisson(index_tab, index_tabij) = 1._f64
98  elseif (k1 == num_cells .and. k2 == num_cells) then! corner top
99  matrix_poisson(index_tab, index_tabij) = 1._f64
100  else if (k2 == num_cells .and. k1 == 0) then! corner top left
101  matrix_poisson(index_tab, index_tabij) = 1._f64
102  else if (k1 == -num_cells .and. k2 == 0) then! corner bottom left
103  matrix_poisson(index_tab, index_tabij) = 1._f64
104  else if (-k1 == num_cells .and. k2 == -num_cells) then!corner bottom
105  matrix_poisson(index_tab, index_tabij) = 1._f64
106  else if (k2 == -num_cells .and. k1 == 0) then! corner bottom right
107  matrix_poisson(index_tab, index_tabij) = 1._f64
108  end if
109 
110  else if (k1*k2 < 0 .and. abs(k1) + abs(k2) == num_cells .or. &
111  abs(k1) == num_cells .or. abs(k2) == num_cells &
112  ) then ! edges
113 
114  if (k1 == num_cells) then! edge top right
115  matrix_poisson(index_tab, index_tabij) = 1._f64
116  elseif (k2 == num_cells) then! edge top left
117  matrix_poisson(index_tab, index_tabij) = 1._f64
118  else if (k1*k2 < 0 .and. -k1 + k2 == num_cells) then! edge left
119  matrix_poisson(index_tab, index_tabij) = 1._f64
120  else if (k1 == -num_cells) then! edge bottom left
121  matrix_poisson(index_tab, index_tabij) = 1._f64
122  else if (k2 == -num_cells) then! edge bottom right
123  matrix_poisson(index_tab, index_tabij) = 1._f64
124  else if (k1*k2 < 0 .and. k1 - k2 == num_cells) then! edge right
125  matrix_poisson(index_tab, index_tabij) = 1._f64
126  end if
127 
128  elseif (abs(k1) == num_cells - 1 .and. abs(k2) == num_cells - 1 .or. &
129  abs(k1) == num_cells - 1 .and. k2 == 0 .or. &
130  abs(k2) == num_cells - 1 .and. k1 == 0 &
131  ) then ! corners
132 
133  if (k1 == num_cells - 1 .and. k2 == 0) then! corner top right
134  matrix_poisson(index_tab, index_tabi_1j_1) = -1._f64
135  matrix_poisson(index_tab, index_tabi_1j) = -1._f64
136  matrix_poisson(index_tab, index_tabij) = 6._f64
137  matrix_poisson(index_tab, index_tabij1) = -1._f64
138  elseif (k1 == num_cells - 1 .and. k2 == num_cells - 1) then! corner top
139  matrix_poisson(index_tab, index_tabi_1j_1) = -1._f64
140  matrix_poisson(index_tab, index_tabi_1j) = -1._f64
141  matrix_poisson(index_tab, index_tabij_1) = -1._f64
142  matrix_poisson(index_tab, index_tabij) = 6._f64
143  else if (k2 == num_cells - 1 .and. k1 == 0) then! corner top left
144  matrix_poisson(index_tab, index_tabi_1j_1) = -1._f64
145  matrix_poisson(index_tab, index_tabij_1) = -1._f64
146  matrix_poisson(index_tab, index_tabij) = 6._f64
147  matrix_poisson(index_tab, index_tabi1j) = -1._f64
148  else if (k1 == -num_cells + 1 .and. k2 == 0) then! corner bottom left
149  matrix_poisson(index_tab, index_tabij_1) = -1._f64
150  matrix_poisson(index_tab, index_tabij) = 6._f64
151  matrix_poisson(index_tab, index_tabi1j) = -1._f64
152  matrix_poisson(index_tab, index_tabi1j1) = -1._f64
153  else if (-k1 == num_cells - 1 .and. k2 == -num_cells + 1) then!corner bottom
154  matrix_poisson(index_tab, index_tabij) = 6._f64
155  matrix_poisson(index_tab, index_tabij1) = -1._f64
156  matrix_poisson(index_tab, index_tabi1j) = -1._f64
157  matrix_poisson(index_tab, index_tabi1j1) = -1._f64
158  else if (k2 == -num_cells + 1 .and. k1 == 0) then! corner bottom right
159  matrix_poisson(index_tab, index_tabi_1j) = -1._f64
160  matrix_poisson(index_tab, index_tabij) = 6._f64
161  matrix_poisson(index_tab, index_tabij1) = -1._f64
162  matrix_poisson(index_tab, index_tabi1j1) = -1._f64
163  end if
164 
165  else if (k1*k2 < 0 .and. abs(k1) + abs(k2) == num_cells - 1 .or. &
166  abs(k1) == num_cells - 1 .or. abs(k2) == num_cells - 1 &
167  ) then ! edges
168 
169  if (k1 == num_cells - 1) then! edge top right
170  matrix_poisson(index_tab, index_tabi_1j_1) = -1._f64
171  matrix_poisson(index_tab, index_tabi_1j) = -1._f64
172  matrix_poisson(index_tab, index_tabij_1) = -1._f64
173  matrix_poisson(index_tab, index_tabij) = 6._f64
174  matrix_poisson(index_tab, index_tabij1) = -1._f64
175  elseif (k2 == num_cells - 1) then! edge top left
176  matrix_poisson(index_tab, index_tabi_1j_1) = -1._f64
177  matrix_poisson(index_tab, index_tabi_1j) = -1._f64
178  matrix_poisson(index_tab, index_tabij_1) = -1._f64
179  matrix_poisson(index_tab, index_tabij) = 6._f64
180  matrix_poisson(index_tab, index_tabi1j) = -1._f64
181  else if (k1*k2 < 0 .and. -k1 + k2 == num_cells - 1) then! edge left
182  matrix_poisson(index_tab, index_tabi_1j_1) = -1._f64
183  matrix_poisson(index_tab, index_tabij_1) = -1._f64
184  matrix_poisson(index_tab, index_tabij) = 6._f64
185  matrix_poisson(index_tab, index_tabi1j) = -1._f64
186  matrix_poisson(index_tab, index_tabi1j1) = -1._f64
187  else if (k1 == -num_cells + 1) then! edge bottom left
188  matrix_poisson(index_tab, index_tabij_1) = -1._f64
189  matrix_poisson(index_tab, index_tabij) = 6._f64
190  matrix_poisson(index_tab, index_tabij1) = -1._f64
191  matrix_poisson(index_tab, index_tabi1j) = -1._f64
192  matrix_poisson(index_tab, index_tabi1j1) = -1._f64
193  else if (k2 == -num_cells + 1) then! edge bottom right
194  matrix_poisson(index_tab, index_tabi_1j) = -1._f64
195  matrix_poisson(index_tab, index_tabij) = 6._f64
196  matrix_poisson(index_tab, index_tabij1) = -1._f64
197  matrix_poisson(index_tab, index_tabi1j) = -1._f64
198  matrix_poisson(index_tab, index_tabi1j1) = -1._f64
199  else if (k1*k2 < 0 .and. k1 - k2 == num_cells - 1) then! edge right
200  matrix_poisson(index_tab, index_tabi_1j_1) = -1._f64
201  matrix_poisson(index_tab, index_tabi_1j) = -1._f64
202  matrix_poisson(index_tab, index_tabij) = 6._f64
203  matrix_poisson(index_tab, index_tabij1) = -1._f64
204  matrix_poisson(index_tab, index_tabi1j1) = -1._f64
205  end if
206 
207  else ! general case
208 
209  matrix_poisson(index_tab, index_tabi_1j_1) = -1._f64
210  matrix_poisson(index_tab, index_tabi_1j) = -1._f64
211  matrix_poisson(index_tab, index_tabij_1) = -1._f64
212  matrix_poisson(index_tab, index_tabij) = 6._f64
213  matrix_poisson(index_tab, index_tabij1) = -1._f64
214  matrix_poisson(index_tab, index_tabi1j) = -1._f64
215  matrix_poisson(index_tab, index_tabi1j1) = -1._f64
216 
217  end if
218  end do
219 
220  end subroutine sll_s_hex_matrix_poisson
221 
222  subroutine sll_s_hex_second_terme_poisson(second_terme, mesh, rho)
223  type(sll_t_hex_mesh_2d), pointer :: mesh
224  sll_real64, dimension(:), intent(out) :: second_terme
225  sll_real64, dimension(:), intent(in) :: rho
226  sll_int32 :: num_cells, i, index_tab, k1, k2
227  sll_real64 :: step
228  sll_real64 :: f1, f2, f3, f4, f5, f6
229  sll_int32 :: global
230 
231  num_cells = mesh%num_cells
232  step = mesh%delta**2*1.5_f64
233 
234  !************************
235  ! general case
236  !************************
237 
238  do global = 1, mesh%num_pts_tot
239 
240  k1 = mesh%hex_coord(1, global)
241  k2 = mesh%hex_coord(2, global)
242 
243  call mesh%index_hex_to_global(k1, k2, index_tab)
244 
245  f1 = value_if_inside_rho(k1 + 1, k2, mesh, rho)
246  f2 = value_if_inside_rho(k1 + 1, k2 + 1, mesh, rho)
247  f3 = value_if_inside_rho(k1, k2 + 1, mesh, rho)
248  f4 = value_if_inside_rho(k1 - 1, k2, mesh, rho)
249  f5 = value_if_inside_rho(k1 - 1, k2 - 1, mesh, rho)
250  f6 = value_if_inside_rho(k1, k2 - 1, mesh, rho)
251 
252  second_terme(index_tab) = (0.75_f64*rho(global) + &
253  (f1 + f2 + f3 + f4 + f5 + f6)/24._f64)*step ! ordre 4
254 
255  !second_terme(index_tab) = rho(global) * step ! order 2
256 
257  end do
258 
259  !************************
260  ! Boundaries
261  !************************
262 
263  ! corners of the hexagon
264 
265  call mesh%index_hex_to_global(num_cells, 0, index_tab)
266  second_terme(index_tab) = 0._f64
267  call mesh%index_hex_to_global(num_cells, num_cells, index_tab)
268  second_terme(index_tab) = 0._f64
269  call mesh%index_hex_to_global(0, num_cells, index_tab)
270  second_terme(index_tab) = 0._f64
271  call mesh%index_hex_to_global(-num_cells, 0, index_tab)
272  second_terme(index_tab) = 0._f64
273  call mesh%index_hex_to_global(-num_cells, -num_cells, index_tab)
274  second_terme(index_tab) = 0._f64
275  call mesh%index_hex_to_global(0, -num_cells, index_tab)
276  second_terme(index_tab) = 0._f64
277 
278  ! edges of the hexagon
279 
280  do i = 1, num_cells - 1 !( 0 and num_cells are the corners )
281 
282  ! top right edge
283  call mesh%index_hex_to_global(num_cells, i, index_tab)
284  second_terme(index_tab) = 0._f64
285  ! top left edge
286  call mesh%index_hex_to_global(num_cells - i, num_cells, index_tab)
287  second_terme(index_tab) = 0._f64
288  ! left edge
289  call mesh%index_hex_to_global(-i, num_cells - i, index_tab)
290  second_terme(index_tab) = 0._f64
291  ! bottom left edge
292  call mesh%index_hex_to_global(-num_cells, -i, index_tab)
293  second_terme(index_tab) = 0._f64
294  ! bottom right edge
295  call mesh%index_hex_to_global(-i, -num_cells, index_tab)
296  second_terme(index_tab) = 0._f64
297  ! right edge
298  call mesh%index_hex_to_global(num_cells - i, -i, index_tab)
299  second_terme(index_tab) = 0._f64
300 
301  end do
302 
303  !************************
304  ! Boundary conditions -> 0 everywhere at the moment
305  !************************
306 
307  ! corners of the hexagon
308 
309  call mesh%index_hex_to_global(num_cells - 1, 1, index_tab)
310  second_terme(index_tab) = second_terme(index_tab) + 0._f64
311  call mesh%index_hex_to_global(num_cells - 1, num_cells - 1, index_tab)
312  second_terme(index_tab) = second_terme(index_tab) + 0._f64
313  call mesh%index_hex_to_global(1, num_cells - 1, index_tab)
314  second_terme(index_tab) = second_terme(index_tab) + 0._f64
315  call mesh%index_hex_to_global(-num_cells + 1, 1, index_tab)
316  second_terme(index_tab) = second_terme(index_tab) + 0._f64
317  call mesh%index_hex_to_global(-num_cells + 1, -num_cells + 1, index_tab)
318  second_terme(index_tab) = second_terme(index_tab) + 0._f64
319  call mesh%index_hex_to_global(1, -num_cells + 1, index_tab)
320  second_terme(index_tab) = second_terme(index_tab) + 0._f64
321 
322  ! edges of the hexagon
323 
324  do i = 2, num_cells - 2 !( 1 and num_cells-1 are the corners )
325 
326  ! top right edge
327  call mesh%index_hex_to_global(num_cells - 1, i, index_tab)
328  second_terme(index_tab) = second_terme(index_tab) + 0._f64
329  ! top left edge
330  call mesh%index_hex_to_global(num_cells - 1 - i, num_cells - 1, index_tab)
331  second_terme(index_tab) = second_terme(index_tab) + 0._f64
332  ! left edge
333  call mesh%index_hex_to_global(-i, num_cells - 1 - i, index_tab)
334  second_terme(index_tab) = second_terme(index_tab) + 0._f64
335  ! bottom left edge
336  call mesh%index_hex_to_global(-num_cells + 1, -i, index_tab)
337  second_terme(index_tab) = second_terme(index_tab) + 0._f64
338  ! bottom right edge
339  call mesh%index_hex_to_global(-i, -num_cells + 1, index_tab)
340  second_terme(index_tab) = second_terme(index_tab) + 0._f64
341  ! right edge
342  call mesh%index_hex_to_global(num_cells - 1 - i, -i, index_tab)
343  second_terme(index_tab) = second_terme(index_tab) + 0._f64
344 
345  end do
346 
347  end subroutine sll_s_hex_second_terme_poisson
348 
350  function value_if_inside_rho(k1, k2, mesh, rho) result(f)
351  type(sll_t_hex_mesh_2d), pointer :: mesh
352  sll_real64, dimension(:) :: rho
353  sll_int32 :: k1, k2, n
354  sll_real64 :: f
355 
356  if (abs(k1) > mesh%num_cells .or. abs(k2) > mesh%num_cells .or. &
357  (k1)*(k2) < 0 .and. (abs(k1) + abs(k2) > mesh%num_cells)) then
358  f = 0._f64 ! null dirichlet boundary condition
359  else
360  n = mesh%hex_to_global(k1, k2)
361  f = rho(n)
362  end if
363 
364  end function value_if_inside_rho
365 
366 ! subroutine to compute the fields and its derivatives from the results
367 ! of solving the poisson equation . Precision is of order 4
368 
369  subroutine sll_s_compute_hex_fields(mesh, uxn, uyn, dxuxn, dyuxn, dxuyn, dyuyn, phi, type)
370  type(sll_t_hex_mesh_2d), pointer :: mesh
371  sll_real64, dimension(:) :: uxn, uyn, phi, dxuxn, dyuxn, dxuyn, dyuyn
372  sll_int32, intent(in) :: type
373  sll_int32 :: i, h1, h2
374  sll_real64 :: r11, r12, r21, r22, det
375  sll_real64 :: phi2j1, phi1j2, phi_2j1, phi_2j_1, phi_1j_2, phi1j_2, phi2j_1, phi_1j2
376  sll_real64 :: phi_2j_2, phi2j_2, phi2j2, phi_2j2
377  sll_real64 :: phi_3j, phij_3, phij, phi1j1, phi_1j_1, phi1j_1, phi_1j1
378  sll_real64 :: phi_2j, phi_1j, phi1j, phi2j, phij_2, phij_1, phij1, phij2
379  sll_real64 :: uh1, uh2, uh1h1, uh2h2, uh1h2
380 
381  det = (mesh%r1_x1*mesh%r2_x2 - mesh%r1_x2*mesh%r2_x1)/mesh%delta
382 
383  r11 = +mesh%r2_x2/det
384  r12 = -mesh%r2_x1/det
385  r21 = -mesh%r1_x2/det
386  r22 = +mesh%r1_x1/det
387 
388  if (type == 1) then
389 
390  do i = 1, mesh%num_pts_tot
391 
392  h1 = mesh%hex_coord(1, i)
393  h2 = mesh%hex_coord(2, i)
394 
395  phi_1j_1 = value_if_inside_phi(h1 - 1, h2 - 1, mesh, phi)
396  phi1j1 = value_if_inside_phi(h1 + 1, h2 + 1, mesh, phi)
397  phi1j_1 = value_if_inside_phi(h1 + 1, h2 - 1, mesh, phi)
398  phi_1j1 = value_if_inside_phi(h1 - 1, h2 + 1, mesh, phi)
399 
400  phi_2j_2 = value_if_inside_phi(h1 - 2, h2 - 2, mesh, phi)
401  phi2j2 = value_if_inside_phi(h1 + 2, h2 + 2, mesh, phi)
402  phi2j_2 = value_if_inside_phi(h1 + 2, h2 - 2, mesh, phi)
403  phi_2j2 = value_if_inside_phi(h1 - 2, h2 + 2, mesh, phi)
404 
405  phi_1j_2 = value_if_inside_phi(h1 - 1, h2 - 2, mesh, phi)
406  phi1j2 = value_if_inside_phi(h1 + 1, h2 + 2, mesh, phi)
407  phi1j_2 = value_if_inside_phi(h1 + 1, h2 - 2, mesh, phi)
408  phi_1j2 = value_if_inside_phi(h1 - 1, h2 + 2, mesh, phi)
409 
410  phi_2j_1 = value_if_inside_phi(h1 - 2, h2 - 1, mesh, phi)
411  phi2j1 = value_if_inside_phi(h1 + 2, h2 + 1, mesh, phi)
412  phi2j_1 = value_if_inside_phi(h1 + 2, h2 - 1, mesh, phi)
413  phi_2j1 = value_if_inside_phi(h1 - 2, h2 + 1, mesh, phi)
414 
415  phi_3j = value_if_inside_phi(h1 - 3, h2, mesh, phi)
416  phi_2j = value_if_inside_phi(h1 - 2, h2, mesh, phi)
417  phi_1j = value_if_inside_phi(h1 - 1, h2, mesh, phi)
418  phi1j = value_if_inside_phi(h1 + 1, h2, mesh, phi)
419  phi2j = value_if_inside_phi(h1 + 2, h2, mesh, phi)
420 
421  phij_3 = value_if_inside_phi(h1, h2 - 3, mesh, phi)
422  phij_2 = value_if_inside_phi(h1, h2 - 2, mesh, phi)
423  phij_1 = value_if_inside_phi(h1, h2 - 1, mesh, phi)
424  phij = value_if_inside_phi(h1, h2, mesh, phi)
425  phij1 = value_if_inside_phi(h1, h2 + 1, mesh, phi)
426  phij2 = value_if_inside_phi(h1, h2 + 2, mesh, phi)
427 
428  ! uh1 = (phii - phii_1)/ (mesh%delta) ! order 1 - very bad
429  ! uh2 = (phij - phij_1)/ (mesh%delta)
430 
431  ! uh1 = ( phii1 - phii_1 ) / (2._f64*mesh%delta) ! order 2
432  ! uh2 = ( phij1 - phij_1 ) / (2._f64*mesh%delta)
433 
434  ! uh1 = ( phii1/3._f64 + phii/2._f64 - phii_1 + phii_2/6._f64 ) / (mesh%delta) ! order 3
435  ! uh2 = ( phij1/3._f64 + phij/2._f64 - phij_1 + phij_2/6._f64 ) / (mesh%delta)
436 
437  uh1 = (phi_2j + 8._f64*(phi1j - phi_1j) - phi2j)/(12._f64*mesh%delta) ! order 4
438  uh2 = (phij_2 + 8._f64*(phij1 - phij_1) - phij2)/(12._f64*mesh%delta)
439 
440  ! uh1 = ( -phii_3/30._f64 + 0.25_f64*phii_2 - phii_1 + phii/3._f64&
441  ! + phii1*0.5_f64 - phii2/20._f64 ) / (mesh%delta) ! order 5
442  ! uh2 = ( -phij_3/30._f64 + 0.25_f64*phij_2 - phij_1 + phij/3._f64&
443  ! + phij1*0.5_f64 - phij2/20._f64 ) / (mesh%delta)
444 
445  ! order 4 approximation made with the values in the x and y directions
446  !-> not as good
447  ! uxn(i) = - ( phiyj_2 + 8._f64 * ( - phiyj_1 + phiyj1 ) &
448  ! - phiyj2 ) / (12._f64*mesh%delta)
449  ! uyn(i) = + ( phixi_2 + 8._f64 * ( - phixi_1 + phixi1 ) &
450  ! - phixi2 ) / (12._f64*sqrt(3._f64)*mesh%delta)
451 
452  !approximation made with the values in the r1 and r2 directions
453 
454  uxn(i) = -(uh1*r12 + uh2*r22)
455  uyn(i) = +(uh1*r11 + uh2*r21)
456 
457  ! order 2 approximation of the second derivatives
458 
459  ! uh1h1 = ( phi1j - 2._f64*phij + phi_1j )/mesh%delta**2
460  ! uh2h2 = ( phij1 - 2._f64*phij + phij_1 )/mesh%delta**2
461  ! uh1h2 = ( phi1j1 - phi1j_1 - phi_1j1 + phi_1j_1 )/(4._f64*mesh%delta**2)
462 
463  ! directe approximation of the second order
464  ! dyuxn(i) = -(phi1j1 - 2._f64*phij + phi_1j_1)/(mesh%delta**2) ! -dyy phi
465  ! dxuyn(i) = +(phi1j_1- 2._f64*phij + phi_1j1 )/(3._f64*mesh%delta**2) ! dxx phi
466 
467  ! order 4 approximation of the second derivatives
468 
469  uh1h1 = (-phi_2j + 16._f64*(phi_1j + phi1j) - 30._f64*phij - phi2j)/(12._f64*mesh%delta**2)
470  uh2h2 = (-phij_2 + 16._f64*(phij_1 + phij1) - 30._f64*phij - phij2)/(12._f64*mesh%delta**2)
471  uh1h2 = (64._f64*(phi1j1 - phi1j_1 - phi_1j1 + phi_1j_1) + &
472  8._f64*(-phi2j1 - phi1j2 - phi_2j_1 - phi_1j_2 + phi_1j2 + phi_2j1 + phi1j_2 + phi2j_1) + &
473  phi2j2 - phi2j_2 - phi_2j2 + phi_2j_2)/(144._f64*mesh%delta**2)
474 
475  dxuxn(i) = -(uh1h1*r11*r12 + uh1h2*(r11*r22 + r12*r21) + uh2h2*r21*r22) ! -dxy
476  !dyuxn(i) = - (uh1h1*r12*r12+2._f64*uh1h2*r12*r22+uh2h2*r22*r22) ! -dyy phi
477  dyuyn(i) = +(uh1h1*r11*r12 + uh1h2*(r11*r22 + r12*r21) + uh2h2*r21*r22) ! dyx
478  dxuyn(i) = +(uh1h1*r11*r11 + 2._f64*uh1h2*r11*r21 + uh2h2*r21*r21) !dxx phi
479 
480  ! directe approximation of the fourth order ->more or less the same results
481  dyuxn(i) = -(-phi_2j_2 + 16._f64*(phi_1j_1 + phi1j1) - 30._f64*phij - phi2j2)/ &
482  (12._f64*mesh%delta**2) ! -dyy phi
483  !dxuyn(i) = +( -phi_2j2 + 16._f64*(phi1j_1+phi_1j1) - 30._f64*phij - phi2j_2)/&
484  ! (3._f64*12._f64*mesh%delta**2) ! +dxx phi
485 
486  ! -> by combining the two approach we get the best results
487 
488  ! advection circulaire
489 
490  ! uxn(i) = + mesh%cartesian_coord(2,i)
491  ! uyn(i) = - mesh%cartesian_coord(1,i)
492 
493  ! dxuxn(i) = 0
494  ! dxuyn(i) = -1
495 
496  ! dyuxn(i) = 1
497  ! dyuyn(i) = 0
498 
499  end do
500  end if
501 
502  end subroutine sll_s_compute_hex_fields
503 
504  function value_if_inside_phi(k1, k2, mesh, phi) result(f)
505  type(sll_t_hex_mesh_2d), pointer :: mesh
506  sll_real64, dimension(:) :: phi
507  sll_int32 :: k1, k2, n
508  sll_real64 :: f
509 
510  if (abs(k1) > mesh%num_cells .or. abs(k2) > mesh%num_cells .or. &
511  (k1)*(k2) < 0 .and. (abs(k1) + abs(k2) > mesh%num_cells)) then
512  f = 0._f64 ! null dirichlet boundary condition
513  else
514  n = mesh%hex_to_global(k1, k2)
515  f = phi(n)
516  end if
517 
518  end function value_if_inside_phi
519 
520 end module sll_m_hex_poisson
sll_m_hex_poisson
Definition: sll_m_hex_poisson.F90:1
sll_m_hex_poisson::sll_s_hex_second_terme_poisson
subroutine, public sll_s_hex_second_terme_poisson(second_terme, mesh, rho)
Definition: sll_m_hex_poisson.F90:223
sll_m_hex_poisson::value_if_inside_phi
real(kind=f64) function value_if_inside_phi(k1, k2, mesh, phi)
Definition: sll_m_hex_poisson.F90:505
sll_m_hex_poisson::sll_s_hex_matrix_poisson
subroutine, public sll_s_hex_matrix_poisson(matrix_poisson, mesh, type)
Definition: sll_m_hex_poisson.F90:34
sll_m_hex_poisson::value_if_inside_rho
real(kind=f64) function value_if_inside_rho(k1, k2, mesh, rho)
Definition: sll_m_hex_poisson.F90:351
sll_m_hex_poisson::sll_s_compute_hex_fields
subroutine, public sll_s_compute_hex_fields(mesh, uxn, uyn, dxuxn, dyuxn, dxuyn, dyuyn, phi, type)
Definition: sll_m_hex_poisson.F90:370
sll_m_hexagonal_meshes
Definition: sll_m_hexagonal_meshes.F90:23
sll_m_hexagonal_meshes::sll_f_new_hex_mesh_2d
type(sll_t_hex_mesh_2d) function, pointer, public sll_f_new_hex_mesh_2d(num_cells, center_x1, center_x2, r11, r12, r21, r22, r31, r32, radius, EXTRA_TABLES)
Creates and initializes a new hexagonal mesh.
Definition: sll_m_hexagonal_meshes.F90:191
sll_m_hexagonal_meshes::sll_t_hex_mesh_2d
2d hexagonal mesh
Definition: sll_m_hexagonal_meshes.F90:63
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1