Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_characteristics_2d_explicit_euler_conservative.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 
21 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22 #include "sll_assert.h"
23 #include "sll_memory.h"
24 #include "sll_working_precision.h"
25 
27  sll_p_periodic, &
28  sll_p_set_to_limit, &
29  sll_p_user_defined
30 
36 
37  implicit none
38 
39  public :: &
41 
42  private
43 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
44 
45  type, extends(sll_c_characteristics_2d_base) :: &
47  sll_int32 :: npts1
48  sll_int32 :: npts2
49  sll_real64 :: eta1_min
50  sll_real64 :: eta1_max
51  sll_real64 :: eta2_min
52  sll_real64 :: eta2_max
53  procedure(sll_i_signature_process_outside_point), pointer, nopass :: &
54  process_outside_point1
55  procedure(sll_i_signature_process_outside_point), pointer, nopass :: &
56  process_outside_point2
57  sll_real64, dimension(:, :, :), pointer :: buf_3d
58  sll_int32 :: bc_type_1
59  sll_int32 :: bc_type_2
60 
61  contains
62  !function, pass(charac) :: new => &
63  ! sll_f_new_explicit_euler_conservative_2d_charac
64  procedure, pass(charac) :: initialize => &
66  procedure, pass(charac) :: compute_characteristics => &
69 
70 contains
72  Npts1, &
73  Npts2, &
74  bc_type_1, &
75  bc_type_2, &
76  eta1_min, &
77  eta1_max, &
78  eta2_min, &
79  eta2_max, &
80  process_outside_point1, &
81  process_outside_point2) &
82  result(charac)
83 
84  type(explicit_euler_conservative_2d_charac_computer), pointer :: charac
85  sll_int32, intent(in) :: npts1
86  sll_int32, intent(in) :: npts2
87  sll_int32, intent(in), optional :: bc_type_1
88  sll_int32, intent(in), optional :: bc_type_2
89  sll_real64, intent(in), optional :: eta1_min
90  sll_real64, intent(in), optional :: eta1_max
91  sll_real64, intent(in), optional :: eta2_min
92  sll_real64, intent(in), optional :: eta2_max
93  procedure(sll_i_signature_process_outside_point), optional :: &
94  process_outside_point1
95  procedure(sll_i_signature_process_outside_point), optional :: &
96  process_outside_point2
97  sll_int32 :: ierr
98 
99  sll_allocate(charac, ierr)
101  charac, &
102  npts1, &
103  npts2, &
104  bc_type_1, &
105  bc_type_2, &
106  eta1_min, &
107  eta1_max, &
108  eta2_min, &
109  eta2_max, &
110  process_outside_point1, &
111  process_outside_point2)
112 
114 
116  charac, &
117  Npts1, &
118  Npts2, &
119  bc_type_1, &
120  bc_type_2, &
121  eta1_min, &
122  eta1_max, &
123  eta2_min, &
124  eta2_max, &
125  process_outside_point1, &
126  process_outside_point2)
127 
129  sll_int32, intent(in) :: npts1
130  sll_int32, intent(in) :: npts2
131  sll_int32, intent(in), optional :: bc_type_1
132  sll_int32, intent(in), optional :: bc_type_2
133  sll_real64, intent(in), optional :: eta1_min
134  sll_real64, intent(in), optional :: eta1_max
135  sll_real64, intent(in), optional :: eta2_min
136  sll_real64, intent(in), optional :: eta2_max
137  procedure(sll_i_signature_process_outside_point), optional :: &
138  process_outside_point1
139  procedure(sll_i_signature_process_outside_point), optional :: &
140  process_outside_point2
141 
142  sll_int32 :: ierr
143 
144  charac%Npts1 = npts1
145  charac%Npts2 = npts2
146 
147  sll_allocate(charac%buf_3d(2, 0:npts1, 0:npts2), ierr)
148 
149  if (present(eta1_min)) then
150  charac%eta1_min = eta1_min
151  else
152  charac%eta1_min = 0._f64
153  end if
154  if (present(eta1_max)) then
155  charac%eta1_max = eta1_max
156  else
157  charac%eta1_max = 1._f64
158  end if
159  if (present(eta2_min)) then
160  charac%eta2_min = eta2_min
161  else
162  charac%eta2_min = 0._f64
163  end if
164 
165  if (present(eta2_max)) then
166  charac%eta2_max = eta2_max
167  else
168  charac%eta2_max = 1._f64
169  end if
170 
171  !charac%process_outside_point1 => process_outside_point1
172  !charac%process_outside_point2 => process_outside_point2
173 
174  if (present(process_outside_point1)) then
175  charac%process_outside_point1 => process_outside_point1
176  charac%bc_type_1 = sll_p_user_defined
177  else if (.not. (present(bc_type_1))) then
178  print *, '#provide boundary condition'
179  print *, '#bc_type_1 or process_outside_point1 function'
180  print *, '#in initialize_explicit_euler_conservative_2d_charac'
181  stop
182  else
183  charac%bc_type_1 = bc_type_1
184  select case (bc_type_1)
185  case (sll_p_periodic)
186  charac%process_outside_point1 => sll_f_process_outside_point_periodic
187  case (sll_p_set_to_limit)
188  charac%process_outside_point1 => sll_f_process_outside_point_set_to_limit
189  case default
190  print *, '#bad value of boundary condition'
191  print *, '#in initialize_explicit_euler_conservative_2d_charac_computer'
192  stop
193  end select
194  end if
195 
196  if ((present(process_outside_point1)) .and. (present(bc_type_1))) then
197  print *, '#provide either process_outside_point1 or bc_type_1'
198  print *, '#and not both'
199  print *, '#in initialize_explicit_euler_conservative_2d_charac_computer'
200  stop
201  end if
202 
203  if (present(process_outside_point2)) then
204  charac%process_outside_point2 => process_outside_point2
205  charac%bc_type_2 = sll_p_user_defined
206  else if (.not. (present(bc_type_2))) then
207  print *, '#provide boundary condition'
208  print *, '#bc_type_2 or process_outside_point1 function'
209  stop
210  else
211  charac%bc_type_2 = bc_type_2
212  select case (bc_type_2)
213  case (sll_p_periodic)
214  charac%process_outside_point2 => sll_f_process_outside_point_periodic
215  case (sll_p_set_to_limit)
216  charac%process_outside_point2 => sll_f_process_outside_point_set_to_limit
217  case default
218  print *, '#bad value of boundary condition'
219  print *, '#in initialize_explicit_euler_conservative_2d_charac_computer'
220  stop
221  end select
222  end if
223 
224  if ((present(process_outside_point2)) .and. (present(bc_type_2))) then
225  print *, '#provide either process_outside_point2 or bc_type_2'
226  print *, '#and not both'
227  print *, '#in initialize_explicit_euler_conservative_2d_charac_computer'
228  stop
229  end if
230 
232 
234  charac, &
235  A1, &
236  A2, &
237  dt, &
238  input1, &
239  input2, &
240  output1, &
241  output2)
242 
244  sll_real64, dimension(:, :), intent(in) :: a1
245  sll_real64, dimension(:, :), intent(in) :: a2
246  sll_real64, intent(in) :: dt
247  sll_real64, dimension(:), intent(in) :: input1
248  sll_real64, dimension(:), intent(in) :: input2
249  sll_real64, dimension(:, :), intent(out) :: output1
250  sll_real64, dimension(:, :), intent(out) :: output2
251  sll_int32 :: i
252  sll_int32 :: j
253  sll_int32 :: npts1
254  sll_int32 :: npts2
255  sll_real64 :: eta1_min
256  sll_real64 :: eta1_max
257  sll_real64 :: eta2_min
258  sll_real64 :: eta2_max
259 
260  npts1 = charac%Npts1
261  npts2 = charac%Npts2
262  eta1_min = charac%eta1_min
263  eta1_max = charac%eta1_max
264  eta2_min = charac%eta2_min
265  eta2_max = charac%eta2_max
266 
267  sll_assert(size(a1, 1) >= charac%Npts1 - 1)
268  sll_assert(size(a1, 2) >= charac%Npts2 - 1)
269  sll_assert(size(a2, 1) >= charac%Npts1 - 1)
270  sll_assert(size(a2, 2) >= charac%Npts2 - 1)
271  sll_assert(size(input1) >= charac%Npts1)
272  sll_assert(size(input2) >= charac%Npts2)
273  sll_assert(size(output1, 1) >= charac%Npts1)
274  sll_assert(size(output1, 2) >= charac%Npts2)
275  sll_assert(size(output2, 1) >= charac%Npts1)
276  sll_assert(size(output2, 2) >= charac%Npts2)
277 
278  do j = 1, npts2 - 1
279  do i = 1, npts1 - 1
280  charac%buf_3d(1, i, j) = 0.5_f64*(input1(i) + input1(i + 1)) - dt*a1(i, j)
281  !if((output1(i,j)<=eta1_min).or.(output1(i,j)>=eta1_max))then
282  ! output1(i,j)=charac%process_outside_point1(output1(i,j),eta1_min,eta1_max)
283  !endif
284  charac%buf_3d(2, i, j) = 0.5_f64*(input2(j) + input2(j + 1)) - dt*a2(i, j)
285  !if((output2(i,j)<=eta2_min).or.(output2(i,j)>=eta2_max))then
286  ! output2(i,j)=charac%process_outside_point2(output2(i,j),eta2_min,eta2_max)
287  !endif
288  end do
289  end do
290 
291  select case (charac%bc_type_1)
292  case (sll_p_periodic)
293  do j = 1, npts2 - 1
294  charac%buf_3d(1, 0, j) = charac%buf_3d(1, npts1 - 1, j) - (eta1_max - eta1_min)
295  charac%buf_3d(2, 0, j) = charac%buf_3d(2, npts1 - 1, j)
296  charac%buf_3d(1, npts1, j) = charac%buf_3d(1, 1, j) + (eta1_max - eta1_min)
297  charac%buf_3d(2, npts1, j) = charac%buf_3d(2, 1, j)
298  end do
299  case (sll_p_set_to_limit)
300  do j = 1, npts2 - 1
301  charac%buf_3d(1, 0, j) = 2._f64*eta1_min - charac%buf_3d(1, 1, j)
302  charac%buf_3d(2, 0, j) = charac%buf_3d(2, 1, j)
303  charac%buf_3d(1, npts1, j) = 2._f64*eta1_max - charac%buf_3d(1, npts1 - 1, j)
304  charac%buf_3d(2, npts1, j) = charac%buf_3d(2, npts1 - 1, j)
305  end do
306  case default
307  print *, '#bad value for charac%bc_type_1'
308  stop
309  end select
310 
311  select case (charac%bc_type_2)
312  case (sll_p_periodic)
313  do i = 0, npts1
314  charac%buf_3d(2, i, 0) = charac%buf_3d(2, i, npts2 - 1) - (eta2_max - eta2_min)
315  charac%buf_3d(1, i, 0) = charac%buf_3d(1, i, npts2 - 1)
316  charac%buf_3d(2, i, npts1) = charac%buf_3d(2, i, 1) + (eta2_max - eta2_min)
317  charac%buf_3d(1, i, npts1) = charac%buf_3d(1, i, 1)
318  end do
319  case (sll_p_set_to_limit)
320  do i = 0, npts1
321  charac%buf_3d(2, i, 0) = 2._f64*eta2_min - charac%buf_3d(2, i, 1)
322  charac%buf_3d(1, i, 0) = charac%buf_3d(1, i, 1)
323  charac%buf_3d(2, i, npts1) = 2._f64*eta2_max - charac%buf_3d(2, i, npts1 - 1)
324  charac%buf_3d(1, i, npts1) = charac%buf_3d(1, i, npts1 - 1)
325  end do
326  case default
327  print *, '#bad value for charac%bc_type_2'
328  stop
329  end select
330 
331  do j = 1, npts2
332  do i = 1, npts1
333  output1(i, j) = &
334  0.25_f64*(charac%buf_3d(1, i, j) &
335  + charac%buf_3d(1, i - 1, j) &
336  + charac%buf_3d(1, i - 1, j - 1) &
337  + charac%buf_3d(1, i, j - 1))
338  output2(i, j) = &
339  0.25_f64*(charac%buf_3d(2, i, j) &
340  + charac%buf_3d(2, i - 1, j) &
341  + charac%buf_3d(2, i - 1, j - 1) &
342  + charac%buf_3d(2, i, j - 1))
343  end do
344  end do
345 
347 
Abstract class to compute the characteristic in two dimensions.
real(kind=f64) function, public sll_f_process_outside_point_periodic(eta, eta_min, eta_max)
real(kind=f64) function, public sll_f_process_outside_point_set_to_limit(eta, eta_min, eta_max)
computes the characteristic with explicit euler conservative scheme
subroutine compute_explicit_euler_conservative_2d_charac(charac, A1, A2, dt, input1, input2, output1, output2)
subroutine initialize_explicit_euler_conservative_2d_charac(charac, Npts1, Npts2, bc_type_1, bc_type_2, eta1_min, eta1_max, eta2_min, eta2_max, process_outside_point1, process_outside_point2)
type(explicit_euler_conservative_2d_charac_computer) function, pointer, public sll_f_new_explicit_euler_conservative_2d_charac(Npts1, Npts2, bc_type_1, bc_type_2, eta1_min, eta1_max, eta2_min, eta2_max, process_outside_point1, process_outside_point2)
    Report Typos and Errors