Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_pic_visu.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 module sll_m_pic_visu
23 
24 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
25 #include "sll_assert.h"
26 #include "sll_working_precision.h"
27 
28  use sll_m_ascii_io, only: &
29  sll_s_ascii_file_create
30 
31  use sll_m_gnuplot, only: &
32  sll_o_gnuplot_2d
33 
34  use sll_m_utilities, only: &
35  sll_s_int2string, &
36  sll_s_new_file_id
37 
38  use sll_m_xdmf, only: &
39  sll_s_xdmf_corect2d_nodes
40 
41  implicit none
42 
43  public :: &
44  sll_s_compute_df_cic, &
45  sll_o_distribution_gnuplot, &
46  sll_s_distribution_m4_gnuplot, &
47  sll_s_distribution_tsc_gnuplot, &
48  sll_s_distribution_xdmf, &
49  sll_s_electricpotential_gnuplot_inline, &
50  sll_s_energies_electrostatic_gnuplot_inline, &
51  sll_o_particles_center_gnuplot, &
52  sll_s_particles_center_gnuplot_inline, &
53  sll_o_plot_format_points3d, &
54  sll_s_plot_format_xmdv
55 
56  private
57 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
58 
60  interface sll_o_particles_center_gnuplot
61  module procedure xv_particles_center_gnuplot
62  end interface sll_o_particles_center_gnuplot
63 
65  interface sll_o_distribution_gnuplot
66  module procedure distribution_xv_gnuplot
67  end interface sll_o_distribution_gnuplot
68 
70  interface sll_o_plot_format_points3d
71  module procedure pq_plot_format_points3d
72  module procedure pqr_plot_format_points3d
73  module procedure pqrs_plot_format_points3d
74  end interface sll_o_plot_format_points3d
75 
76 contains
77 
79  subroutine xv_particles_center_gnuplot(plot_name, &
80  x, v, xmin, xmax, vmin, vmax, iplot, time)
81 
82  character(len=*), intent(in) :: plot_name
83  sll_real64, dimension(:), intent(in) :: x
84  sll_real64, dimension(:), intent(in) :: v
85  sll_real64, intent(in) :: xmin
86  sll_real64, intent(in) :: xmax
87  sll_real64, intent(in) :: vmin
88  sll_real64, intent(in) :: vmax
89  sll_int32, intent(in) :: iplot
90  sll_real64, optional :: time
91 
92  character(len=4) :: fin
93  sll_int32 :: file_id
94  sll_int32 :: nbpart
95  sll_int32 :: k
96 
97  call sll_s_int2string(iplot, fin)
98 
99  open (newunit=file_id, file=plot_name//'.gnu', position="append")
100  if (iplot <= 1) then
101  rewind(file_id)
102  write (file_id, "('set xr[',g15.3,':',g15.3,']')") xmin, xmax
103  write (file_id, "('set yr[',g15.3,':',g15.3,']')") vmin, vmax
104  end if
105  if (present(time)) then
106  write (file_id, "(A18,G10.3,A1)") "set title 'Time = ", time, "'"
107  end if
108  write (file_id, "(a)") "plot '"//plot_name//"_"//fin//".dat' w d "
109  close (file_id)
110 
111  nbpart = size(x)
112  sll_assert(nbpart == size(v))
113  open (newunit=file_id, file=plot_name//"_"//fin//'.dat')
114  do k = 1, nbpart
115  write (file_id, *) sngl(x(k)), sngl(v(k))
116  end do
117  close (file_id)
118 
119  end subroutine xv_particles_center_gnuplot
120 
121 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
122  subroutine compute_df(df, x, v, xmin, xmax, nx, vmin, vmax, nv)
123  sll_real64, dimension(:), intent(in) :: x, v
124  sll_real64 :: delta_x, delta_v, xmin, xmax, vmin, vmax
125  sll_int32, intent(in) :: nx
126  sll_int32, intent(in) :: nv
127  sll_real64, dimension(nx, nv) :: df
128  sll_real64 :: weight
129  sll_int32 :: i, j, k
130 
131  delta_x = (xmax - xmin)/(nx - 1)
132  delta_v = (vmax - vmin)/(nv - 1)
133 
134  weight = 1._f64!/(delta_x*delta_v) ! needs improvement
135 
136  df = 0.0_f64
137  do k = 1, size(x)
138  i = floor((x(k) - xmin)/delta_x) + 1
139  j = floor((v(k) - vmin)/delta_v) + 1
140  df(i, j) = df(i, j) + weight
141  end do
142 
143  end subroutine compute_df
144 
145  subroutine distribution_xv_gnuplot(plot_name, x, v, xmin, xmax, nx, &
146  vmin, vmax, nv, iplot, time)
147 
148  character(len=*), intent(in) :: plot_name
149  sll_int32, intent(in) :: nx
150  sll_int32, intent(in) :: nv
151  sll_real64, dimension(:), intent(in) :: x, v
152  sll_int32 :: iplot
153  sll_real64, dimension(nx, nv) :: df
154  sll_real64 :: time
155  sll_real64 :: delta_x, delta_v, xmin, xmax, vmin, vmax
156  character(len=4) :: fin
157 !sll_int32 :: file_id
158  sll_int32 :: error
159 
160  delta_x = (xmax - xmin)/nx
161  delta_v = (vmax - vmin)/nv
162  call sll_s_int2string(iplot, fin)
163 
164  sll_assert(size(x) == size(v))
165 
166  write (*, "(A7,G10.3)") "Time = ", time
167  call compute_df(df, x, v, xmin, xmax, nx, vmin, vmax, nv)
168 
169 !call sll_s_new_file_id(file_id, error)
170 !open( file_id, file = plot_name//'.gnu', position="append" )
171 !if ( iplot == 1 ) rewind(file_id)
172 !write(file_id,"(A18,G10.3,A1)")"set title 'Time = ",time,"'"
173 !write(file_id,*)"splot '"//plot_name//"_"//fin//".dat' w l"
174 !close(file_id)
175 
176  call sll_o_gnuplot_2d(xmin, xmax, nx, vmin, vmax, nv, df, plot_name, iplot, error)
177 
178  end subroutine distribution_xv_gnuplot
179 
180 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
181 
183  subroutine sll_s_particles_center_gnuplot_inline(x, v, xmin, xmax, ymin, ymax, time)
184  sll_real64, dimension(:), intent(in) :: x, v
185  sll_real64, intent(in) :: time
186  sll_real64 :: xmin, xmax, ymin, ymax
187  sll_int32 :: nbpart
188  sll_int32 :: k
189 
190  print *, 'set term x11 2'
191  print *, "set title 'time=", time, "'"
192  print *, "unset logscale"
193  nbpart = size(x)
194  sll_assert(nbpart == size(v))
195 
196 100 format('plot [', f7.3, ':', f7.3, '][', f7.3, ':', f7.3, '] ''-'' w d')
197 
198  write (*, 100) xmin, xmax, ymin, ymax
199  do k = 1, nbpart
200  print *, x(k), v(k)
201  end do
202  print *, 'e'
203 
204  end subroutine sll_s_particles_center_gnuplot_inline
205 
210  subroutine pq_plot_format_points3d(plot_name, x, v, iplot)
211  character(len=*), intent(in) :: plot_name
212  sll_real64, dimension(:), intent(in) :: x, v
213  sll_int32, intent(in) :: iplot
214  character(len=4) :: fin
215  sll_int32 :: file_id, error
216  sll_int32 :: nbpart
217  sll_int32 :: k
218 
219  call sll_s_int2string(iplot, fin)
220 
221  call sll_s_ascii_file_create(plot_name//"_"//fin//".3D", file_id, error)
222  nbpart = size(x)
223  sll_assert(size(v) == nbpart)
224  write (file_id, "(a)") 'x v zero zero'
225  do k = 1, nbpart
226  write (file_id, "(4e15.3)") x(k), v(k), 0., 0.
227  end do
228  close (file_id)
229 
230  end subroutine pq_plot_format_points3d
231 
235  subroutine pqr_plot_format_points3d(plot_name, x, y, z, iplot)
236  character(len=*), intent(in) :: plot_name
237  sll_real64, dimension(:), intent(in) :: x, y, z
238  sll_int32, intent(in) :: iplot
239  character(len=4) :: fin
240  sll_int32 :: file_id, error
241  sll_int32 :: nbpart
242  sll_int32 :: k
243 
244  call sll_s_int2string(iplot, fin)
245 
246  call sll_s_ascii_file_create(plot_name//"_"//fin//".3D", file_id, error)
247  nbpart = size(x)
248  sll_assert(size(y) == nbpart)
249  sll_assert(size(z) == nbpart)
250  write (file_id, "(a)") 'x y z zero'
251  do k = 1, nbpart
252  write (file_id, "(4e15.3)") x(k), y(k), z(k), 0.
253  end do
254  close (file_id)
255 
256  end subroutine pqr_plot_format_points3d
257 
262  subroutine pqrs_plot_format_points3d(plot_name, p, q, r, s, iplot)
263  character(len=*), intent(in) :: plot_name
264  sll_real64, dimension(:), intent(in) :: p, q, r, s
265  sll_int32, intent(in) :: iplot
266  character(len=4) :: fin
267  sll_int32 :: file_id, error
268  sll_int32 :: nbpart
269  sll_int32 :: k
270 
271  call sll_s_int2string(iplot, fin)
272 
273  call sll_s_ascii_file_create(plot_name//"_"//fin//".3D", file_id, error)
274  nbpart = size(p)
275  sll_assert(size(q) == nbpart)
276  sll_assert(size(r) == nbpart)
277  sll_assert(size(s) == nbpart)
278  write (file_id, "(a)") 'x y z val'
279  do k = 1, nbpart
280  write (file_id, "(4e15.3)") p(k), q(k), r(k), s(k)
281  end do
282  close (file_id)
283 
284  end subroutine pqrs_plot_format_points3d
285 
286 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
287 
292  subroutine sll_s_plot_format_xmdv(plot_name, x, v, iplot, xmin, xmax, vmin, vmax)
293  character(len=*), intent(in) :: plot_name
294  sll_real64, dimension(:), intent(in) :: x
295  sll_real64, dimension(:), intent(in) :: v
296  sll_int32, intent(in) :: iplot
297  character(len=4) :: fin
298  sll_int32 :: file_id, error
299  sll_real64 :: xmin, xmax, vmin, vmax
300  sll_int32 :: nbpart
301  sll_int32 :: k
302 
303  nbpart = size(x)
304 
305  sll_assert(size(v) == nbpart)
306  call sll_s_int2string(iplot, fin)
307 
308  call sll_s_ascii_file_create(plot_name//"_"//fin//".okc", file_id, error)
309 
310  write (file_id, "(2i7)") 2, nbpart, 1 ! two characteristics, the number 1 is unused
311  write (file_id, "(a)") 'x'
312  write (file_id, "(a)") 'v'
313  write (file_id, "(2f8.3,1x,i7)") xmin, xmax, 1
314  write (file_id, "(2f8.3,1x,i7)") vmin, vmax, 1
315  do k = 1, nbpart
316  write (file_id, "(2e15.3)") x(k), v(k)
317  end do
318  close (file_id)
319 
320  end subroutine sll_s_plot_format_xmdv
321 
322 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
323 
326  subroutine sll_s_distribution_xdmf(plot_name, x, v, w, &
327  xmin, xmax, nx, &
328  vmin, vmax, nv, iplot)
329  character(len=*), intent(in) :: plot_name
330  sll_real64, dimension(:), intent(in) :: x
331  sll_real64, dimension(:), intent(in) :: v
332  sll_real64, dimension(:), intent(in) :: w
333  sll_int32, intent(in) :: nx
334  sll_int32, intent(in) :: nv
335  sll_int32 :: iplot
336  sll_real64, dimension(nx, nv) :: df
337  sll_real64 :: delta_x, delta_v, xmin, xmax, vmin, vmax
338  character(len=4) :: fin
339 
340  call sll_s_int2string(iplot, fin)
341 
342  call sll_s_compute_df_cic(x, v, w, xmin, xmax, nx, vmin, vmax, nv, df)
343  delta_x = (xmax - xmin)/(nx - 1)
344  delta_v = (vmax - vmin)/(nv - 1)
345  call sll_s_xdmf_corect2d_nodes(plot_name//'_'//fin, df, "density", xmin, delta_x, vmin, delta_v)
346 
347  end subroutine sll_s_distribution_xdmf
348 
351  subroutine compute_df_ngp(xp, yp, wp, xmin, xmax, nx, ymin, ymax, ny, df)
352  sll_real64, dimension(:), intent(in) :: xp, yp, wp
353  sll_real64, intent(in) :: xmin, xmax, ymin, ymax
354  sll_int32 :: ip, jp, kp
355  sll_real64 :: xt, yt
356  sll_int32 :: nx, ny
357  sll_real64, dimension(nx, ny), intent(out) :: df
358  sll_int32 :: nbpart
359 
360  nbpart = size(xp)
361  sll_assert(nbpart == size(yp))
362  sll_assert(nbpart == size(wp))
363  df = 0.0_f64
364 
365  do kp = 1, nbpart
366 
367  xt = (xp(kp) - xmin)/(xmax - xmin)*nx
368  yt = (yp(kp) - ymin)/(ymax - ymin)*ny
369 
370  ip = floor(xt)
371  jp = floor(yt)
372 
373  sll_assert(ip >= 1 .and. ip <= nx .and. jp >= 1 .and. jp <= ny)
374 
375  df(ip, jp) = df(ip, jp) + wp(kp)
376 
377  end do
378 
379  end subroutine compute_df_ngp
380 
382 !Cell)
383  subroutine sll_s_compute_df_cic(xp, yp, wp, xmin, xmax, nx, ymin, ymax, ny, df)
384  sll_real64, dimension(:), intent(in) :: xp, yp, wp
385  sll_real64, intent(in) :: xmin, xmax, ymin, ymax
386  sll_int32 :: ip, jp, kp
387  sll_real64 :: a1, a2, a3, a4, xt, yt
388  sll_int32 :: nx, ny
389  sll_real64, dimension(nx, ny), intent(out) :: df
390  sll_int32 :: nbpart
391 
392  nbpart = size(xp)
393  sll_assert(nbpart == size(yp))
394  sll_assert(nbpart == size(wp))
395  df = 0.0_f64
396 
397  do kp = 1, nbpart
398 
399  xt = (xp(kp) - xmin)/(xmax - xmin)*nx
400  yt = (yp(kp) - ymin)/(ymax - ymin)*ny
401 
402  ip = floor(xt)
403  jp = floor(yt)
404 
405  a1 = (ip + 1 - xt)*(jp + 1 - yt)
406  a2 = (xt - ip)*(jp + 1 - yt)
407  a3 = (ip + 1 - xt)*(yt - jp)
408  a4 = (xt - ip)*(yt - jp)
409 
410  ip = ip + 1
411  jp = jp + 1
412  if (ip == nx + 1) ip = nx
413  if (jp == ny + 1) jp = ny
414 !print *,ip, xt, jp, floor(yt), nx, ny
415  sll_assert(ip > 0 .and. ip <= nx .and. jp > 0 .and. jp <= ny)
416 
417 ! df(ip ,jp )=df(ip ,jp )+a1*wp(kp)
418 ! df(ip+1,jp )=df(ip+1,jp )+a2*wp(kp)
419 ! df(ip ,jp+1)=df(ip ,jp+1)+a3*wp(kp)
420 ! df(ip+1,jp+1)=df(ip+1,jp+1)+a4*wp(kp)
421 
422  if (jp > 0 .and. jp <= ny) then
423  if (ip > 0 .and. ip <= nx) df(ip, jp) = df(ip, jp) + a1*wp(kp)
424  if (ip + 1 > 0 .and. ip + 1 <= nx) df(ip + 1, jp) = df(ip + 1, jp) + a2*wp(kp)
425  end if
426  if (jp + 1 > 0 .and. jp + 1 <= ny) then
427  if (ip > 0 .and. ip <= nx) df(ip, jp + 1) = df(ip, jp + 1) + a3*wp(kp)
428  if (ip + 1 > 0 .and. ip + 1 <= nx) df(ip + 1, jp + 1) = df(ip + 1, jp + 1) + a4*wp(kp)
429  end if
430  end do
431 
432  end subroutine sll_s_compute_df_cic
433 
434  subroutine compute_df_m4(xp, yp, wp, xmin, xmax, nx, ymin, ymax, ny, df)
435  sll_real64, intent(in) :: xp(:)
436  sll_real64, intent(in) :: yp(:)
437  sll_real64, intent(in) :: wp(:)
438  sll_real64, intent(in) :: xmin
439  sll_real64, intent(in) :: xmax
440  sll_int32, intent(in) :: nx
441  sll_real64, intent(in) :: ymin
442  sll_real64, intent(in) :: ymax
443  sll_int32, intent(in) :: ny
444  sll_real64, intent(out) :: df(nx, ny)
445 
446  sll_int32 :: ip, jp, kp
447  sll_real64 :: xt, x, cx, cm1x, cm2x, cp1x, cp2x
448  sll_real64 :: yt, y, cy, cm1y, cm2y, cp1y, cp2y
449 
450 ! Initialize output array to zero
451  df = 0.0_f64
452 
453  do kp = 1, size(xp)
454 
455  xt = (xp(kp) - xmin)/(xmax - xmin)*nx
456  yt = (yp(kp) - ymin)/(ymax - ymin)*ny
457 
458  ip = floor(xt)
459  jp = floor(yt)
460 
461  x = xt - ip
462  y = yt - jp
463 
464  cm2x = f_m4(2.+x)
465  cp2x = f_m4(2.-x)
466  cm1x = f_m4(1.+x)
467  cp1x = f_m4(1.-x)
468  cx = f_m4(x)
469  cy = f_m4(y)
470  cm2y = f_m4(2.+y)
471  cp2y = f_m4(2.-y)
472  cm1y = f_m4(1.+y)
473  cp1y = f_m4(1.-y)
474 
475  if (ip > 2 .and. jp > 2 .and. ip < nx - 1 .and. jp < ny - 1) then
476 
477  df(ip - 2, jp - 2) = df(ip - 2, jp - 2) + cm2x*cm2y*wp(kp)
478  df(ip - 2, jp - 1) = df(ip - 2, jp - 1) + cm2x*cm1y*wp(kp)
479  df(ip - 2, jp) = df(ip - 2, jp) + cm2x*cy*wp(kp)
480  df(ip - 2, jp + 1) = df(ip - 2, jp + 1) + cm2x*cp1y*wp(kp)
481  df(ip - 2, jp + 2) = df(ip - 2, jp + 2) + cm2x*cp2y*wp(kp)
482 
483  df(ip - 1, jp - 2) = df(ip - 1, jp - 2) + cm1x*cm2y*wp(kp)
484  df(ip - 1, jp - 1) = df(ip - 1, jp - 1) + cm1x*cm1y*wp(kp)
485  df(ip - 1, jp) = df(ip - 1, jp) + cm1x*cy*wp(kp)
486  df(ip - 1, jp + 1) = df(ip - 1, jp + 1) + cm1x*cp1y*wp(kp)
487  df(ip - 1, jp + 2) = df(ip - 1, jp + 2) + cm1x*cp2y*wp(kp)
488 
489  df(ip, jp - 2) = df(ip, jp - 2) + cx*cm2y*wp(kp)
490  df(ip, jp - 1) = df(ip, jp - 1) + cx*cm1y*wp(kp)
491  df(ip, jp) = df(ip, jp) + cx*cy*wp(kp)
492  df(ip, jp + 1) = df(ip, jp + 1) + cx*cp1y*wp(kp)
493  df(ip, jp + 2) = df(ip, jp + 2) + cx*cp2y*wp(kp)
494 
495  df(ip + 1, jp - 2) = df(ip + 1, jp - 2) + cp1x*cm2y*wp(kp)
496  df(ip + 1, jp - 1) = df(ip + 1, jp - 1) + cp1x*cm1y*wp(kp)
497  df(ip + 1, jp) = df(ip + 1, jp) + cp1x*cy*wp(kp)
498  df(ip + 1, jp + 1) = df(ip + 1, jp + 1) + cp1x*cp1y*wp(kp)
499  df(ip + 1, jp + 2) = df(ip + 1, jp + 2) + cp1x*cp2y*wp(kp)
500 
501  df(ip + 2, jp - 2) = df(ip + 2, jp - 2) + cp2x*cm2y*wp(kp)
502  df(ip + 2, jp - 1) = df(ip + 2, jp - 1) + cp2x*cm1y*wp(kp)
503  df(ip + 2, jp) = df(ip + 2, jp) + cp2x*cy*wp(kp)
504  df(ip + 2, jp + 1) = df(ip + 2, jp + 1) + cp2x*cp1y*wp(kp)
505  df(ip + 2, jp + 2) = df(ip + 2, jp + 2) + cp2x*cp2y*wp(kp)
506 
507  end if
508 
509  end do
510 
511  end subroutine compute_df_m4
512 
514  sll_real64 function f_m4(x)
515  sll_real64, intent(in) :: x
516  if (x .gt. 2.) then
517  f_m4 = 0._f64
518  else if (x .ge. 1. .and. x .le. 2.) then
519  f_m4 = 0.5*(2.-x)**2*(1.-x)
520  else if (x .le. 1.) then
521  f_m4 = 1.-2.5*x**2 + 1.5*(abs(x))**3
522  end if
523 
524  return
525  end function f_m4
526 
528  subroutine sll_s_distribution_m4_gnuplot(plot_name, x, v, w, &
529  xmin, xmax, nx, &
530  vmin, vmax, nv, iplot)
531  character(len=*), intent(in) :: plot_name
532  sll_real64, intent(in) :: x(:)
533  sll_real64, intent(in) :: v(:)
534  sll_real64, intent(in) :: w(:)
535  sll_real64, intent(in) :: xmin
536  sll_real64, intent(in) :: xmax
537  sll_int32, intent(in) :: nx
538  sll_real64, intent(in) :: vmin
539  sll_real64, intent(in) :: vmax
540  sll_int32, intent(in) :: nv
541  sll_int32, intent(in) :: iplot
542 
543  sll_int32 :: error
544  sll_real64 :: df(nx, nv)
545  character(len=4) :: fin
546 
547  call sll_s_int2string(iplot, fin)
548 
549  call compute_df_m4(x, v, w, xmin, xmax, nx, vmin, vmax, nv, df)
550 
551  call sll_o_gnuplot_2d(xmin, xmax, nx, vmin, vmax, nv, df, plot_name, iplot, error)
552 
553  end subroutine sll_s_distribution_m4_gnuplot
554 
555 #define FONCTION1( X ) (0.75_f64-(X)*(X))
556 
557 #define FONCTION2( X ) (0.5_f64 * ( 1.5_f64 - (X) )**2)
558 
559 #define BSPLINES(X,Y) \
560  c_1x = fonction2(1 + x); \
561  c1x = fonction1(x); \
562  c2x = fonction2(1 - x); \
563  c_1y = fonction2(1 + y); \
564  c1y = fonction1(y); \
565  c2y = fonction2(1 - y)
566 
567  subroutine compute_df_tsc(xp, yp, wp, xmin, xmax, nx, ymin, ymax, ny, df)
568  sll_real64, dimension(:), intent(in) :: xp, yp, wp
569  sll_real64, intent(in) :: xmin, xmax, ymin, ymax
570  sll_int32 :: ip, jp, kp
571  sll_real64 :: xt, yt
572  sll_int32 :: nx, ny
573  sll_real64, dimension(nx, ny), intent(out) :: df
574  sll_real64 :: x, c1x, c_1x, c2x
575  sll_real64 :: y, c1y, c_1y, c2y
576 
577  df = 0.0_f64
578 
579  do kp = 1, size(xp)
580 
581  xt = (xp(kp) - xmin)/(xmax - xmin)*nx
582  yt = (yp(kp) - ymin)/(ymax - ymin)*ny
583 
584  ip = floor(xt)
585  jp = floor(yt)
586 
587  x = xt - ip
588  y = yt - jp
589 
590  bsplines(x, y)
591 
592  if (ip > 1 .and. ip < nx .and. jp > 1 .and. jp < ny) then
593 
594  df(ip, jp) = df(ip, jp) + c1x*c1y*wp(kp)
595  df(ip, jp - 1) = df(ip, jp - 1) + c1x*c_1y*wp(kp)
596  df(ip, jp + 1) = df(ip, jp + 1) + c1x*c2y*wp(kp)
597  df(ip - 1, jp) = df(ip - 1, jp) + c_1x*c1y*wp(kp)
598  df(ip - 1, jp - 1) = df(ip - 1, jp - 1) + c_1x*c_1y*wp(kp)
599  df(ip - 1, jp + 1) = df(ip - 1, jp + 1) + c_1x*c2y*wp(kp)
600  df(ip + 1, jp) = df(ip + 1, jp) + c2x*c1y*wp(kp)
601  df(ip + 1, jp - 1) = df(ip + 1, jp - 1) + c2x*c_1y*wp(kp)
602  df(ip + 1, jp + 1) = df(ip + 1, jp + 1) + c2x*c2y*wp(kp)
603 
604  end if
605 
606  end do
607 
608  end subroutine compute_df_tsc
609 
611  subroutine sll_s_distribution_tsc_gnuplot(plot_name, x, v, w, &
612  xmin, xmax, nx, &
613  vmin, vmax, nv, iplot)
614  character(len=*), intent(in) :: plot_name
615  sll_real64, dimension(:), intent(in) :: x
616  sll_real64, dimension(:), intent(in) :: v
617  sll_real64, dimension(:), intent(in) :: w
618  sll_int32, intent(in) :: nx
619  sll_int32, intent(in) :: nv
620  sll_int32 :: iplot, error
621  sll_real64, dimension(nx, nv) :: df
622  sll_real64 :: xmin, xmax, vmin, vmax
623  character(len=4) :: fin
624 
625  call sll_s_int2string(iplot, fin)
626 
627  call compute_df_tsc(x, v, w, xmin, xmax, nx, vmin, vmax, nv, df)
628 
629  call sll_o_gnuplot_2d(xmin, xmax, nx, vmin, vmax, nv, df, plot_name, iplot, error)
630 
631  end subroutine sll_s_distribution_tsc_gnuplot
632 
633 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
634 
638  subroutine sll_s_energies_electrostatic_gnuplot_inline(kinetic_energy, &
639  electrostatic_energy, impulse, timestepwidth)
640  sll_real64, dimension(:), intent(in) :: kinetic_energy, electrostatic_energy, impulse
641  sll_int32 :: timesteps, k
642  sll_real64, dimension(size(kinetic_energy)) :: total_energy
643  sll_real64, intent(in) :: timestepwidth
644 
645  timesteps = size(kinetic_energy)
646  sll_assert(timesteps == size(electrostatic_energy))
647  sll_assert(timesteps == size(total_energy))
648 
649  total_energy = electrostatic_energy + kinetic_energy
650 
651  if (timesteps > 2) then
652  print *, " "
653 !!Plot Energy errors
654  print *, "set term x11 3"
655  print *, "set logscale y"
656  print *, "set title 'Relative Energies'"
657  print *, "plot '-' title 'kinetic' with lines,\ "
658  if (maxval(electrostatic_energy) /= 0.0_f64) print *, "'-' title 'electrostatic' with lines,\ "
659  print *, "'-' title 'total' with lines;"
660  do k = 2, timesteps
661  print *, (k - 1)*timestepwidth, kinetic_energy(k)/maxval(kinetic_energy)
662  end do
663  print *, "e"
664  if (maxval(electrostatic_energy) /= 0.0_f64) then
665  do k = 2, timesteps
666  print *, (k - 1)*timestepwidth, electrostatic_energy(k)/maxval(electrostatic_energy)
667  end do
668  end if
669  print *, "e"
670  do k = 2, timesteps
671  print *, (k - 1)*timestepwidth, total_energy(k)/maxval(total_energy)
672  end do
673  print *, "e"
674  print *, " "
675 
676 !Plot Energy errors
677  print *, "set term x11 4"
678  print *, "set autoscale y"
679  print *, "set title 'Total Energy Error'"
680  print *, "set logscale y"
681 
682  print *, "plot '-' using 1:2 title 'Relative Error' with lines"
683  do k = 2, timesteps
684  print *, (k - 1)*timestepwidth, abs(((total_energy(k) - total_energy(1))/total_energy(1)))
685  end do
686  print *, "e"
687  print *, " "
688 
689 !!plot Impulse Error
690  print *, "set term x11 6"
691  print *, "set autoscale y"
692  print *, "set title 'Relative Impulse Error'"
693  print *, "set logscale y"
694  print *, "plot '-' using 1:2 title 'Relative Impulse Error' with lines"
695  do k = 2, timesteps
696  if (impulse(1) /= 0) then
697  print *, (k - 1)*timestepwidth, abs(impulse(k) - impulse(1))/abs(impulse(1))
698  else if (maxval(abs(impulse)) /= 0.0_f64) then
699  print *, (k - 1)*timestepwidth, abs(impulse(k) - impulse(1))/maxval(abs(impulse))
700  else
701  print *, (k - 1)*timestepwidth, abs(impulse(k) - impulse(1))
702  end if
703  end do
704  print *, "e"
705 
706 !!!plot Electrostatic energy to see landau damping
707 !print*, "set term x11 5"
708 !print*, "set autoscale y"
709 !print*, "set title 'Relative Electrostatic Energy'"
710 !print*, "set logscale y"
711 !print*, "plot '-' using 1:2 title 'Electrostatic Energy' with lines"
712 !do k = 2, timesteps
713 ! print*, k, electrostatic_energy(k)
714 !end do
715 !print*, "e"
716 
717  end if
718  end subroutine sll_s_energies_electrostatic_gnuplot_inline
719 
723  subroutine sll_s_electricpotential_gnuplot_inline(values, eval_knots)
724  sll_real64, dimension(:), intent(in) :: values, eval_knots
725  sll_int32 :: nknots, k
726 
727  nknots = size(eval_knots)
728  sll_assert(nknots == size(eval_knots))
729 
730 !Plot Energy errors
731  print *, "set term x11 7"
732  print *, "set autoscale y"
733  print *, "set title 'Electric Potential'"
734  print *, "unset logscale y"
735  print *, "plot '-' using 1:2 title 'Potential' with lines"
736  do k = 1, nknots
737  print *, eval_knots(k), values(k)
738  end do
739  print *, "e"
740 
741  end subroutine sll_s_electricpotential_gnuplot_inline
742 
743 end module sll_m_pic_visu
sll_m_gnuplot::sll_o_gnuplot_2d
Write file for gnuplot to display 2d field.
Definition: sll_m_gnuplot.F90:55
sll_m_pic_visu::sll_o_distribution_gnuplot
plot particles distribution with gnuplot
Definition: sll_m_pic_visu.F90:65
sll_m_pic_visu::sll_o_particles_center_gnuplot
plot particles centers with gnuplot
Definition: sll_m_pic_visu.F90:60
sll_m_pic_visu::sll_o_plot_format_points3d
write point3d file to plot particles characteristics
Definition: sll_m_pic_visu.F90:70
sll_m_ascii_io
Module that contains routines to write data in ASCII format file.
Definition: sll_m_ascii_io.F90:25
sll_m_ascii_io::sll_s_ascii_file_create
subroutine, public sll_s_ascii_file_create(file_name, file_id, error)
Create ASCII file.
Definition: sll_m_ascii_io.F90:56
sll_m_gnuplot
Implements the functions to write data file plotable by GNUplot.
Definition: sll_m_gnuplot.F90:27
sll_m_pic_visu
This module provides some routines for plotting during PIC simulations.
Definition: sll_m_pic_visu.F90:22
sll_m_pic_visu::xv_particles_center_gnuplot
subroutine xv_particles_center_gnuplot(plot_name, x, v, xmin, xmax, vmin, vmax, iplot, time)
To plot particles run : gnuplot -persitant plot_name.gnu.
Definition: sll_m_pic_visu.F90:81
sll_m_pic_visu::pqrs_plot_format_points3d
subroutine pqrs_plot_format_points3d(plot_name, p, q, r, s, iplot)
point3D format http://www.visitusers.org/index.php?title=Reading_point_data This format is designed t...
Definition: sll_m_pic_visu.F90:263
sll_m_pic_visu::sll_s_distribution_m4_gnuplot
subroutine, public sll_s_distribution_m4_gnuplot(plot_name, x, v, w, xmin, xmax, nx, vmin, vmax, nv, iplot)
GNUplot readable output for particles density.
Definition: sll_m_pic_visu.F90:531
sll_m_pic_visu::sll_s_energies_electrostatic_gnuplot_inline
subroutine, public sll_s_energies_electrostatic_gnuplot_inline(kinetic_energy, electrostatic_energy, impulse, timestepwidth)
Call this function inline with | gnuplot Plots different energies for Electrostatic Particle-in-Cell ...
Definition: sll_m_pic_visu.F90:640
sll_m_pic_visu::sll_s_distribution_xdmf
subroutine, public sll_s_distribution_xdmf(plot_name, x, v, w, xmin, xmax, nx, vmin, vmax, nv, iplot)
VisIt readable output for particles density Data file format could be XML, HDF5 or Binary (not fully ...
Definition: sll_m_pic_visu.F90:329
sll_m_pic_visu::distribution_xv_gnuplot
subroutine distribution_xv_gnuplot(plot_name, x, v, xmin, xmax, nx, vmin, vmax, nv, iplot, time)
Definition: sll_m_pic_visu.F90:147
sll_m_pic_visu::sll_s_particles_center_gnuplot_inline
subroutine, public sll_s_particles_center_gnuplot_inline(x, v, xmin, xmax, ymin, ymax, time)
Call this function inline with | gnuplot.
Definition: sll_m_pic_visu.F90:184
sll_m_pic_visu::compute_df_m4
subroutine compute_df_m4(xp, yp, wp, xmin, xmax, nx, ymin, ymax, ny, df)
Definition: sll_m_pic_visu.F90:435
sll_m_pic_visu::sll_s_plot_format_xmdv
subroutine, public sll_s_plot_format_xmdv(plot_name, x, v, iplot, xmin, xmax, vmin, vmax)
xmdv format http://davis.wpi.edu/xmdv/fileformats.html This format is designed for particles with a l...
Definition: sll_m_pic_visu.F90:293
sll_m_pic_visu::pqr_plot_format_points3d
subroutine pqr_plot_format_points3d(plot_name, x, y, z, iplot)
point3D format http://www.visitusers.org/index.php?title=Reading_point_data This format is designed t...
Definition: sll_m_pic_visu.F90:236
sll_m_pic_visu::compute_df_ngp
subroutine compute_df_ngp(xp, yp, wp, xmin, xmax, nx, ymin, ymax, ny, df)
Compute grid field from particles distribution with the NGP scheme (Nearest Grid Point)
Definition: sll_m_pic_visu.F90:352
sll_m_pic_visu::sll_s_distribution_tsc_gnuplot
subroutine, public sll_s_distribution_tsc_gnuplot(plot_name, x, v, w, xmin, xmax, nx, vmin, vmax, nv, iplot)
GNUplot readable output for particles density.
Definition: sll_m_pic_visu.F90:614
sll_m_pic_visu::sll_s_electricpotential_gnuplot_inline
subroutine, public sll_s_electricpotential_gnuplot_inline(values, eval_knots)
Call this function inline with | gnuplot Plots different energies for Electrostatic Particle-in-Cell ...
Definition: sll_m_pic_visu.F90:724
sll_m_pic_visu::compute_df
subroutine compute_df(df, x, v, xmin, xmax, nx, vmin, vmax, nv)
Definition: sll_m_pic_visu.F90:123
sll_m_pic_visu::pq_plot_format_points3d
subroutine pq_plot_format_points3d(plot_name, x, v, iplot)
point3D format http://www.visitusers.org/index.php?title=Reading_point_data This format is designed t...
Definition: sll_m_pic_visu.F90:211
sll_m_pic_visu::f_m4
real(kind=f64) function f_m4(x)
M4 function from Monhagan (SPH method)
Definition: sll_m_pic_visu.F90:515
sll_m_pic_visu::sll_s_compute_df_cic
subroutine, public sll_s_compute_df_cic(xp, yp, wp, xmin, xmax, nx, ymin, ymax, ny, df)
Compute grid field from particles distribution with the CIC scheme (Cloud In.
Definition: sll_m_pic_visu.F90:384
sll_m_pic_visu::compute_df_tsc
subroutine compute_df_tsc(xp, yp, wp, xmin, xmax, nx, ymin, ymax, ny, df)
Definition: sll_m_pic_visu.F90:568
sll_m_utilities
Some common numerical utilities.
Definition: sll_m_utilities.F90:20
sll_m_utilities::sll_s_int2string
subroutine, public sll_s_int2string(istep, cstep)
Convert an integer < 9999 to a 4 characters string.
Definition: sll_m_utilities.F90:170
sll_m_utilities::sll_s_new_file_id
subroutine, public sll_s_new_file_id(file_id, error)
Get a file unit number free before creating a file.
Definition: sll_m_utilities.F90:191
sll_m_xdmf
Implements the functions to write xdmf file plotable by VisIt.
Definition: sll_m_xdmf.F90:24
sll_m_xdmf::sll_s_xdmf_corect2d_nodes
subroutine, public sll_s_xdmf_corect2d_nodes(file_name, array, array_name, eta1_min, delta_eta1, eta2_min, delta_eta2, file_format, iplot, time)
Subroutine to write a 2D array in xdmf format The field is describe on a cartesian mesh Axis are perp...
Definition: sll_m_xdmf.F90:265
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:41 for Semi-Lagrangian Library by   doxygen 1.9.1