Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_maxwell_3d_trafo.F90
Go to the documentation of this file.
1 
8 
9 module sll_m_maxwell_3d_trafo
10  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
11 #include "sll_assert.h"
12 #include "sll_memory.h"
13 #include "sll_working_precision.h"
14 
15  use sll_m_collective, only: &
16  sll_f_get_collective_rank, &
17  sll_v_world_collective
18 
19  use sll_m_gauss_legendre_integration, only: &
20  sll_f_gauss_legendre_points_and_weights
21 
22  use sll_m_linear_operator_block, only : &
23  sll_t_linear_operator_block
24 
25  use sll_m_linear_operator_penalized, only : &
26  sll_t_linear_operator_penalized
27 
28  use sll_m_linear_operator_poisson_3d, only : &
29  sll_t_linear_operator_poisson_3d
30 
31  use sll_m_linear_operator_schur_eb_3d, only : &
32  sll_t_linear_operator_schur_eb_3d
33 
34  use sll_m_linear_solver_cg, only : &
35  sll_t_linear_solver_cg
36 
37  use sll_m_low_level_bsplines, only: &
38  sll_s_uniform_bsplines_eval_basis
39 
40  use sll_m_matrix_csr, only: &
41  sll_t_matrix_csr
42 
43  use sll_m_mapping_3d, only: &
44  sll_t_mapping_3d
45 
46  use sll_m_maxwell_3d_base, only: &
47  sll_i_function_3d_real64, &
48  sll_c_maxwell_3d_base
49 
50  use sll_m_preconditioner_fft, only : &
51  sll_t_preconditioner_fft
52 
53  use sll_m_profile_functions, only: &
54  sll_t_profile_functions
55 
56  use sll_m_spline_fem_utilities, only : &
57  sll_s_multiply_g, &
58  sll_s_multiply_gt, &
59  sll_s_multiply_c, &
60  sll_s_multiply_ct
61 
62  use sll_m_spline_fem_utilities_3d, only: &
63  sll_s_spline_fem_mass3d, &
64  sll_s_spline_fem_mixedmass3d
65 
66  implicit none
67 
68  public :: &
69  sll_t_maxwell_3d_trafo
70 
71  private
72  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
73 
74  type, extends(sll_c_maxwell_3d_base) :: sll_t_maxwell_3d_trafo
75 
76  sll_real64, allocatable :: work0(:)
77  sll_real64, allocatable :: work(:)
78  sll_real64, allocatable :: work2(:)
79  type(sll_t_matrix_csr) :: mass0
80  type(sll_t_matrix_csr) :: mass1(3,3)
81  type(sll_t_matrix_csr) :: mass2(3,3)
82  type(sll_t_matrix_csr) :: mixed_mass(6)
83  type(sll_t_linear_operator_block) :: mass1_operator
84  type(sll_t_linear_operator_block) :: mass2_operator
85  type(sll_t_linear_solver_cg) :: mass0_solver
86  type(sll_t_linear_solver_cg) :: mass1_solver
87  type(sll_t_linear_solver_cg) :: mass2_solver
88  type(sll_t_preconditioner_fft) :: preconditioner_fft
89 
90  type(sll_t_linear_operator_poisson_3d) :: poisson_matrix
91  type(sll_t_linear_operator_penalized) :: poisson_operator
92  type(sll_t_linear_solver_cg) :: poisson_solver
93  type( sll_t_linear_operator_schur_eb_3d ) :: linear_op_schur_eb
94  type( sll_t_linear_solver_cg ) :: linear_solver_schur_eb
95  type(sll_t_mapping_3d), pointer :: map
96 
97  logical :: adiabatic_electrons = .false.
98 
99  contains
100 
101  procedure :: &
102  compute_e_from_b => sll_s_compute_e_from_b_3d_trafo
103  procedure :: &
104  compute_b_from_e => sll_s_compute_b_from_e_3d_trafo
105  procedure :: &
106  compute_curl_part => sll_s_compute_curl_part_3d_trafo
107  procedure :: &
108  compute_e_from_rho => sll_s_compute_e_from_rho_3d_trafo
109  procedure :: &
110  compute_rho_from_e => sll_s_compute_rho_from_e_3d_trafo
111  procedure :: &
112  compute_e_from_j => sll_s_compute_e_from_j_3d_trafo
113  procedure :: &
114  compute_phi_from_rho => sll_s_compute_phi_from_rho_3d_trafo
115  procedure :: &
116  compute_phi_from_j => sll_s_compute_phi_from_j_3d_trafo
117  procedure :: &
118  compute_rhs_from_function => sll_s_compute_rhs_trafo
119  procedure :: &
120  l2projection => l2projection_3d_trafo
121  procedure :: &
122  l2norm_squared => l2norm_squared_3d_trafo
123  procedure :: &
124  inner_product => inner_product_3d_trafo
125  procedure :: &
126  init => init_3d_trafo
127  procedure :: &
128  init_from_file => init_from_file_3d_trafo
129  procedure :: &
130  free => free_3d_trafo
131  procedure :: &
132  multiply_g
133  procedure :: &
134  multiply_gt
135  procedure :: &
136  multiply_c
137  procedure :: &
138  multiply_ct
139  procedure :: &
140  multiply_mass => multiply_mass_trafo
141  procedure :: &
142  multiply_mass_inverse => multiply_mass_inverse_trafo
143 
144  end type sll_t_maxwell_3d_trafo
145 
146 contains
147 
148 
150  subroutine sll_s_compute_e_from_b_3d_trafo( self, delta_t, field_in, field_out )
151  class(sll_t_maxwell_3d_trafo) :: self
152  sll_real64, intent( in ) :: delta_t
153  sll_real64, intent( in ) :: field_in(:)
154  sll_real64, intent( inout ) :: field_out(:)
155 
156  call multiply_mass_trafo( self, [2], field_in, self%work )
157 
158  call multiply_ct( self, self%work, self%work2 )
159 
160  call self%mass1_solver%solve( self%work2, self%work )
161  ! Update e from self value
162  field_out = field_out + delta_t*self%work
163 
164  end subroutine sll_s_compute_e_from_b_3d_trafo
165 
166 
168  subroutine sll_s_compute_b_from_e_3d_trafo( self, delta_t, field_in, field_out )
169  class(sll_t_maxwell_3d_trafo) :: self
170  sll_real64, intent( in ) :: delta_t
171  sll_real64, intent( in ) :: field_in(:)
172  sll_real64, intent( inout ) :: field_out(:)
173 
174  call multiply_c( self, field_in, self%work )
175  ! Update b from self value
176  field_out = field_out - delta_t * self%work
177 
178  end subroutine sll_s_compute_b_from_e_3d_trafo
179 
180 
182  subroutine sll_s_compute_curl_part_3d_trafo( self, delta_t, efield, bfield, betar )
183  class(sll_t_maxwell_3d_trafo) :: self
184  sll_real64, intent( in ) :: delta_t
185  sll_real64, intent( inout ) :: efield(:)
186  sll_real64, intent( inout ) :: bfield(:)
187  sll_real64, optional :: betar
188  !local variables
189  sll_real64 :: factor
190 
191  if( present(betar) ) then
192  factor = betar
193  else
194  factor = 1._f64
195  end if
196 
197  ! Compute C^T M2 b
198  call self%mass2_operator%dot( bfield, self%work )
199  call self%multiply_ct( self%work, self%work2 )
200 
201  self%linear_op_schur_eb%sign = -delta_t**2*factor*0.25_f64
202  call self%linear_op_schur_eb%dot( efield, self%work )
203  self%work = self%work + delta_t*factor*self%work2
204 
205  ! Save efield dofs from previous time step for B field update
206  self%work2 = efield
207 
208  ! Invert Schur complement matrix
209  self%linear_op_schur_eb%sign = delta_t**2*factor*0.25_f64
210  call self%linear_solver_schur_eb%set_guess( efield )
211  call self%linear_solver_schur_eb%solve( self%work, efield)
212 
213  ! Update B field
214  self%work2 = self%work2 + efield
215  call self%compute_b_from_e( delta_t*0.5_f64, self%work2, bfield)
216 
217  end subroutine sll_s_compute_curl_part_3d_trafo
218 
219 
221  subroutine sll_s_compute_e_from_rho_3d_trafo( self, field_in, field_out )
222  class(sll_t_maxwell_3d_trafo) :: self
223  sll_real64, intent( in ) :: field_in(:)
224  sll_real64, intent( out ) :: field_out(:)
225 
226  call self%poisson_solver%solve( field_in, self%work0 )
227  call multiply_g(self, self%work0, field_out)
228  field_out = -field_out
229 
230  end subroutine sll_s_compute_e_from_rho_3d_trafo
231 
232 
234  subroutine sll_s_compute_rho_from_e_3d_trafo( self, field_in, field_out )
235  class(sll_t_maxwell_3d_trafo) :: self
236  sll_real64, intent( in ) :: field_in(:)
237  sll_real64, intent( out ) :: field_out(:)
238 
239  call multiply_mass_trafo( self, [1], field_in, self%work )
240 
241  call multiply_gt( self, self%work, field_out )
242 
243  field_out = - field_out
244 
245  end subroutine sll_s_compute_rho_from_e_3d_trafo
246 
247 
249  subroutine sll_s_compute_e_from_j_3d_trafo( self, current, E, component )
250  class(sll_t_maxwell_3d_trafo) :: self
251  sll_real64, intent( in ) :: current(:)
252  sll_real64, intent( inout ) :: e(:)
253  sll_int32, intent( in ), optional :: component
254 
255  call self%mass1_solver%solve( current, self%work )
256  e = e - self%work
257 
258  end subroutine sll_s_compute_e_from_j_3d_trafo
259 
260 
262  subroutine sll_s_compute_phi_from_rho_3d_trafo( self, field_in, field_out, efield_dofs )
263  class(sll_t_maxwell_3d_trafo) :: self
264  sll_real64, intent( in ) :: field_in(:)
265  sll_real64, intent( inout ) :: field_out(:)
266  sll_real64, intent( out ) :: efield_dofs(:)
267 
268  call self%mass0_solver%solve( field_in, field_out )
269  call self%multiply_g( field_out, efield_dofs )
270  efield_dofs = -efield_dofs
271 
272  end subroutine sll_s_compute_phi_from_rho_3d_trafo
273 
274 
276  subroutine sll_s_compute_phi_from_j_3d_trafo( self, field_in, field_out, efield_dofs )
277  class(sll_t_maxwell_3d_trafo) :: self
278  sll_real64, intent( in ) :: field_in(:)
279  sll_real64, intent( inout ) :: field_out(:)
280  sll_real64, intent( out ) :: efield_dofs(:)
281 
282  call self%multiply_gt( field_in, self%work0 )
283  call self%mass0_solver%solve( self%work0, self%work2(1:self%n_total) )
284  field_out = field_out + self%work2(1:self%n_total)
285 
286  call self%multiply_g( field_out, efield_dofs )
287  efield_dofs = -efield_dofs
288 
289  end subroutine sll_s_compute_phi_from_j_3d_trafo
290 
291 
294  subroutine sll_s_compute_rhs_trafo( self, form, component, coefs_dofs, func1, func2, func3 )
295  class(sll_t_maxwell_3d_trafo) :: self
296  sll_int32, intent( in ) :: form
297  sll_int32, intent( in ) :: component
298  sll_real64, intent( out ) :: coefs_dofs(:)
299  procedure(sll_i_function_3d_real64) :: func1
300  procedure(sll_i_function_3d_real64), optional :: func2
301  procedure(sll_i_function_3d_real64), optional :: func3
302  ! local variables
303  sll_int32 :: degree(3)
304 
305  ! Define the spline degree in the 3 dimensions, depending on form and component of the form
306  if ( form == 0 ) then
307  call rhs_zeroform( self, self%s_deg_0, func1, coefs_dofs )
308  elseif (form == 1 ) then
309  degree = self%s_deg_0
310  degree(component) = self%s_deg_1(component)
311  call rhs_oneform( self, degree, func1, func2, func3, component, coefs_dofs )
312  elseif( form == 2) then
313  degree = self%s_deg_1
314  degree(component) = self%s_deg_0(component)
315  call rhs_twoform( self, degree, func1, func2, func3, component, coefs_dofs )
316  elseif( form == 3) then
317  call rhs_threeform( self, self%s_deg_1, func1, coefs_dofs )
318  else
319  print*, 'Wrong form.'
320  end if
321 
322  end subroutine sll_s_compute_rhs_trafo
323 
324 
326  subroutine rhs_zeroform( self, deg, func, coefs_dofs )
327  class(sll_t_maxwell_3d_trafo) :: self
328  sll_int32, intent( in ) :: deg(3)
329  procedure(sll_i_function_3d_real64) :: func
330  sll_real64, intent( out ) :: coefs_dofs(:)
331  ! local variables
332  sll_int32 :: i1, i2, i3, j1, j2, j3, k1, k2, k3, q(3), counter
333  sll_real64 :: coef, jacobian, c(3)
334  sll_real64, allocatable :: xw_gauss_d1(:,:), xw_gauss_d2(:,:), xw_gauss_d3(:,:)
335  sll_real64, allocatable :: bspl_d1(:,:), bspl_d2(:,:), bspl_d3(:,:)
336  sll_real64 :: scratch(maxval(deg+1))
337 
338  ! take enough Gauss points so that projection is exact for splines of degree deg
339  q = deg+1
340  ! rescale on [0,1] for compatibility with B-splines
341  allocate(xw_gauss_d1(2, q(1)))
342  allocate(bspl_d1(q(1), deg(1)+1 ))
343  xw_gauss_d1 = sll_f_gauss_legendre_points_and_weights(q(1), 0.0_f64, 1.0_f64)
344  ! Compute bsplines at gauss_points
345  do k1 = 1, q(1)
346  call sll_s_uniform_bsplines_eval_basis(deg(1),xw_gauss_d1(1,k1), scratch(1:deg(1)+1))
347  bspl_d1(k1,:) = scratch(1:deg(1)+1)
348  end do
349 
350  allocate(xw_gauss_d2(2, q(2)))
351  allocate(bspl_d2(q(2), deg(2)+1))
352  xw_gauss_d2 = sll_f_gauss_legendre_points_and_weights(q(2), 0.0_f64, 1.0_f64)
353  ! Compute bsplines at gauss_points
354  do k2 = 1, q(2)
355  call sll_s_uniform_bsplines_eval_basis(deg(2),xw_gauss_d2(1,k2), scratch(1:deg(2)+1))
356  bspl_d2(k2,:) = scratch(1:deg(2)+1)
357  end do
358 
359  allocate(xw_gauss_d3(2,q(3)))
360  allocate(bspl_d3(q(3), deg(3)+1 ))
361  xw_gauss_d3 = sll_f_gauss_legendre_points_and_weights(q(3), 0.0_f64, 1.0_f64)
362  ! Compute bsplines at gauss_points
363  do k3 = 1, q(3)
364  call sll_s_uniform_bsplines_eval_basis(deg(3),xw_gauss_d3(1,k3), scratch(1:deg(3)+1) )
365  bspl_d3(k3,:) = scratch(1:deg(3)+1)
366  end do
367 
368 
369  counter = 1
370  ! Compute coefs_dofs = int f(\xi)N_i(\xi) J_F(\xi) d\xi
371  do i3 = 1, self%n_dofs(3)
372  do i2 = 1, self%n_dofs(2)
373  do i1 = 1, self%n_dofs(1)
374  coef = 0.0_f64
375  ! loop over support of B spline
376  do j3 = 1, deg(3)+1
377  do j2 = 1, deg(2)+1
378  do j1 = 1, deg(1)+1
379  ! loop over Gauss points
380  do k3 = 1, q(3)
381  c(3) = self%delta_x(3) * (xw_gauss_d3(1,k3) + real(i3 + j3 - 2-((i3 + j3 - 2)/self%n_dofs(3)) * self%n_dofs(3),f64))
382  do k2 = 1, q(2)
383  c(2) = self%delta_x(2) * (xw_gauss_d2(1,k2) + real(i2 + j2 - 2-((i2 + j2 - 2)/self%n_dofs(2)) * self%n_dofs(2),f64))
384  do k1 = 1, q(1)
385  c(1) = self%delta_x(1) * (xw_gauss_d1(1,k1) + real(i1 + j1 - 2-((i1 + j1 - 2)/self%n_dofs(1)) * self%n_dofs(1),f64))
386  jacobian = self%map%jacobian( c )
387  coef = coef + xw_gauss_d1(2,k1) * xw_gauss_d2(2,k2) * xw_gauss_d3(2,k3) *&
388  func(self%map%get_x(c) ) * &
389  bspl_d1(k1,deg(1)+2-j1) *&
390  bspl_d2(k2,deg(2)+2-j2) *&
391  bspl_d3(k3,deg(3)+2-j3) * abs(jacobian)
392  enddo
393  enddo
394  end do
395  end do
396  end do
397  end do
398  ! rescale by cell size
399  coefs_dofs(counter) = coef*self%volume
400  counter = counter+1
401  enddo
402  end do
403  end do
404 
405  end subroutine rhs_zeroform
406 
407 
409  subroutine rhs_oneform( self, deg, func1, func2, func3, component, coefs_dofs )
410  class(sll_t_maxwell_3d_trafo) :: self
411  sll_int32, intent( in ) :: deg(3)
412  procedure(sll_i_function_3d_real64) :: func1
413  procedure(sll_i_function_3d_real64) :: func2
414  procedure(sll_i_function_3d_real64) :: func3
415  sll_int32, intent( in ) :: component
416  sll_real64, intent( out ) :: coefs_dofs(:)
417  ! local variables
418  sll_int32 :: i1, i2, i3,j1, j2, j3, k1, k2, k3, q(3), counter
419  sll_real64 :: c(3)
420  sll_real64 :: coef, jacobian, jmatrix(3,3)
421  sll_real64, allocatable :: xw_gauss_d1(:,:), xw_gauss_d2(:,:), xw_gauss_d3(:,:)
422  sll_real64, allocatable :: bspl_d1(:,:), bspl_d2(:,:), bspl_d3(:,:)
423  sll_real64 :: scratch(maxval(deg+1))
424 
425  ! take enough Gauss points so that projection is exact for splines of degree deg
426  q = deg+1
427  ! rescale on [0,1] for compatibility with B-splines
428  allocate(xw_gauss_d1(2, q(1)))
429  allocate(bspl_d1(q(1), deg(1)+1 ))
430  xw_gauss_d1 = sll_f_gauss_legendre_points_and_weights(q(1), 0.0_f64, 1.0_f64)
431  ! Compute bsplines at gauss_points
432  do k1 = 1, q(1)
433  call sll_s_uniform_bsplines_eval_basis(deg(1),xw_gauss_d1(1,k1), scratch(1:deg(1)+1))
434  bspl_d1(k1,:) = scratch(1:deg(1)+1)
435  end do
436 
437  allocate(xw_gauss_d2(2, q(2)))
438  allocate(bspl_d2(q(2), deg(2)+1))
439  xw_gauss_d2 = sll_f_gauss_legendre_points_and_weights(q(2), 0.0_f64, 1.0_f64)
440  ! Compute bsplines at gauss_points
441  do k2 = 1, q(2)
442  call sll_s_uniform_bsplines_eval_basis(deg(2),xw_gauss_d2(1,k2), scratch(1:deg(2)+1))
443  bspl_d2(k2,:) = scratch(1:deg(2)+1)
444  end do
445 
446  allocate(xw_gauss_d3(2,q(3)))
447  allocate(bspl_d3(q(3), deg(3)+1 ))
448  xw_gauss_d3 = sll_f_gauss_legendre_points_and_weights(q(3), 0.0_f64, 1.0_f64)
449  ! Compute bsplines at gauss_points
450  do k3 = 1, q(3)
451  call sll_s_uniform_bsplines_eval_basis(deg(3),xw_gauss_d3(1,k3), scratch(1:deg(3)+1) )
452  bspl_d3(k3,:) = scratch(1:deg(3)+1)
453  end do
454 
455 
456  counter = 1
457  ! Compute coefs_dofs = int f(\xi)N_i(\xi) J_F(\xi) d\xi
458  do i3 = 1, self%n_dofs(3)
459  do i2 = 1, self%n_dofs(2)
460  do i1 = 1, self%n_dofs(1)
461  coef = 0.0_f64
462  ! loop over support of B spline
463  do j3 = 1, deg(3)+1
464  do j2 = 1, deg(2)+1
465  do j1 = 1, deg(1)+1
466  ! loop over Gauss points
467  do k3 = 1, q(3)
468  c(3) = self%delta_x(3) * (xw_gauss_d3(1,k3) + real(i3 + j3 - 2, f64))
469  do k2 = 1, q(2)
470  c(2) = self%delta_x(2) * (xw_gauss_d2(1,k2) + real(i2 + j2 - 2, f64))
471  do k1 = 1, q(1)
472  c(1) = self%delta_x(1) * (xw_gauss_d1(1,k1) + real(i1 + j1 - 2, f64))
473  jacobian = self%map%jacobian( c )
474  jmatrix = self%map%jacobian_matrix_inverse_transposed( c )
475  coef = coef + xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* xw_gauss_d3(2,k3) *&
476  ( func1( self%map%get_x(c) )*jmatrix(1,component) + &
477  func2( self%map%get_x(c) )*jmatrix(2,component) + &
478  func3( self%map%get_x(c) )*jmatrix(3,component)) * &
479  (bspl_d1(k1,deg(1)+2-j1)*&
480  bspl_d2(k2,deg(2)+2-j2)*&
481  bspl_d3(k3,deg(3)+2-j3)) * abs(jacobian)
482  enddo
483  enddo
484  end do
485  end do
486  end do
487  end do
488  ! rescale by cell size
489  coefs_dofs(counter) = coef*self%volume
490  counter = counter+1
491  enddo
492  end do
493  end do
494 
495  end subroutine rhs_oneform
496 
497 
499  subroutine rhs_twoform( self, deg, func1, func2, func3, component, coefs_dofs )
500  class(sll_t_maxwell_3d_trafo) :: self
501  sll_int32, intent( in ) :: deg(3)
502  procedure(sll_i_function_3d_real64) :: func1
503  procedure(sll_i_function_3d_real64) :: func2
504  procedure(sll_i_function_3d_real64) :: func3
505  sll_int32, intent( in ) :: component
506  sll_real64, intent( out ) :: coefs_dofs(:)
507  ! local variables
508  sll_int32 :: i1, i2, i3,j1, j2, j3, k1, k2, k3, q(3), counter
509  sll_real64 :: c(3)
510  sll_real64 :: coef, jmatrix(3,3)
511  sll_real64, allocatable :: xw_gauss_d1(:,:), xw_gauss_d2(:,:), xw_gauss_d3(:,:)
512  sll_real64, allocatable :: bspl_d1(:,:), bspl_d2(:,:), bspl_d3(:,:)
513  sll_real64 :: scratch(maxval(deg+1))
514 
515  ! take enough Gauss points so that projection is exact for splines of degree deg
516  q = deg+1
517  ! rescale on [0,1] for compatibility with B-splines
518  allocate(xw_gauss_d1(2, q(1)))
519  allocate(bspl_d1(q(1), deg(1)+1 ))
520  xw_gauss_d1 = sll_f_gauss_legendre_points_and_weights(q(1), 0.0_f64, 1.0_f64)
521  ! Compute bsplines at gauss_points
522  do k1 = 1, q(1)
523  call sll_s_uniform_bsplines_eval_basis(deg(1),xw_gauss_d1(1,k1), scratch(1:deg(1)+1))
524  bspl_d1(k1,:) = scratch(1:deg(1)+1)
525  end do
526 
527  allocate(xw_gauss_d2(2, q(2)))
528  allocate(bspl_d2(q(2), deg(2)+1))
529  xw_gauss_d2 = sll_f_gauss_legendre_points_and_weights(q(2), 0.0_f64, 1.0_f64)
530  ! Compute bsplines at gauss_points
531  do k2 = 1, q(2)
532  call sll_s_uniform_bsplines_eval_basis(deg(2),xw_gauss_d2(1,k2), scratch(1:deg(2)+1))
533  bspl_d2(k2,:) = scratch(1:deg(2)+1)
534  end do
535 
536  allocate(xw_gauss_d3(2,q(3)))
537  allocate(bspl_d3(q(3), deg(3)+1 ))
538  xw_gauss_d3 = sll_f_gauss_legendre_points_and_weights(q(3), 0.0_f64, 1.0_f64)
539  ! Compute bsplines at gauss_points
540  do k3 = 1, q(3)
541  call sll_s_uniform_bsplines_eval_basis(deg(3),xw_gauss_d3(1,k3), scratch(1:deg(3)+1) )
542  bspl_d3(k3,:) = scratch(1:deg(3)+1)
543  end do
544 
545 
546  counter = 1
547  ! Compute coefs_dofs = int f(\xi)N_i(\xi) J_F(\xi) d\xi
548  do i3 = 1, self%n_dofs(3)
549  do i2 = 1, self%n_dofs(2)
550  do i1 = 1, self%n_dofs(1)
551  coef = 0.0_f64
552  ! loop over support of B spline
553  do j3 = 1, deg(3)+1
554  do j2 = 1, deg(2)+1
555  do j1 = 1, deg(1)+1
556  ! loop over Gauss points
557  do k3 = 1, q(3)
558  c(3) = self%delta_x(3) * (xw_gauss_d3(1,k3) + real(i3 + j3 - 2, f64))
559  do k2 = 1, q(2)
560  c(2) = self%delta_x(2) * (xw_gauss_d2(1,k2) + real(i2 + j2 - 2, f64))
561  do k1 = 1, q(1)
562  c(1) = self%delta_x(1) * (xw_gauss_d1(1,k1) + real(i1 + j1 - 2, f64))
563 
564  jmatrix = self%map%jacobian_matrix( c ) * sign( 1._f64, self%map%jacobian( c ) )
565  coef = coef + xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* xw_gauss_d3(2,k3) *&
566  (func1( self%map%get_x(c) ) * jmatrix(1,component) + &
567  func2( self%map%get_x(c) ) * jmatrix(2,component) + &
568  func3( self%map%get_x(c) ) * jmatrix(3,component)) *&
569  (bspl_d1(k1,deg(1)+2-j1)*&
570  bspl_d2(k2,deg(2)+2-j2)*&
571  bspl_d3(k3,deg(3)+2-j3))
572  enddo
573  enddo
574  end do
575  end do
576  end do
577  end do
578  ! rescale by cell size
579  coefs_dofs(counter) = coef*self%volume
580  counter = counter+1
581  enddo
582  end do
583  end do
584 
585  end subroutine rhs_twoform
586 
587 
589  subroutine rhs_threeform( self, deg, func, coefs_dofs )
590  class(sll_t_maxwell_3d_trafo) :: self
591  sll_int32, intent( in ) :: deg(3)
592  procedure(sll_i_function_3d_real64) :: func
593  sll_real64, intent( out ) :: coefs_dofs(:)
594 
595  ! local variables
596  sll_int32 :: i1, i2, i3,j1, j2, j3, k1, k2, k3, q(3), counter
597  sll_real64 :: coef, c(3)
598  sll_real64, allocatable :: xw_gauss_d1(:,:), xw_gauss_d2(:,:), xw_gauss_d3(:,:)
599  sll_real64, allocatable :: bspl_d1(:,:), bspl_d2(:,:), bspl_d3(:,:)
600 
601  ! take enough Gauss points so that projection is exact for splines of degree deg
602  q = deg+1
603  ! rescale on [0,1] for compatibility with B-splines
604  allocate(xw_gauss_d1(2, q(1)))
605  allocate(bspl_d1(q(1), deg(1)+1 ))
606  xw_gauss_d1 = sll_f_gauss_legendre_points_and_weights(q(1), 0.0_f64, 1.0_f64)
607  ! Compute bsplines at gauss_points
608  do k1 = 1, q(1)
609  call sll_s_uniform_bsplines_eval_basis(deg(1),xw_gauss_d1(1,k1), bspl_d1(k1,:))
610  end do
611 
612  allocate(xw_gauss_d2(2, q(2)))
613  allocate(bspl_d2(q(2), deg(2)+1))
614  xw_gauss_d2 = sll_f_gauss_legendre_points_and_weights(q(2), 0.0_f64, 1.0_f64)
615  ! Compute bsplines at gauss_points
616  do k2 = 1, q(2)
617  call sll_s_uniform_bsplines_eval_basis(deg(2),xw_gauss_d2(1,k2), bspl_d2(k2,:))
618  end do
619 
620  allocate(xw_gauss_d3(2,q(3)))
621  allocate(bspl_d3(q(3), deg(3)+1 ))
622  xw_gauss_d3 = sll_f_gauss_legendre_points_and_weights(q(3), 0.0_f64, 1.0_f64)
623  ! Compute bsplines at gauss_points
624  do k3 = 1, q(3)
625  call sll_s_uniform_bsplines_eval_basis(deg(3),xw_gauss_d3(1,k3), bspl_d3(k3,:))
626  end do
627 
628 
629  counter = 1
630  ! Compute coefs_dofs = int f(\xi)N_i(\xi) J_F(\xi) d\xi
631  do i3 = 1, self%n_dofs(3)
632  do i2 = 1, self%n_dofs(2)
633  do i1 = 1, self%n_dofs(1)
634  coef = 0.0_f64
635  ! loop over support of B spline
636  do j3 = 1, deg(3)+1
637  do j2 = 1, deg(2)+1
638  do j1 = 1, deg(1)+1
639  ! loop over Gauss points
640  do k3 = 1, q(3)
641  c(3) = self%delta_x(3) * (xw_gauss_d3(1,k3) + real(i3 + j3 - 2, f64))
642  do k2 = 1, q(2)
643  c(2) = self%delta_x(2) * (xw_gauss_d2(1,k2) + real(i2 + j2 - 2, f64))
644  do k1 = 1, q(1)
645  c(1) = self%delta_x(1) * (xw_gauss_d1(1,k1) + real(i1 + j1 - 2, f64))
646  coef = coef + xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* xw_gauss_d3(2,k3) *&
647  func(self%map%get_x(c) ) * &
648  bspl_d1(k1,deg(1)+2-j1)*&
649  bspl_d2(k2,deg(2)+2-j2)*&
650  bspl_d3(k3,deg(3)+2-j3)
651  enddo
652  enddo
653  end do
654  end do
655  end do
656  end do
657  ! rescale by cell size
658  coefs_dofs(counter) = coef*self%volume
659  counter = counter+1
660  enddo
661  end do
662  end do
663 
664  end subroutine rhs_threeform
665 
666 
668  subroutine l2projection_3d_trafo( self, form, component, coefs_dofs, func1, func2, func3 )
669  class(sll_t_maxwell_3d_trafo) :: self
670  sll_int32, intent( in ) :: form
671  sll_int32, intent( in ) :: component
672  sll_real64, intent( out ) :: coefs_dofs(:)
673  procedure(sll_i_function_3d_real64) :: func1
674  procedure(sll_i_function_3d_real64), optional :: func2
675  procedure(sll_i_function_3d_real64), optional :: func3
676 
677  if( present(func2) .and. present(func3) ) then
678  ! Compute right-hand-side
679  call sll_s_compute_rhs_trafo( self, form, 1, self%work(1:self%n_total), func1, func2, func3 )
680  call sll_s_compute_rhs_trafo( self, form, 2, self%work(1+self%n_total:2*self%n_total), func1, func2, func3 )
681  call sll_s_compute_rhs_trafo( self, form, 3, self%work(1+2*self%n_total:3*self%n_total), func1, func2, func3 )
682 
683  select case( form )
684  case( 1 )
685  call self%mass1_solver%solve( self%work, coefs_dofs )
686  case( 2 )
687  call self%mass2_solver%solve( self%work, coefs_dofs )
688  case default
689  print*, 'L2projection for', form, '-form not implemented.'
690  end select
691  else
692  select case( form )
693  case(0)
694  call sll_s_compute_rhs_trafo( self, form, 0, self%work0, func1 )
695  call self%mass0_solver%solve(self%work0, coefs_dofs)
696  case default
697  print*,'l2 projection not for single function or component implemented'
698  end select
699  end if
700 
701 
702  end subroutine l2projection_3d_trafo
703 
704 
706  function l2norm_squared_3d_trafo( self, coefs, form, component ) result ( r )
707  class(sll_t_maxwell_3d_trafo) :: self
708  sll_real64 :: coefs(:)
709  sll_int32 :: form
710  sll_int32 :: component
711  sll_real64 :: r
712 
713  r = inner_product_3d_trafo( self, coefs, coefs, form, component )
714 
715  end function l2norm_squared_3d_trafo
716 
717 
719  function inner_product_3d_trafo( self, coefs1, coefs2, form, component ) result ( r )
720  class(sll_t_maxwell_3d_trafo) :: self
721  sll_real64 :: coefs1(:)
722  sll_real64 :: coefs2(:)
723  sll_int32 :: form
724  sll_int32, optional :: component
725  sll_real64 :: r
726 
727  !local variables
728  sll_int32 :: istart, iend
729 
730  if ( form == 0 ) then
731  call self%multiply_mass( [form], coefs2, self%work0)
732  r = sum( coefs1*self%work0)
733  else
734  if (present(component)) then
735  istart = 1+(component-1)*self%n_total
736  iend = component*self%n_total
737  else
738  istart=1
739  iend=3*self%n_total
740  end if
741  call multiply_mass_trafo( self, [form], coefs2, self%work(1:3*self%n_total) )
742  r = sum(coefs1(istart:iend)*self%work(istart:iend))
743  end if
744 
745  end function inner_product_3d_trafo
746 
747 
749  subroutine init_3d_trafo( self, domain, n_dofs, s_deg_0, map, mass_tolerance, poisson_tolerance, solver_tolerance, adiabatic_electrons, profile )
750  class(sll_t_maxwell_3d_trafo), intent( inout ) :: self
751  sll_real64, intent(in) :: domain(3,2)
752  sll_int32, intent( in ) :: n_dofs(3)
753  sll_int32, intent( in ) :: s_deg_0(3)
754  type(sll_t_mapping_3d), target, intent( inout ) :: map
755  sll_real64, intent(in), optional :: mass_tolerance
756  sll_real64, intent(in), optional :: poisson_tolerance
757  sll_real64, intent(in), optional :: solver_tolerance
758  logical, intent(in), optional :: adiabatic_electrons
759  type(sll_t_profile_functions), intent(in), optional :: profile
760  ! local variables
761  sll_int32 :: j
762  sll_int32 :: deg1(3), deg2(3)
763  sll_real64, allocatable :: nullspace(:,:)!, ones(:)
764 
765  self%Lx = map%Lx
766 
767  if (present( mass_tolerance) ) then
768  self%mass_solver_tolerance = mass_tolerance
769  else
770  self%mass_solver_tolerance = 1d-12
771  end if
772  if (present( poisson_tolerance) ) then
773  self%poisson_solver_tolerance = poisson_tolerance
774  else
775  self%poisson_solver_tolerance = 1d-12
776  end if
777  if (present( solver_tolerance) ) then
778  self%solver_tolerance = solver_tolerance
779  else
780  self%solver_tolerance = 1d-12
781  end if
782 
783  if( present( adiabatic_electrons ) ) then
784  self%adiabatic_electrons = adiabatic_electrons
785  end if
786 
787  if( present( profile ) ) then
788  self%profile = profile
789  end if
790 
791  self%map=>map
792  self%n_cells = n_dofs
793  self%n_dofs = n_dofs
794  self%n_total = product(n_dofs)
795  self%n_total0 = self%n_total
796  self%n_total1 = self%n_total
797 
798  self%delta_x = 1._f64 / real(n_dofs,f64)
799  self%s_deg_0 = s_deg_0
800  self%s_deg_1 = s_deg_0 - 1
801  self%volume = product(self%delta_x)
802 
803  ! Allocate scratch data
804  allocate( self%work0(1:self%n_total) )
805  allocate( self%work(1:self%n_total*3) )
806  allocate( self%work2(1:self%n_total*3) )
807 
808  ! Assemble the sparse mass matrices
809  !0-form
810  deg1=self%s_deg_0
811  if(self%adiabatic_electrons) then
812 
813  call sll_s_spline_fem_mass3d( n_dofs, deg1, -1, self%mass0, profile_m0 )
814  else
815 
816  call sll_s_spline_fem_mass3d( n_dofs, deg1, 0, self%mass0, profile_0 )
817  end if
818  do j=1, 3
819  deg1=self%s_deg_0
820  deg1(j)=self%s_deg_1(j)
821  call sll_s_spline_fem_mass3d( n_dofs, deg1, j, self%mass1(j,j), profile_1 )
822 
823  deg1 = self%s_deg_1
824  deg1(j) = self%s_deg_0(j)
825  call sll_s_spline_fem_mass3d( n_dofs, deg1, j, self%mass2(j,j), profile_2 )
826 
827  end do
828  if(self%map%flag2d)then
829  deg1=self%s_deg_0
830  deg1(1)=self%s_deg_1(1)
831  deg2=self%s_deg_0
832  deg2(2) = self%s_deg_1(2)
833  call sll_s_spline_fem_mixedmass3d( n_dofs, deg1, deg2, [1,2], self%mass1(1,2), profile_m1 )
834  call sll_s_spline_fem_mixedmass3d( n_dofs, deg2, deg1, [2,1], self%mass1(2,1), profile_m1 )
835 
836  deg1 = self%s_deg_1
837  deg1(1) = self%s_deg_0(1)
838  deg2 = self%s_deg_1
839  deg2(2) = self%s_deg_0(2)
840  call sll_s_spline_fem_mixedmass3d( n_dofs, deg1, deg2, [1,2], self%mass2(1,2), profile_m2 )
841  call sll_s_spline_fem_mixedmass3d( n_dofs, deg2, deg1, [2,1], self%mass2(2,1), profile_m2 )
842 
843  if(self%map%flag3d)then
844  do j = 2, 3
845  deg1=self%s_deg_0
846  deg1(j)=self%s_deg_1(j)
847  deg2=self%s_deg_0
848  deg2(modulo(j,3)+1) = self%s_deg_1(modulo(j,3)+1)
849  call sll_s_spline_fem_mixedmass3d( n_dofs, deg1, deg2, [j,modulo(j,3)+1], self%mass1(j,modulo(j,3)+1), profile_m1 )
850  call sll_s_spline_fem_mixedmass3d( n_dofs, deg2, deg1, [modulo(j,3)+1,j], self%mass1(modulo(j,3)+1,j), profile_m1 )
851 
852 
853  deg1 = self%s_deg_1
854  deg1(j) = self%s_deg_0(j)
855  deg2 = self%s_deg_1
856  deg2(modulo(j,3)+1) = self%s_deg_0(modulo(j,3)+1)
857  call sll_s_spline_fem_mixedmass3d( n_dofs, deg1, deg2, [j,modulo(j,3)+1], self%mass2(j,modulo(j,3)+1), profile_m2 )
858  call sll_s_spline_fem_mixedmass3d( n_dofs, deg2, deg1, [modulo(j,3)+1,j], self%mass2(modulo(j,3)+1,j), profile_m2 )
859  end do
860  end if
861  end if
862 
863  call self%mass1_operator%create( 3, 3 )
864  call self%mass2_operator%create( 3, 3 )
865 
866  do j=1,3
867  call self%mass1_operator%set( j, j, self%mass1(j,j) )
868  call self%mass2_operator%set( j, j, self%mass2(j,j) )
869  end do
870  if(self%map%flag2d)then
871  call self%mass1_operator%set( 1, 2, self%mass1(1, 2) )
872  call self%mass1_operator%set( 2, 1, self%mass1(2, 1) )
873  call self%mass2_operator%set( 1, 2, self%mass2(1, 2) )
874  call self%mass2_operator%set( 2, 1, self%mass2(2, 1) )
875  if(self%map%flag3d)then
876  do j = 2, 3
877  call self%mass1_operator%set( j, modulo(j,3)+1, self%mass1(j, modulo(j,3)+1) )
878  call self%mass1_operator%set( modulo(j,3)+1, j, self%mass1(modulo(j,3)+1, j) )
879  call self%mass2_operator%set( j, modulo(j,3)+1, self%mass2(j, modulo(j,3)+1) )
880  call self%mass2_operator%set( modulo(j,3)+1, j, self%mass2(modulo(j,3)+1, j) )
881  end do
882  end if
883  end if
884 
885  call self%preconditioner_fft%init( self%Lx, n_dofs, s_deg_0 )
886  call self%mass0_solver%create( self%mass0 )
887  call self%mass1_solver%create( self%mass1_operator, self%preconditioner_fft%inverse_mass1_3d)
888  call self%mass2_solver%create( self%mass2_operator, self%preconditioner_fft%inverse_mass2_3d)
889  self%mass1_solver%atol = self%mass_solver_tolerance/maxval(self%Lx)
890  self%mass2_solver%atol = self%mass_solver_tolerance/maxval(self%Lx)**2
891  !self%mass1_solver%verbose = .true.
892  !self%mass2_solver%verbose = .true.
893 
894  call self%poisson_matrix%create( self%mass1_operator, self%n_dofs, self%delta_x )
895  ! Penalized Poisson operator
896  allocate(nullspace(1,self%n_total))
897  nullspace(1,:) = 1.0_f64
898  call self%poisson_operator%create( self%poisson_matrix, nullspace(:,1:self%n_total), 1 )
899  ! Poisson solver
900  call self%poisson_solver%create( self%poisson_operator )
901  self%poisson_solver%null_space = .true.
902  self%poisson_solver%atol = self%poisson_solver_tolerance
903  !self%poisson_solver%verbose = .true.
904  self%poisson_solver%n_maxiter=40000
905 
906  ! Only for Schur complement eb solver
907  call self%linear_op_schur_eb%create( self%mass1_operator, self%mass2_operator, self%n_total, self%n_dofs, self%delta_x )
908  call self%linear_solver_schur_eb%create( self%linear_op_schur_eb, self%preconditioner_fft%inverse_mass1_3d)
909  self%linear_solver_schur_eb%atol = self%solver_tolerance/maxval(self%Lx)
910  !self%linear_solver_schur_eb%verbose = .true.
911  !self%linear_solver_schur_eb%n_maxiter = 2000
912 
913  contains
914  function profile_m0( x, component)
915  sll_real64 :: profile_m0
916  sll_real64, intent(in) :: x(3)
917  sll_int32, optional, intent(in) :: component(:)
918 
919  profile_m0 = self%map%jacobian( x ) * self%profile%rho_0( x(1) )/ self%profile%T_e( x(1) )
920 
921  end function profile_m0
922 
923  function profile_0(x, component)
924  sll_real64 :: profile_0
925  sll_real64, intent(in) :: x(3)
926  sll_int32, optional, intent(in) :: component(:)
927 
928  profile_0 = self%map%jacobian( x )
929 
930  end function profile_0
931 
932  function profile_1(x, component)
933  sll_real64 :: profile_1
934  sll_real64, intent(in) :: x(3)
935  sll_int32, optional, intent(in) :: component(:)
936 
937  profile_1 = self%map%metric_inverse_single_jacobian( x, component(1), component(1) )
938 
939  end function profile_1
940 
941  function profile_2(x, component)
942  sll_real64 :: profile_2
943  sll_real64, intent(in) :: x(3)
944  sll_int32, optional, intent(in) :: component(:)
945 
946  profile_2 = self%map%metric_single_jacobian( x, component(1), component(1) )
947 
948  end function profile_2
949 
950  function profile_m1(x, component)
951  sll_real64 :: profile_m1
952  sll_real64, intent(in) :: x(3)
953  sll_int32, optional, intent(in) :: component(:)
954 
955  profile_m1 = self%map%metric_inverse_single_jacobian( x, component(1), component(2) )
956  end function profile_m1
957 
958  function profile_m2(x, component)
959  sll_real64 :: profile_m2
960  sll_real64, intent(in) :: x(3)
961  sll_int32, optional, intent(in) :: component(:)
962 
963  profile_m2 = self%map%metric_single_jacobian( x, component(1), component(2) )
964 
965  end function profile_m2
966 
967  end subroutine init_3d_trafo
968 
969 
971  subroutine init_from_file_3d_trafo( self, domain, n_dofs, s_deg_0, map, nml_file, adiabatic_electrons, profile )
972  class(sll_t_maxwell_3d_trafo), intent( inout ) :: self
973  sll_real64, intent(in) :: domain(3,2)
974  sll_int32, intent( in ) :: n_dofs(3)
975  sll_int32, intent( in ) :: s_deg_0(3)
976  type(sll_t_mapping_3d), target, intent( inout ) :: map
977  character(len=*), intent(in) :: nml_file
978  logical, intent(in), optional :: adiabatic_electrons
979  type(sll_t_profile_functions), intent(in), optional :: profile
980  ! local variables
981  character(len=256) :: file_prefix
982  sll_int32 :: input_file
983  sll_int32 :: io_stat, io_stat0, io_stat1, rank, file_id
984  sll_real64 :: mass_tolerance
985  sll_real64 :: poisson_tolerance
986  sll_real64 :: maxwell_tolerance
987 
988  namelist /output/ file_prefix
989  namelist /maxwell_solver/ mass_tolerance, poisson_tolerance
990  namelist /time_solver/ maxwell_tolerance
991 
992  rank = sll_f_get_collective_rank(sll_v_world_collective)
993 
994  if( present( adiabatic_electrons ) ) then
995  self%adiabatic_electrons = adiabatic_electrons
996  end if
997 
998  if( present( profile ) ) then
999  self%profile = profile
1000  end if
1001 
1002  ! Read in solver tolerance
1003  open(newunit = input_file, file=nml_file, status='old',iostat=io_stat)
1004  if (io_stat /= 0) then
1005  if (rank == 0 ) then
1006  print*, 'sll_m_maxwell_3d_trafo: Input file does not exist. Set default tolerance.'
1007  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1008  write(file_id, *) 'mass solver tolerance:', 1d-12
1009  write(file_id, *) 'poisson solver tolerance:', 1d-12
1010  close(file_id)
1011  end if
1012  call self%init( domain, n_dofs, s_deg_0, map )
1013  else
1014  read(input_file, output,iostat=io_stat0)
1015  read(input_file, maxwell_solver,iostat=io_stat)
1016  read(input_file, time_solver,iostat=io_stat1)
1017  if (io_stat /= 0 .and. io_stat1 /= 0) then
1018  if (rank == 0 ) then
1019  print*, 'sll_m_maxwell_3d_trafo: Input parameter does not exist. Set default tolerance.'
1020  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1021  write(file_id, *) 'mass solver tolerance:', 1d-12
1022  write(file_id, *) 'poisson solver tolerance:', 1d-12
1023  close(file_id)
1024  end if
1025  call self%init( domain, n_dofs, s_deg_0, map )
1026  else if (io_stat /= 0 .and. io_stat1 == 0) then
1027  call self%init( domain, n_dofs, s_deg_0, map, solver_tolerance=maxwell_tolerance)
1028  else if (io_stat == 0 .and. io_stat1 /= 0) then
1029  if (rank == 0 ) then
1030  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1031  write(file_id, *) 'mass solver tolerance:', mass_tolerance
1032  write(file_id, *) 'poisson solver tolerance:', poisson_tolerance
1033  close(file_id)
1034  end if
1035  call self%init( domain, n_dofs, s_deg_0, map, mass_tolerance, poisson_tolerance )
1036  else
1037  if (rank == 0 ) then
1038  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1039  write(file_id, *) 'mass solver tolerance:', mass_tolerance
1040  write(file_id, *) 'poisson solver tolerance:', poisson_tolerance
1041  close(file_id)
1042  end if
1043  call self%init( domain, n_dofs, s_deg_0, map, mass_tolerance, poisson_tolerance, maxwell_tolerance )
1044  end if
1045  close(input_file)
1046  end if
1047 
1048  end subroutine init_from_file_3d_trafo
1049 
1050 
1052  subroutine free_3d_trafo( self )
1053  class(sll_t_maxwell_3d_trafo) :: self
1054  self%map => null()
1055  call self%preconditioner_fft%free()
1056  call self%mass0_solver%free()
1057  call self%mass1_solver%free()
1058  call self%mass2_solver%free()
1059  call self%mass1_operator%free()
1060  call self%mass2_operator%free()
1061  call self%poisson_matrix%free()
1062  call self%poisson_solver%free()
1063  call self%poisson_operator%free()
1064  call self%linear_solver_schur_eb%free()
1065  call self%linear_op_schur_eb%free()
1066  deallocate(self%work0)
1067  deallocate(self%work)
1068  deallocate(self%work2)
1069 
1070  end subroutine free_3d_trafo
1071 
1072 
1074  subroutine multiply_g( self, field_in, field_out )
1075  class(sll_t_maxwell_3d_trafo) :: self
1076  sll_real64, intent( in ) :: field_in(:)
1077  sll_real64, intent( out ) :: field_out(:)
1078 
1079  call sll_s_multiply_g(self%n_dofs, self%delta_x, field_in, field_out)
1080 
1081  end subroutine multiply_g
1082 
1083 
1085  subroutine multiply_gt(self, field_in, field_out)
1086  class(sll_t_maxwell_3d_trafo) :: self
1087  sll_real64, intent( in ) :: field_in(:)
1088  sll_real64, intent( out ) :: field_out(:)
1089 
1090  call sll_s_multiply_gt(self%n_dofs, self%delta_x, field_in, field_out)
1091 
1092  end subroutine multiply_gt
1093 
1094 
1096  subroutine multiply_c(self, field_in, field_out)
1097  class(sll_t_maxwell_3d_trafo) :: self
1098  sll_real64, intent( in ) :: field_in(:)
1099  sll_real64, intent( out ) :: field_out(:)
1100 
1101  call sll_s_multiply_c(self%n_dofs, self%delta_x, field_in, field_out)
1102 
1103  end subroutine multiply_c
1104 
1105 
1107  subroutine multiply_ct(self, field_in, field_out)
1108  class(sll_t_maxwell_3d_trafo) :: self
1109  sll_real64, intent( in ) :: field_in(:)
1110  sll_real64, intent( out ) :: field_out(:)
1111 
1112  call sll_s_multiply_ct(self%n_dofs, self%delta_x, field_in, field_out)
1113 
1114  end subroutine multiply_ct
1115 
1116 
1118  subroutine multiply_mass_trafo( self, deg, coefs_in, coefs_out )
1119  class(sll_t_maxwell_3d_trafo) :: self
1120  sll_int32, intent( in ) :: deg(:)
1121  sll_real64, intent( in ) :: coefs_in(:)
1122  sll_real64, intent( out ) :: coefs_out(:)
1123 
1124  if( size(deg) == 1 )then
1125  sll_assert(deg(1)>=0 .and. deg(1)<=2)
1126  select case(deg(1))
1127  case(0)
1128  call self%mass0%dot( coefs_in, coefs_out )
1129  case(1)
1130  call self%mass1_operator%dot( coefs_in, coefs_out )
1131  case(2)
1132  call self%mass2_operator%dot( coefs_in, coefs_out )
1133  case default
1134  print*, 'multiply mass for other form not yet implemented'
1135  stop
1136  end select
1137  else if( size(deg) == 3 ) then
1138  coefs_out = 0._f64
1139  end if
1140 
1141  end subroutine multiply_mass_trafo
1142 
1143 
1145  subroutine multiply_mass_inverse_trafo( self, form, coefs_in, coefs_out )
1146  class(sll_t_maxwell_3d_trafo) :: self
1147  sll_int32, intent( in ) :: form
1148  sll_real64, intent( in ) :: coefs_in(:)
1149  sll_real64, intent( out ) :: coefs_out(:)
1150 
1151 
1152  sll_assert(form>=0 .and. form<=2)
1153  select case(form)
1154  case(0)
1155  call self%mass0_solver%solve( coefs_in, coefs_out )
1156  case(1)
1157  call self%mass1_solver%solve( coefs_in, coefs_out )
1158  case(2)
1159  call self%mass2_solver%solve( coefs_in, coefs_out )
1160  case default
1161  print*, 'multiply inverse mass for other form not yet implemented'
1162  stop
1163  end select
1164  end subroutine multiply_mass_inverse_trafo
1165 
1166 
1167 end module sll_m_maxwell_3d_trafo
sll_m_maxwell_3d_base::inner_product
Definition: sll_m_maxwell_3d_base.F90:192
sll_m_maxwell_3d_base::multiply_mass_inverse
Definition: sll_m_maxwell_3d_base.F90:232
sll_m_maxwell_3d_base::multiply_mass
Definition: sll_m_maxwell_3d_base.F90:220
sll_m_maxwell_3d_base::sll_i_function_3d_real64
3d real function
Definition: sll_m_maxwell_3d_base.F90:147
sll_m_collective
Parallelizing facility.
Definition: sll_m_collective.F90:128
sll_m_collective::sll_f_get_collective_rank
integer(kind=i32) function, public sll_f_get_collective_rank(col)
Determines the rank of the calling process in the communicator.
Definition: sll_m_collective.F90:533
sll_m_collective::sll_v_world_collective
type(sll_t_collective_t), pointer, public sll_v_world_collective
Control of subset assignment. Processes with the same color are in the same new communicator.
Definition: sll_m_collective.F90:255
sll_m_gauss_legendre_integration
Gauss-Legendre integration.
Definition: sll_m_gauss_legendre_integration.F90:8
sll_m_gauss_legendre_integration::sll_f_gauss_legendre_points_and_weights
real(kind=f64) function, dimension(2, 1:npoints), public sll_f_gauss_legendre_points_and_weights(npoints, a, b)
Returns a 2d array of size (2,npoints) containing gauss-legendre points and weights in the interval [...
Definition: sll_m_gauss_legendre_integration.F90:265
sll_m_linear_operator_block
module for a block linear operator
Definition: sll_m_linear_operator_block.F90:11
sll_m_linear_operator_penalized
module for a penalized linear operator
Definition: sll_m_linear_operator_penalized.F90:11
sll_m_linear_operator_poisson_3d
Definition: sll_m_linear_operator_poisson_3d.F90:1
sll_m_linear_operator_schur_eb_3d
Definition: sll_m_linear_operator_schur_eb_3d.F90:1
sll_m_linear_solver_cg
module for conjugate gradient method in pure form
Definition: sll_m_linear_solver_cg.F90:11
sll_m_low_level_bsplines
Low level arbitrary degree splines.
Definition: sll_m_low_level_bsplines.F90:27
sll_m_mapping_3d
Module interfaces for coordinate transformation.
Definition: sll_m_mapping_3d.F90:13
sll_m_matrix_csr
module for Compressed Sparse Row Matrix (CSR)
Definition: sll_m_matrix_csr.F90:12
sll_m_maxwell_3d_base
Module interface to solve Maxwell's equations in 3D.
Definition: sll_m_maxwell_3d_base.F90:7
sll_m_maxwell_3d_trafo
Module interface to solve Maxwell's equations with coordinate transformation in 3D.
Definition: sll_m_maxwell_3d_trafo.F90:9
sll_m_maxwell_3d_trafo::rhs_zeroform
subroutine rhs_zeroform(self, deg, func, coefs_dofs)
Compute $\int f(F(\xi)) N_i(\xi) J_F(\xi) d\xi$ for scalar function f, where $N_i$ is the B-spline st...
Definition: sll_m_maxwell_3d_trafo.F90:327
sll_m_maxwell_3d_trafo::multiply_c
subroutine multiply_c(self, field_in, field_out)
Multiply by discrete curl matrix.
Definition: sll_m_maxwell_3d_trafo.F90:1097
sll_m_maxwell_3d_trafo::sll_s_compute_rhs_trafo
subroutine sll_s_compute_rhs_trafo(self, form, component, coefs_dofs, func1, func2, func3)
Compute the FEM right-hand-side for a given function f and periodic splines of given degree Its compo...
Definition: sll_m_maxwell_3d_trafo.F90:295
sll_m_maxwell_3d_trafo::rhs_twoform
subroutine rhs_twoform(self, deg, func1, func2, func3, component, coefs_dofs)
Compute $\int f(F(\xi))\cdot DF(\xi) N_i(\xi) d\xi$ for vector function f, where $N_i$ is the B-splin...
Definition: sll_m_maxwell_3d_trafo.F90:500
sll_m_maxwell_3d_trafo::sll_s_compute_phi_from_rho_3d_trafo
subroutine sll_s_compute_phi_from_rho_3d_trafo(self, field_in, field_out, efield_dofs)
Compute phi from rho_i integrated over the time interval.
Definition: sll_m_maxwell_3d_trafo.F90:263
sll_m_maxwell_3d_trafo::rhs_oneform
subroutine rhs_oneform(self, deg, func1, func2, func3, component, coefs_dofs)
Compute $\int f(F(\xi))\cdot DF^{-T}(\xi) N_i(\xi) J_F(\xi) d\xi$ for vector function f,...
Definition: sll_m_maxwell_3d_trafo.F90:410
sll_m_maxwell_3d_trafo::free_3d_trafo
subroutine free_3d_trafo(self)
Finalization.
Definition: sll_m_maxwell_3d_trafo.F90:1053
sll_m_maxwell_3d_trafo::inner_product_3d_trafo
real(kind=f64) function inner_product_3d_trafo(self, coefs1, coefs2, form, component)
Compute inner product.
Definition: sll_m_maxwell_3d_trafo.F90:720
sll_m_maxwell_3d_trafo::rhs_threeform
subroutine rhs_threeform(self, deg, func, coefs_dofs)
Compute $\int f(F(\xi)) N_i(\xi) d\xi$ for scalar function f, where $N_i$ is the B-spline starting at...
Definition: sll_m_maxwell_3d_trafo.F90:590
sll_m_maxwell_3d_trafo::sll_s_compute_rho_from_e_3d_trafo
subroutine sll_s_compute_rho_from_e_3d_trafo(self, field_in, field_out)
compute rho from e using weak Gauss law ( rho = G^T M_1 e )
Definition: sll_m_maxwell_3d_trafo.F90:235
sll_m_maxwell_3d_trafo::multiply_g
subroutine multiply_g(self, field_in, field_out)
Multiply by dicrete gradient matrix.
Definition: sll_m_maxwell_3d_trafo.F90:1075
sll_m_maxwell_3d_trafo::sll_s_compute_curl_part_3d_trafo
subroutine sll_s_compute_curl_part_3d_trafo(self, delta_t, efield, bfield, betar)
Solve curl part of Maxwell's equations.
Definition: sll_m_maxwell_3d_trafo.F90:183
sll_m_maxwell_3d_trafo::multiply_gt
subroutine multiply_gt(self, field_in, field_out)
Multiply by transpose of dicrete gradient matrix.
Definition: sll_m_maxwell_3d_trafo.F90:1086
sll_m_maxwell_3d_trafo::l2norm_squared_3d_trafo
real(kind=f64) function l2norm_squared_3d_trafo(self, coefs, form, component)
Compute square of the L2norm.
Definition: sll_m_maxwell_3d_trafo.F90:707
sll_m_maxwell_3d_trafo::sll_s_compute_e_from_b_3d_trafo
subroutine sll_s_compute_e_from_b_3d_trafo(self, delta_t, field_in, field_out)
compute E from B using weak Ampere formulation
Definition: sll_m_maxwell_3d_trafo.F90:151
sll_m_maxwell_3d_trafo::sll_s_compute_e_from_j_3d_trafo
subroutine sll_s_compute_e_from_j_3d_trafo(self, current, E, component)
Compute E_i from j_i integrated over the time interval using weak Ampere formulation,...
Definition: sll_m_maxwell_3d_trafo.F90:250
sll_m_maxwell_3d_trafo::multiply_mass_inverse_trafo
subroutine multiply_mass_inverse_trafo(self, form, coefs_in, coefs_out)
Multiply by the inverse mass matrix.
Definition: sll_m_maxwell_3d_trafo.F90:1146
sll_m_maxwell_3d_trafo::l2projection_3d_trafo
subroutine l2projection_3d_trafo(self, form, component, coefs_dofs, func1, func2, func3)
Compute the L2 projection of a given function f on periodic splines of given degree.
Definition: sll_m_maxwell_3d_trafo.F90:669
sll_m_maxwell_3d_trafo::init_3d_trafo
subroutine init_3d_trafo(self, domain, n_dofs, s_deg_0, map, mass_tolerance, poisson_tolerance, solver_tolerance, adiabatic_electrons, profile)
Initialization.
Definition: sll_m_maxwell_3d_trafo.F90:750
sll_m_maxwell_3d_trafo::sll_s_compute_b_from_e_3d_trafo
subroutine sll_s_compute_b_from_e_3d_trafo(self, delta_t, field_in, field_out)
Compute B from E using strong 3D Faraday equation for spline coefficients.
Definition: sll_m_maxwell_3d_trafo.F90:169
sll_m_maxwell_3d_trafo::sll_s_compute_e_from_rho_3d_trafo
subroutine sll_s_compute_e_from_rho_3d_trafo(self, field_in, field_out)
Compute E_i from rho_i integrated over the time interval using weak Poisson's equation.
Definition: sll_m_maxwell_3d_trafo.F90:222
sll_m_maxwell_3d_trafo::init_from_file_3d_trafo
subroutine init_from_file_3d_trafo(self, domain, n_dofs, s_deg_0, map, nml_file, adiabatic_electrons, profile)
Initialization from nml file.
Definition: sll_m_maxwell_3d_trafo.F90:972
sll_m_maxwell_3d_trafo::sll_s_compute_phi_from_j_3d_trafo
subroutine sll_s_compute_phi_from_j_3d_trafo(self, field_in, field_out, efield_dofs)
Compute phi from j_i integrated over the time interval, delta_t is already included.
Definition: sll_m_maxwell_3d_trafo.F90:277
sll_m_maxwell_3d_trafo::multiply_mass_trafo
subroutine multiply_mass_trafo(self, deg, coefs_in, coefs_out)
Multiply by the mass matrix.
Definition: sll_m_maxwell_3d_trafo.F90:1119
sll_m_maxwell_3d_trafo::multiply_ct
subroutine multiply_ct(self, field_in, field_out)
Multiply by transpose of discrete curl matrix.
Definition: sll_m_maxwell_3d_trafo.F90:1108
sll_m_preconditioner_fft
Module interface to solve Maxwell's equations with coordinate transformation in 3D The linear systems...
Definition: sll_m_preconditioner_fft.F90:10
sll_m_profile_functions
functions for initial profile of the particle distribution function
Definition: sll_m_profile_functions.F90:4
sll_m_spline_fem_utilities_3d
Utilites for 3D Maxwell solvers with spline finite elements.
Definition: sll_m_spline_fem_utilities_3d.F90:9
sll_m_spline_fem_utilities_3d::sll_s_spline_fem_mass3d
subroutine, public sll_s_spline_fem_mass3d(n_cells, deg, component, matrix, profile, n_cells_min, n_cells_max)
Set up 3d mass matrix for specific spline degrees and profile function.
Definition: sll_m_spline_fem_utilities_3d.F90:56
sll_m_spline_fem_utilities_3d::sll_s_spline_fem_mixedmass3d
subroutine, public sll_s_spline_fem_mixedmass3d(n_cells, deg1, deg2, component, matrix, profile, n_cells_min, n_cells_max)
Set up 3d mixed mass matrix for specific spline degree and profile function.
Definition: sll_m_spline_fem_utilities_3d.F90:240
sll_m_spline_fem_utilities
Utilites for Maxwell solver's with spline finite elements.
Definition: sll_m_spline_fem_utilities.F90:9
sll_m_spline_fem_utilities::sll_s_multiply_g
subroutine, public sll_s_multiply_g(n_dofs, delta_x, field_in, field_out)
Multiply by dicrete gradient matrix (3d version)
Definition: sll_m_spline_fem_utilities.F90:772
sll_m_spline_fem_utilities::sll_s_multiply_c
subroutine, public sll_s_multiply_c(n_dofs, delta_x, field_in, field_out)
Multiply by discrete curl matrix.
Definition: sll_m_spline_fem_utilities.F90:928
sll_m_spline_fem_utilities::sll_s_multiply_ct
subroutine, public sll_s_multiply_ct(n_dofs, delta_x, field_in, field_out)
Multiply by transpose of discrete curl matrix.
Definition: sll_m_spline_fem_utilities.F90:1058
sll_m_spline_fem_utilities::sll_s_multiply_gt
subroutine, public sll_s_multiply_gt(n_dofs, delta_x, field_in, field_out)
Multiply by transpose of dicrete gradient matrix (3d version)
Definition: sll_m_spline_fem_utilities.F90:846
profile_0
real(kind=f64) function profile_0(x)
Definition: sll_m_maxwell_1d_trafo.F90:525
profile_m0
real(kind=f64) function profile_m0(x)
Definition: sll_m_maxwell_1d_trafo.F90:541
profile_1
real(kind=f64) function profile_1(x)
Definition: sll_m_maxwell_1d_trafo.F90:533
profile_m2
real(kind=f64) function profile_m2(x, component)
Definition: sll_m_maxwell_3d_trafo.F90:959
profile_2
real(kind=f64) function profile_2(x, component)
Definition: sll_m_maxwell_3d_trafo.F90:942
profile_m1
real(kind=f64) function profile_m1(x, component)
Definition: sll_m_maxwell_3d_trafo.F90:951
sll_m_linear_operator_block::sll_t_linear_operator_block
class for a linear operator_block
Definition: sll_m_linear_operator_block.F90:34
sll_m_linear_operator_penalized::sll_t_linear_operator_penalized
class for a linear operator_penalized
Definition: sll_m_linear_operator_penalized.F90:27
sll_m_linear_operator_poisson_3d::sll_t_linear_operator_poisson_3d
Definition: sll_m_linear_operator_poisson_3d.F90:20
sll_m_linear_operator_schur_eb_3d::sll_t_linear_operator_schur_eb_3d
Definition: sll_m_linear_operator_schur_eb_3d.F90:18
sll_m_linear_solver_cg::sll_t_linear_solver_cg
class for the cg linear solver
Definition: sll_m_linear_solver_cg.F90:38
sll_m_mapping_3d::sll_t_mapping_3d
type collecting functions for analytical coordinate mapping
Definition: sll_m_mapping_3d.F90:68
sll_m_matrix_csr::sll_t_matrix_csr
class for a csr matrix
Definition: sll_m_matrix_csr.F90:33
sll_m_maxwell_3d_base::sll_c_maxwell_3d_base
Definition: sll_m_maxwell_3d_base.F90:27
sll_m_maxwell_3d_trafo::sll_t_maxwell_3d_trafo
Definition: sll_m_maxwell_3d_trafo.F90:74
sll_m_preconditioner_fft::sll_t_preconditioner_fft
Definition: sll_m_preconditioner_fft.F90:35
sll_m_profile_functions::sll_t_profile_functions
Definition: sll_m_profile_functions.F90:19
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1