Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_advection_2d_CSL.F90
Go to the documentation of this file.
1 !**************************************************************
2 ! Copyright INRIA
3 ! Authors :
4 ! CALVI project team
5 !
6 ! This code SeLaLib (for Semi-Lagrangian-Library)
7 ! is a parallel library for simulating the plasma turbulence
8 ! in a tokamak.
9 !
10 ! This software is governed by the CeCILL-B license
11 ! under French law and abiding by the rules of distribution
12 ! of free software. You can use, modify and redistribute
13 ! the software under the terms of the CeCILL-B license as
14 ! circulated by CEA, CNRS and INRIA at the following URL
15 ! "http://www.cecill.info".
16 !**************************************************************
17 
18 !in development
19 
21 module sll_m_advection_2d_csl
22 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23 #include "sll_assert.h"
24 #include "sll_memory.h"
25 #include "sll_working_precision.h"
26 
27  use sll_m_advection_2d_base, only: &
28  sll_c_advector_2d
29 
30  use sll_m_characteristics_2d_base, only: &
31  sll_c_characteristics_2d_base
32 
33  implicit none
34 
35  public :: &
36  sll_f_new_csl_2d_advector, &
37  sll_t_csl_2d_advector
38 
39  private
40 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
41 
42  type, extends(sll_c_advector_2d) :: sll_t_csl_2d_advector
43 
44  class(sll_c_characteristics_2d_base), pointer :: charac
45  sll_real64, dimension(:), pointer :: eta1_coords
46  sll_real64, dimension(:), pointer :: eta2_coords
47  sll_real64, dimension(:, :), pointer :: charac_feet1
48  sll_real64, dimension(:, :), pointer :: charac_feet2
49  sll_int32 :: npts1
50  sll_int32 :: npts2
51  sll_int32 :: nbmax
52  sll_real64, dimension(:, :), pointer :: tt
53  sll_real64, dimension(:, :), pointer :: tcell
54  sll_real64, dimension(:), pointer :: intx
55  sll_real64, dimension(:), pointer :: inty
56  sll_real64, dimension(:, :), pointer :: dir
57  sll_int32, dimension(:, :, :), pointer :: cell
58 
59  contains
60 
61  procedure, pass(adv) :: init => initialize_csl_2d_advector
62  procedure, pass(adv) :: advect_2d => csl_advect_2d
63 
64  end type sll_t_csl_2d_advector
65 
66 contains
67 
68  function sll_f_new_csl_2d_advector( &
69  charac, &
70  Npts1, &
71  Npts2, &
72  eta1_min, &
73  eta1_max, &
74  eta2_min, &
75  eta2_max, &
76  eta1_coords, &
77  eta2_coords) &
78  result(adv)
79 
80  type(sll_t_csl_2d_advector), pointer :: adv
81  class(sll_c_characteristics_2d_base), pointer :: charac
82  sll_int32, intent(in) :: npts1
83  sll_int32, intent(in) :: npts2
84  sll_real64, intent(in), optional :: eta1_min
85  sll_real64, intent(in), optional :: eta1_max
86  sll_real64, intent(in), optional :: eta2_min
87  sll_real64, intent(in), optional :: eta2_max
88  sll_real64, dimension(:), pointer, optional :: eta1_coords
89  sll_real64, dimension(:), pointer, optional :: eta2_coords
90  sll_int32 :: ierr
91 
92  sll_allocate(adv, ierr)
93 
94  call initialize_csl_2d_advector( &
95  adv, &
96  charac, &
97  npts1, &
98  npts2, &
99  eta1_min, &
100  eta1_max, &
101  eta2_min, &
102  eta2_max, &
103  eta1_coords, &
104  eta2_coords)
105 
106  end function sll_f_new_csl_2d_advector
107 
108  subroutine initialize_csl_2d_advector( &
109  adv, &
110  charac, &
111  Npts1, &
112  Npts2, &
113  eta1_min, &
114  eta1_max, &
115  eta2_min, &
116  eta2_max, &
117  eta1_coords, &
118  eta2_coords)
119  class(sll_t_csl_2d_advector) :: adv
120  class(sll_c_characteristics_2d_base), pointer :: charac
121  sll_int32, intent(in) :: npts1
122  sll_int32, intent(in) :: npts2
123  sll_real64, intent(in), optional :: eta1_min
124  sll_real64, intent(in), optional :: eta1_max
125  sll_real64, intent(in), optional :: eta2_min
126  sll_real64, intent(in), optional :: eta2_max
127  sll_real64, dimension(:), pointer, optional :: eta1_coords
128  sll_real64, dimension(:), pointer, optional :: eta2_coords
129  sll_int32 :: ierr
130  sll_int32 :: i
131  sll_real64 :: delta_eta1
132  sll_real64 :: delta_eta2
133 
134  adv%Npts1 = npts1
135  adv%Npts2 = npts2
136  adv%charac => charac
137  !SLL_ALLOCATE(adv%x1_mesh(Npts1),ierr)
138  !SLL_ALLOCATE(adv%x2_mesh(Npts2),ierr)
139  sll_allocate(adv%eta1_coords(npts1), ierr)
140  sll_allocate(adv%eta2_coords(npts2), ierr)
141 
142  sll_allocate(adv%charac_feet1(npts1, npts2), ierr)
143  sll_allocate(adv%charac_feet2(npts1, npts2), ierr)
144 
145  if (present(eta1_min) .and. present(eta1_max)) then
146  if (present(eta1_coords)) then
147  print *, '#provide either eta1_coords or eta1_min and eta1_max'
148  print *, '#and not both in subroutine initialize_csl_2d_advector'
149  stop
150  else
151  delta_eta1 = (eta1_max - eta1_min)/real(npts1 - 1, f64)
152  do i = 1, npts1
153  adv%eta1_coords(i) = eta1_min + real(i - 1, f64)*delta_eta1
154  end do
155  end if
156  else if (present(eta1_coords)) then
157  if (size(eta1_coords, 1) < npts1) then
158  print *, '#bad size for eta1_coords in initialize_csl_2d_advector'
159  stop
160  else
161  adv%eta1_coords(1:npts1) = eta1_coords(1:npts1)
162  end if
163  else
164  print *, '#Warning, we assume eta1_min = 0._f64 eta1_max = 1._f64'
165  delta_eta1 = 1._f64/real(npts1 - 1, f64)
166  do i = 1, npts1
167  adv%eta1_coords(i) = real(i - 1, f64)*delta_eta1
168  end do
169  end if
170 
171  if (present(eta2_min) .and. present(eta2_max)) then
172  if (present(eta2_coords)) then
173  print *, '#provide either eta2_coords or eta2_min and eta2_max'
174  print *, '#and not both in subroutine initialize_csl_2d_advector'
175  stop
176  else
177  delta_eta2 = (eta2_max - eta2_min)/real(npts2 - 1, f64)
178  do i = 1, npts2
179  adv%eta2_coords(i) = eta2_min + real(i - 1, f64)*delta_eta2
180  end do
181  end if
182  else if (present(eta2_coords)) then
183  if (size(eta2_coords, 1) < npts2) then
184  print *, '#bad size for eta2_coords in initialize_csl_2d_advector'
185  stop
186  else
187  adv%eta2_coords(1:npts2) = eta2_coords(1:npts2)
188  end if
189  else
190  print *, '#Warning, we assume eta2_min = 0._f64 eta2_max = 1._f64'
191  delta_eta2 = 1._f64/real(npts2 - 1, f64)
192  do i = 1, npts2
193  adv%eta2_coords(i) = real(i - 1, f64)*delta_eta2
194  end do
195  end if
196 
197  adv%nbmax = 30
198  !print*,f(1,33)
199  !return
200  sll_allocate(adv%tt(adv%nbmax, 2), ierr)
201  sll_allocate(adv%cell(2, adv%nbmax, 4), ierr)
202  sll_allocate(adv%tcell(adv%nbmax, 4), ierr)
203  sll_allocate(adv%intx(0:adv%nbmax), ierr)
204  sll_allocate(adv%inty(0:adv%nbmax), ierr)
205  sll_allocate(adv%dir(adv%nbmax, 2), ierr)
206 
207  end subroutine initialize_csl_2d_advector
208 
209  subroutine csl_advect_2d( &
210  adv, &
211  A1, &
212  A2, &
213  dt, &
214  input, &
215  output)
216  class(sll_t_csl_2d_advector) :: adv
217 
218  sll_real64, dimension(:, :), intent(in) :: a1
219  sll_real64, dimension(:, :), intent(in) :: a2
220  sll_real64, intent(in) :: dt
221  sll_real64, dimension(:, :), intent(in) :: input
222  sll_real64, dimension(:, :), intent(out) :: output
223  sll_int32 :: i
224  sll_int32 :: j
225  sll_int32 :: nc_x1
226  sll_int32 :: nc_x2
227  sll_real64 :: x1_min
228  sll_real64 :: x1_max
229  sll_real64 :: x2_min
230  sll_real64 :: x2_max
231  sll_real64 :: xx(4)
232  sll_real64 :: yy(4)
233  sll_int32 :: ii(4)
234  sll_int32 :: jj(4)
235  sll_int32 :: imin
236  sll_int32 :: jmin
237  sll_int32 :: imax
238  sll_int32 :: jmax
239  sll_int32 :: ell
240  sll_int32 :: ell1
241  sll_real64 :: xa
242  sll_real64 :: xb
243  sll_real64 :: ya
244  sll_real64 :: yb
245  sll_int32 :: i0
246  sll_int32 :: j0
247  sll_int32 :: i1
248  sll_int32 :: j1
249  sll_int32 :: s
250  sll_int32 :: nbx(4)
251  sll_int32 :: nby(4)
252  sll_int32 :: dirx
253  sll_int32 :: diry
254  sll_int32 :: k
255  sll_int32 :: sx
256  sll_int32 :: sy
257  sll_real64 :: xa_loc
258  sll_real64 :: ya_loc
259  sll_real64 :: xb_loc
260  sll_real64 :: yb_loc
261  sll_int32 :: i0_loc
262  sll_int32 :: j0_loc
263  sll_real64 :: res
264  sll_int32 :: minfl
265  sll_int32 :: maxfl
266 
267  !PN ADDED TO AVOID WARNING
268  sll_assert(size(input, 1) > 0)
269 
270  nc_x1 = adv%Npts1 - 1
271  nc_x2 = adv%Npts2 - 1
272 
273  x1_min = adv%eta1_coords(1)
274  x1_max = adv%eta1_coords(adv%Npts1)
275  x2_min = adv%eta2_coords(1)
276  x2_max = adv%eta2_coords(adv%Npts2)
277 
278  call adv%charac%compute_characteristics( &
279  a1, &
280  a2, &
281  dt, &
282  adv%eta1_coords, &
283  adv%eta2_coords, &
284  adv%charac_feet1, &
285  adv%charac_feet2)
286 
287  do j = 0, nc_x2 - 1
288  do i = 0, nc_x1 - 1
289 
290  output(i + 1, j + 1) = 0._f64
291 
292  xx(1) = adv%charac_feet1(i + 1, j + 1)
293  xx(2) = adv%charac_feet1(i + 2, j + 1)
294  xx(3) = adv%charac_feet1(i + 2, j + 2)
295  xx(4) = adv%charac_feet1(i + 1, j + 2)
296 
297  yy(1) = adv%charac_feet2(i + 1, j + 1)
298  yy(2) = adv%charac_feet2(i + 2, j + 1)
299  yy(3) = adv%charac_feet2(i + 2, j + 2)
300  yy(4) = adv%charac_feet2(i + 1, j + 2)
301 
302  xx(1) = (xx(1) - x1_min)/(x1_max - x1_min)*real(nc_x1, f64)
303  xx(2) = (xx(2) - x1_min)/(x1_max - x1_min)*real(nc_x1, f64)
304  xx(3) = (xx(3) - x1_min)/(x1_max - x1_min)*real(nc_x1, f64)
305  xx(4) = (xx(4) - x1_min)/(x1_max - x1_min)*real(nc_x1, f64)
306 
307  yy(1) = (yy(1) - x2_min)/(x2_max - x2_min)*real(nc_x2, f64)
308  yy(2) = (yy(2) - x2_min)/(x2_max - x2_min)*real(nc_x2, f64)
309  yy(3) = (yy(3) - x2_min)/(x2_max - x2_min)*real(nc_x2, f64)
310  yy(4) = (yy(4) - x2_min)/(x2_max - x2_min)*real(nc_x2, f64)
311 
312  xx = xx + 0.5_f64
313  yy = yy + 0.5_f64
314 
315  ii(1) = floor(xx(1))
316  ii(2) = floor(xx(2))
317  ii(3) = floor(xx(3))
318  ii(4) = floor(xx(4))
319 
320  jj(1) = floor(yy(1))
321  jj(2) = floor(yy(2))
322  jj(3) = floor(yy(3))
323  jj(4) = floor(yy(4))
324 
325  imin = min(ii(1), ii(2), ii(3), ii(4))
326  jmin = min(jj(1), jj(2), jj(3), jj(4))
327  imax = max(ii(1), ii(2), ii(3), ii(4))
328  jmax = max(jj(1), jj(2), jj(3), jj(4))
329 
330  do ell = 1, 4
331 
332  ell1 = ell + 1
333  if (ell1 == 5) then
334  ell1 = 1
335  end if
336  xa = xx(ell)
337  ya = yy(ell)
338  xb = xx(ell1)
339  yb = yy(ell1)
340  i0 = ii(ell)
341  j0 = jj(ell)
342  i1 = ii(ell1)
343  j1 = jj(ell1)
344 
345  s = 1
346  if (i0 < i1) then
347  do k = i0 + 1, i1
348  adv%tt(s, 1) = (real(k, f64) - xa)/(xb - xa)
349  s = s + 1
350  end do
351  dirx = 1
352  end if
353  if (i0 > i1) then
354  do k = i0, i1 + 1, -1
355  adv%tt(s, 1) = (real(k, f64) - xa)/(xb - xa)
356  s = s + 1
357  end do
358  dirx = -1
359  end if
360  nbx(ell) = s - 1
361  s = 1
362  if (j0 < j1) then
363  do k = j0 + 1, j1
364  adv%tt(s, 2) = (real(k, f64) - ya)/(yb - ya)
365  s = s + 1
366  end do
367  diry = 1
368  end if
369  if (j0 > j1) then
370  do k = j0, j1 + 1, -1
371  adv%tt(s, 2) = (real(k, f64) - ya)/(yb - ya)
372  s = s + 1
373  end do
374  diry = -1
375  end if
376  nby(ell) = s - 1
377 
378  adv%cell(1, 1, ell) = i0
379  adv%cell(2, 1, ell) = j0
380  adv%tcell(1, ell) = 0._f64
381  sx = 1
382  sy = 1
383  s = 1
384 
385  do while ((sx <= nbx(ell)) .and. (sy <= nby(ell)))
386  if (adv%tt(sx, 1) < adv%tt(sy, 2)) then
387  s = s + 1
388  adv%cell(1, s, ell) = adv%cell(1, s - 1, ell) + dirx
389  adv%cell(2, s, ell) = adv%cell(2, s - 1, ell)
390  adv%tcell(s, ell) = adv%tt(sx, 1)
391  adv%intx(2*adv%cell(1, s - 1, ell) + (dirx - 1)/2 - 2*imin) = ya + adv%tt(sx, 1)*(yb - ya)
392  adv%dir(s, 1) = 1.0_f64
393  adv%dir(s, 2) = real(dirx, f64)
394  sx = sx + 1
395  else
396  s = s + 1
397  adv%cell(1, s, ell) = adv%cell(1, s - 1, ell)
398  adv%cell(2, s, ell) = adv%cell(2, s - 1, ell) + diry
399  adv%tcell(s, ell) = adv%tt(sy, 2)
400  adv%inty(2*adv%cell(2, s - 1, ell) + (diry - 1)/2 - 2*jmin) = xa + adv%tt(sy, 2)*(xb - xa)
401  adv%dir(s, 1) = 2.0_f64
402  adv%dir(s, 2) = real(diry, f64)
403  sy = sy + 1
404  end if
405  end do
406  do while (sx <= nbx(ell))
407  s = s + 1
408  adv%cell(1, s, ell) = adv%cell(1, s - 1, ell) + dirx
409  adv%cell(2, s, ell) = adv%cell(2, s - 1, ell)
410  adv%tcell(s, ell) = adv%tt(sx, 1)
411  adv%intx(2*adv%cell(1, s - 1, ell) + (dirx - 1)/2 - 2*imin) = ya + adv%tt(sx, 1)*(yb - ya)
412  adv%dir(s, 1) = 1.0_f64
413  adv%dir(s, 2) = real(dirx, f64)
414  sx = sx + 1
415  end do
416  do while (sy <= nby(ell))
417  s = s + 1
418  adv%cell(1, s, ell) = adv%cell(1, s - 1, ell)
419  adv%cell(2, s, ell) = adv%cell(2, s - 1, ell) + diry
420  adv%tcell(s, ell) = adv%tt(sy, 2)
421  adv%inty(2*adv%cell(2, s - 1, ell) + (diry - 1)/2 - 2*jmin) = xa + adv%tt(sy, 2)*(xb - xa)
422  adv%dir(s, 1) = 2.0_f64
423  adv%dir(s, 2) = real(diry, f64)
424  sy = sy + 1
425  end do
426 
427 !Computation of extern edges
428 
429  if ((ell == 1) .or. (ell == 4) .or. &
430  ((ell == 2) .and. (i == nc_x1)) .or. &
431  ((ell == 3) .and. (j == nc_x2))) then
432  xb_loc = xx(ell)
433  yb_loc = yy(ell)
434  do k = 2, nbx(ell) + nby(ell) + 2
435  xa_loc = xb_loc
436  ya_loc = yb_loc
437  i0_loc = adv%cell(1, k - 1, ell)
438  j0_loc = adv%cell(2, k - 1, ell)
439  if (k == nbx(ell) + nby(ell) + 2) then
440  xb_loc = xx(ell1)
441  yb_loc = yy(ell1)
442  else
443  if (adv%dir(k, 1) == 1) then
444  xb_loc = real(adv%cell(1, k, ell) + (1 - adv%dir(k, 2))/2, f64)
445  yb_loc = yy(ell) + adv%tcell(k, ell)*(yy(ell1) - yy(ell))
446  else
447  xb_loc = xx(ell) + adv%tcell(k, ell)*(xx(ell1) - xx(ell))
448  yb_loc = real(adv%cell(2, k, ell) + (1 - adv%dir(k, 2))/2, f64)
449  end if
450  end if
451 
452  res = 0._f64
453  !call calcule_coeff(N0,N1,buf2d,i0_loc,j0_loc,xA_loc,yA_loc,xB_loc,yB_loc,res,aretesh,aretesv,sommets,areteshb,areteshh,&
454  ! aretesvg,aretesvd,sommetsbg,sommetsbd,sommetshg,sommetshd,interp_case)
455 
456  output(i + 1, j + 1) = output(i + 1, j + 1) + res
457 
458  !if ((i-1>=0) .and. (ell==4)) then
459  ! output(i,j+1)=output(i,j+1)-res
460  !endif
461  !if ((j-1>=0) .and. (ell==1)) then
462  ! output(i+1,j)=output(i+1,j)-res
463  !endif
464  end do
465  end if
466  end do
467 
468 !#endif
469 !
470 !Computation of vertical intern edges
471 !
472 !#ifdef COMPUTE_V
473  do ell = 0, imax - imin - 1
474  minfl = min(floor(adv%intx(2*ell)), floor(adv%intx(2*ell + 1)))
475  maxfl = max(floor(adv%intx(2*ell)), floor(adv%intx(2*ell + 1)))
476 
477  i0_loc = imin + ell
478  yb_loc = min(adv%intx(2*ell), adv%intx(2*ell + 1))
479  do k = 0, maxfl - minfl
480  ya_loc = yb_loc
481  j0_loc = minfl + k
482  if (k == maxfl - minfl) then
483  yb_loc = max(adv%intx(2*ell), adv%intx(2*ell + 1))
484  else
485  yb_loc = real(minfl + k + 1, f64)
486  end if
487 
488  !call calcule_coeffv(N0,N1,buf2d,i0_loc,j0_loc,yA_loc,yB_loc,res,aretesh,aretesv,sommets,areteshb,areteshh,&
489  !aretesvg,aretesvd,sommetsbg,sommetsbd,sommetshg,sommetshd,interp_case)
490  output(i + 1, j + 1) = output(i + 1, j + 1) + res
491  end do
492  end do
493 !
494 !#endif
495 !
496 !Computation of horizontal intern edges
497 !
498 !#ifdef COMPUTE_X
499  do ell = 0, jmax - jmin - 1
500  minfl = min(floor(adv%inty(2*ell)), floor(adv%inty(2*ell + 1)))
501  maxfl = max(floor(adv%inty(2*ell)), floor(adv%inty(2*ell + 1)))
502 
503  j0_loc = jmin + ell
504  xa_loc = min(adv%inty(2*ell), adv%inty(2*ell + 1))
505  do k = 0, maxfl - minfl
506  xb_loc = xa_loc
507  i0_loc = minfl + k
508  if (k == maxfl - minfl) then
509  xa_loc = max(adv%inty(2*ell), adv%inty(2*ell + 1))
510  else
511  xa_loc = real(minfl + k + 1, f64)
512  end if
513 
514  !call calcule_coeffh(N0,N1,buf2d,i0_loc,j0_loc,xA_loc,xB_loc,res,aretesh,aretesv,sommets,areteshb,areteshh,&
515  !aretesvg,aretesvd,sommetsbg,sommetsbd,sommetshg,sommetshd,interp_case)
516  output(i + 1, j + 1) = output(i + 1, j + 1) + res
517 
518  end do
519  end do
520 
521 !#endif
522 ! deallocate(tnbr)
523 ! deallocate(tpts)
524 
525 ! enddo
526 !
527 ! enddo
528 
529 !#ifdef DIAG_TIME
530 ! f=buf2d
531 !#endif
532 ! deallocate(tt,cell,tcell,intx,inty,dir)
533 !
534 ! if ((interp_case==3) .or. (interp_case==6) .or. (interp_case==7)) then
535 ! deallocate(aretesh,aretesv,sommets)
536 ! deallocate(areteshb,areteshh,aretesvg,aretesvd,sommetsbg,sommetsbd,sommetshg,sommetshd)
537 ! endif
538 ! if (interp_case==4) then
539 ! deallocate(areteshb,areteshh,aretesvg,aretesvd,sommetsbg,sommetsbd,sommetshg,sommetshd)
540 ! endif
541 
542  end do
543  end do
544 
545  end subroutine csl_advect_2d
546 
547  subroutine aux(N0, N1, f, aretesh, aretesv, sommets, ordre, dom) !used in PPM case
548 
549  sll_int32, intent(in)::n0, n1, ordre
550  sll_real64, dimension(0:1, 0:1), intent(in)::dom
551  real(f64), dimension(0:N0, 0:N1 - 1), intent(in)::f
552  ! sll_real64,dimension(2,-1:N0+1,-1:N1+1),intent(in)::carac
553  sll_int32::i, j, im3, im2, im1, ib, ip1, ip2, jm3, jm2, jm1, jb, jp1, jp2
554  real(f64), dimension(0:N0, 0:N1 - 1), intent(inout)::aretesh, aretesv, sommets
555 ! print*,N0,N1,ordre
556 !stop
557 
558  if (dom(1, 1) == 1e-23) then
559  print *, '#is this possible?'
560  end if
561 
562  do i = 0, n0
563  do j = 0, n1 - 1
564  im3 = i - 3
565  im2 = i - 2
566  im1 = i - 1
567  ip1 = i + 1
568  ip2 = i + 2
569  if (i - 3 <= 0) im3 = 0;
570  if (i - 2 <= 0) im2 = 0;
571  if (i - 1 <= 0) im1 = 0;
572  ib = i;
573  if (i + 1 >= n0) ip1 = n0;
574  if (i + 2 >= n0) ip2 = n0;
575  jm3 = modulo(j - 3, n1)
576  jm2 = modulo(j - 2, n1)
577  jm1 = modulo(j - 1, n1)
578  jb = modulo(j, n1)
579  jp1 = modulo(j + 1, n1)
580  jp2 = modulo(j + 2, n1)
581  if (ordre == 1) then !PPM1
582  aretesv(i, j) = 7._f64/12._f64*(f(im1, jb) + f(ib, jb)) &
583  - 1._f64/12._f64*(f(im2, jb) + f(ip1, jb))
584  aretesh(i, j) = 7._f64/12._f64*(f(ib, jm1) + f(ib, jb)) &
585  - 1._f64/12._f64*(f(ib, jm2) + f(ib, jp1))
586  else if (ordre == 2) then !PPM2
587  aretesv(i, j) = 1._f64/60._f64*(f(ip2, jb) + f(im3, jb)) &
588  - 8._f64/60._f64*(f(ip1, jb) + f(im2, jb)) &
589  + 37._f64/60._f64*(f(ib, jb) + f(im1, jb))
590  aretesh(i, j) = 1._f64/60._f64*(f(ib, jp2) + f(ib, jm3)) &
591  - 8._f64/60._f64*(f(ib, jp1) + f(ib, jm2)) &
592  + 37._f64/60._f64*(f(ib, jb) + f(ib, jm1))
593  else if (ordre == 0) then !PPM0
594  aretesv(i, j) = 1._f64/2._f64*(f(ib, jb) + f(im1, jb))
595  aretesh(i, j) = 1._f64/2._f64*(f(ib, jb) + f(ib, jm1))
596  end if
597  end do
598  end do
599 
600  do i = 0, n0
601  do j = 0, n1 - 1
602  im3 = i - 3
603  im2 = i - 2
604  im1 = i - 1
605  ip1 = i + 1
606  ip2 = i + 2
607  if (i - 3 <= 0) im3 = 0;
608  if (i - 2 <= 0) im2 = 0;
609  if (i - 1 <= 0) im1 = 0;
610  ib = i;
611  if (i + 1 >= n0) ip1 = n0;
612  if (i + 2 >= n0) ip2 = n0;
613  jm3 = modulo(j - 3, n1)
614  jm2 = modulo(j - 2, n1)
615  jm1 = modulo(j - 1, n1)
616  jb = modulo(j, n1)
617  jp1 = modulo(j + 1, n1)
618  jp2 = modulo(j + 2, n1)
619 
620  if (ordre == 1) then !PPM1
621  sommets(i, j) = 7._f64/12._f64*(aretesv(ib, jm1) + aretesv(ib, jb)) &
622  - 1._f64/12._f64*(aretesv(ib, jm2) + aretesv(ib, jp1))
623  else if (ordre == 2) then !PPM2
624  sommets(i, j) = 1._f64/60._f64*(aretesv(ib, jp2) + aretesv(ib, jm3)) &
625  - 8._f64/60._f64*(aretesv(ib, jp1) + aretesv(ib, jm2)) &
626  + 37._f64/60._f64*(aretesv(ib, jb) + aretesv(ib, jm1))
627  else if (ordre == 0) then !PPM0
628  sommets(i, j) = 1._f64/2._f64*(aretesv(ib, jb) + aretesv(ib, jm1))
629  end if
630  end do
631  end do
632 
633  end subroutine aux
634 
635  subroutine compute_aux2_new( &
636  N_x1, &
637  N_x2, &
638  r_x1, &
639  s_x1, &
640  r_x2, &
641  s_x2, &
642  input, &
643  output)
644  sll_int32, intent(in) :: n_x1
645  sll_int32, intent(in) :: n_x2
646  sll_int32, intent(in) :: r_x1
647  sll_int32, intent(in) :: s_x1
648  sll_int32, intent(in) :: r_x2
649  sll_int32, intent(in) :: s_x2
650  sll_real64, dimension(:, :), intent(in) :: input
651  sll_real64, dimension(:, :, :, :), intent(out) :: output
652  !sll_real64, dimension(:), allocatable :: w_x1
653  sll_real64, dimension(:), allocatable :: ww_x1
654  !sll_real64, dimension(:), allocatable :: w_x2
655  sll_real64, dimension(:), allocatable :: ww_x2
656  sll_int32 :: ierr
657 
658  !PN ADDED TO AVOID WARNING
659  output = 0.0_f64
660  sll_assert(n_x1 > 0)
661  sll_assert(n_x2 > 0)
662  sll_assert(size(input, 1) > 0)
663  sll_assert(size(input, 2) > 0)
664 
665  !SLL_ALLOCATE(w_x1(r_x1,s_x1),ierr)
666  !SLL_ALLOCATE(w_x2(r_x2,s_x2),ierr)
667  sll_allocate(ww_x1(r_x1:s_x1 - 1), ierr)
668  sll_allocate(ww_x2(r_x2:s_x2 - 1), ierr)
669 
670  call compute_ww(ww_x1, r_x1, s_x1)
671  call compute_ww(ww_x2, r_x2, s_x2)
672 
673 ! tmp=0._f64
674 ! do ii=r,s-1
675 ! i3=i+ii-1;if(i3<=0)i3=0;if(i3>=N0)i3=N0
676 ! tmp=tmp+ww(r+s-1-ii)*f(i3,j)
677 ! enddo
678 ! aretesvg(i,j)=tmp
679 ! tmp=0._f64
680 ! do ii=r,s-1
681 ! i3=i+ii;if(i3<=0)i3=0;if(i3>=N0)i3=N0
682 ! tmp=tmp+ww(ii)*f(i3,j)
683 ! enddo
684 ! aretesvd(i,j)=tmp
685 ! tmp=0._f64
686 ! do ii=r,s-1
687 ! tmp=tmp+ww(ii)*f(i,modulo(j+ii,N1))
688 ! enddo
689 ! areteshh(i,j)=tmp
690 ! tmp=0._f64
691 ! do ii=r,s-1
692 ! tmp=tmp+ww(r+s-1-ii)*f(i,modulo(j+ii-1,N1))
693 ! enddo
694 ! areteshb(i,j)=tmp
695 
696  end subroutine compute_aux2_new
697 
698  subroutine compute_ww(ww, r, s)
699  sll_real64, dimension(:), intent(out) :: ww
700  sll_int32, intent(in) :: r
701  sll_int32, intent(in) :: s
702  sll_real64, dimension(:), allocatable :: w
703  sll_real64 :: tmp
704  sll_int32 :: i
705  sll_int32 :: j
706  sll_int32 :: ierr
707 
708  sll_allocate(w(r:s), ierr)
709 
710  do i = r, -1
711  tmp = 1._f64
712  do j = r, i - 1
713  tmp = tmp*real(i - j, f64)
714  end do
715  do j = i + 1, s
716  tmp = tmp*real(i - j, f64)
717  end do
718  tmp = 1._f64/tmp
719  do j = r, i - 1
720  tmp = tmp*real(-j, f64)
721  end do
722  do j = i + 1, -1
723  tmp = tmp*real(-j, f64)
724  end do
725  do j = 1, s
726  tmp = tmp*real(-j, f64)
727  end do
728  w(i) = tmp
729  end do
730 
731  do i = 1, s
732  tmp = 1._f64
733  do j = r, i - 1
734  tmp = tmp*real(i - j, f64)
735  end do
736  do j = i + 1, s
737  tmp = tmp*real(i - j, f64)
738  end do
739  tmp = 1._f64/tmp
740  do j = r, -1
741  tmp = tmp*real(-j, f64)
742  end do
743  do j = 1, i - 1
744  tmp = tmp*real(-j, f64)
745  end do
746  do j = i + 1, s
747  tmp = tmp*real(-j, f64)
748  end do
749  w(i) = tmp
750  end do
751 
752  tmp = 0._f64
753  do i = r, -1
754  tmp = tmp + w(i)
755  end do
756  do i = 1, s
757  tmp = tmp + w(i)
758  end do
759  w(0) = -tmp
760 
761  tmp = 0._f64
762  do i = r, -1
763  tmp = tmp + w(i)
764  ww(i) = -tmp
765  end do
766  tmp = 0._f64
767  do i = s, 1, -1
768  tmp = tmp + w(i)
769  ww(i - 1) = tmp
770  end do
771 
772  sll_deallocate_array(w, ierr)
773 
774  end subroutine compute_ww
775 
776  subroutine aux2(N0, N1, f, areteshb, areteshh, aretesvg, aretesvd, sommetsbg, sommetsbd, &
777  sommetshg, sommetshd, ordre, carac, dom) !used in PPM case
778  integer, intent(in)::N0, N1, ordre
779  real(f64), dimension(0:1, 0:1), intent(in)::dom
780  real(f64), dimension(0:N0, 0:N1 - 1), intent(in)::f
781  real(f64), dimension(2, -1:N0 + 1, -1:N1 + 1), intent(in)::carac
782  integer::i, j, i3
783  real(f64), dimension(0:N0, 0:N1 - 1), intent(inout)::areteshb, areteshh, aretesvg, aretesvd, &
784  sommetsbg, sommetsbd, sommetshg, sommetshd
785  real(f64) ::w(-ordre:ordre + 1), tmp, ww(-ordre:ordre)
786  integer::r, s, ii, d
787  !integer :: ib,im3,im2,im1,ip1,ip2,jm3,jm2,jm1,jb,jp1,jp2
788 
789  d = ordre
790  r = -d
791  s = d + 1
792  if (carac(1, 1, 1) == 1e-23) then
793  print *, '#is this possible?'
794  end if
795  if (dom(1, 1) == 1e-23) then
796  print *, '#is this possible?'
797  end if
798  !maple code for generation of w
799  !for k from r to -1 do
800  ! C[k]:=product((k-j),j=r..k-1)*product((k-j),j=k+1..s):
801  ! C[k]:=1/C[k]*product((-j),j=r..k-1)*product((-j),j=k+1..-1)*product((-j),j=1..s):
802  !od:
803  !for k from 1 to s do
804  ! C[k]:=product((k-j),j=r..k-1)*product((k-j),j=k+1..s):
805  ! C[k]:=1/C[k]*product((-j),j=r..-1)*product((-j),j=1..k-1)*product((-j),j=k+1..s):
806  !od:
807  !C[0]:=-add(C[k],k=r..-1)-add(C[k],k=1..s):
808 
809  do i = r, -1
810  tmp = 1._f64
811  do j = r, i - 1
812  tmp = tmp*real(i - j, f64)
813  end do
814  do j = i + 1, s
815  tmp = tmp*real(i - j, f64)
816  end do
817  tmp = 1._f64/tmp
818  do j = r, i - 1
819  tmp = tmp*real(-j, f64)
820  end do
821  do j = i + 1, -1
822  tmp = tmp*real(-j, f64)
823  end do
824  do j = 1, s
825  tmp = tmp*real(-j, f64)
826  end do
827  w(i) = tmp
828  end do
829 
830 ! do i=r,-1
831 ! tmp=1._f64
832 ! !do j=r,i-1
833 ! ! tmp=tmp*real(i-j,f64)
834 ! !enddo
835 ! !do j=i+1,s
836 ! ! tmp=tmp*real(i-j,f64)
837 ! !enddo
838 ! !tmp=1._f64/tmp
839 ! do j=r,i-1 !-j/(i-j)=j/(j-i)=1/(1-i/j)
840 ! tmp=tmp*(1._f64-real(i,f64)/real(j,f64))
841 ! enddo
842 ! do j=i+1,-1
843 ! tmp=tmp*(1._f64-real(i,f64)/real(j,f64))
844 ! enddo
845 ! do j=1,s
846 ! tmp=tmp*(1._f64-real(i,f64)/real(j,f64))
847 ! enddo
848 ! tmp=tmp*real(i,f64)
849 ! w(i)=1._f64/tmp
850 ! enddo
851 
852  do i = 1, s
853  tmp = 1._f64
854  do j = r, i - 1
855  tmp = tmp*real(i - j, f64)
856  end do
857  do j = i + 1, s
858  tmp = tmp*real(i - j, f64)
859  end do
860  tmp = 1._f64/tmp
861  do j = r, -1
862  tmp = tmp*real(-j, f64)
863  end do
864  do j = 1, i - 1
865  tmp = tmp*real(-j, f64)
866  end do
867  do j = i + 1, s
868  tmp = tmp*real(-j, f64)
869  end do
870  w(i) = tmp
871  end do
872 
873  tmp = 0._f64
874  do i = r, -1
875  tmp = tmp + w(i)
876  end do
877  do i = 1, s
878  tmp = tmp + w(i)
879  end do
880  w(0) = -tmp
881 
882  !print *,'w',w
883  !do ii=r,s
884  ! print *,ii,w(r+s-ii)
885  !enddo
886 
887  !compute now ww
888  !maple code
889  !#for conservative formulation
890  !tmp:=0:
891  !for k from r to -1 do
892  !tmp:=tmp+C[k]:
893  !CC[k]:=-tmp:
894  !od:
895  !tmp:=0:
896  !for k from s to 1 by -1 do
897  ! tmp:=tmp+C[k]:
898  ! CC[k-1]:=tmp:
899  !od:
900  !seq(CC[k],k=r..s-1);
901  !evalf(%);
902 
903  tmp = 0._f64
904  do i = r, -1
905  tmp = tmp + w(i)
906  ww(i) = -tmp
907  end do
908  tmp = 0._f64
909  do i = s, 1, -1
910  tmp = tmp + w(i)
911  ww(i - 1) = tmp
912  end do
913 
914  !print *,'ww',ww
915  !do ii=r,s-1
916  ! print *,ii,ww(r+s-1-ii)
917  !enddo
918  !stop
919  do j = 0, n1 - 1
920  do i = 0, n0
921  tmp = 0._f64
922  do ii = r, s - 1
923  i3 = i + ii - 1; if (i3 <= 0) i3 = 0; if (i3 >= n0) i3 = n0
924  tmp = tmp + ww(r + s - 1 - ii)*f(i3, j)
925  end do
926  aretesvg(i, j) = tmp
927  tmp = 0._f64
928  do ii = r, s - 1
929  i3 = i + ii; if (i3 <= 0) i3 = 0; if (i3 >= n0) i3 = n0
930  tmp = tmp + ww(ii)*f(i3, j)
931  end do
932  aretesvd(i, j) = tmp
933  tmp = 0._f64
934  do ii = r, s - 1
935  tmp = tmp + ww(ii)*f(i, modulo(j + ii, n1))
936  end do
937  areteshh(i, j) = tmp
938  tmp = 0._f64
939  do ii = r, s - 1
940  tmp = tmp + ww(r + s - 1 - ii)*f(i, modulo(j + ii - 1, n1))
941  end do
942  areteshb(i, j) = tmp
943  end do
944  end do
945 
946  do j = 0, n1 - 1
947  do i = 0, n0
948  tmp = 0._f64
949  do ii = r, s - 1
950  i3 = i + ii; if (i3 <= 0) i3 = 0; if (i3 >= n0) i3 = n0
951  tmp = tmp + ww(ii)*areteshh(i3, j)
952  end do
953  sommetshd(i, j) = tmp
954  tmp = 0._f64
955  do ii = r, s - 1
956  i3 = i + ii - 1; if (i3 <= 0) i3 = 0; if (i3 >= n0) i3 = n0
957  tmp = tmp + ww(r + s - 1 - ii)*areteshh(i3, j)
958  end do
959  sommetshg(i, j) = tmp
960  tmp = 0._f64
961  do ii = r, s - 1
962  i3 = i + ii; if (i3 <= 0) i3 = 0; if (i3 >= n0) i3 = n0
963  tmp = tmp + ww(ii)*areteshb(i3, j)
964  end do
965  sommetsbd(i, j) = tmp
966  tmp = 0._f64
967  do ii = r, s - 1
968  i3 = i + ii - 1; if (i3 <= 0) i3 = 0; if (i3 >= n0) i3 = n0
969  tmp = tmp + ww(r + s - 1 - ii)*areteshb(i3, j)
970  end do
971  sommetsbg(i, j) = tmp
972  end do
973  end do
974 
975  end subroutine aux2
976 
977 end module sll_m_advection_2d_csl
sll_m_advection_2d_base
Definition: sll_m_advection_2d_base.F90:19
sll_m_advection_2d_csl
Definition: sll_m_advection_2d_CSL.F90:21
sll_m_advection_2d_csl::sll_f_new_csl_2d_advector
type(sll_t_csl_2d_advector) function, pointer, public sll_f_new_csl_2d_advector(charac, Npts1, Npts2, eta1_min, eta1_max, eta2_min, eta2_max, eta1_coords, eta2_coords)
Definition: sll_m_advection_2d_CSL.F90:79
sll_m_advection_2d_csl::initialize_csl_2d_advector
subroutine initialize_csl_2d_advector(adv, charac, Npts1, Npts2, eta1_min, eta1_max, eta2_min, eta2_max, eta1_coords, eta2_coords)
Definition: sll_m_advection_2d_CSL.F90:119
sll_m_advection_2d_csl::compute_ww
subroutine compute_ww(ww, r, s)
Definition: sll_m_advection_2d_CSL.F90:699
sll_m_advection_2d_csl::aux2
subroutine aux2(N0, N1, f, areteshb, areteshh, aretesvg, aretesvd, sommetsbg, sommetsbd, sommetshg, sommetshd, ordre, carac, dom)
Definition: sll_m_advection_2d_CSL.F90:778
sll_m_advection_2d_csl::aux
subroutine aux(N0, N1, f, aretesh, aretesv, sommets, ordre, dom)
Definition: sll_m_advection_2d_CSL.F90:548
sll_m_advection_2d_csl::compute_aux2_new
subroutine compute_aux2_new(N_x1, N_x2, r_x1, s_x1, r_x2, s_x2, input, output)
Definition: sll_m_advection_2d_CSL.F90:644
sll_m_advection_2d_csl::csl_advect_2d
subroutine csl_advect_2d(adv, A1, A2, dt, input, output)
Definition: sll_m_advection_2d_CSL.F90:216
sll_m_characteristics_2d_base
Abstract class to compute the characteristic in two dimensions.
Definition: sll_m_characteristics_2d_base.F90:20
sll_m_advection_2d_base::sll_c_advector_2d
Definition: sll_m_advection_2d_base.F90:37
sll_m_advection_2d_csl::sll_t_csl_2d_advector
Definition: sll_m_advection_2d_CSL.F90:42
sll_m_characteristics_2d_base::sll_c_characteristics_2d_base
Definition: sll_m_characteristics_2d_base.F90:37
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:41 for Semi-Lagrangian Library by   doxygen 1.9.1