Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_maxwell_3d_fem_fft.F90
Go to the documentation of this file.
1 
9 ! TODO: Write FFT-based mass solver: There is such a solver already defined as linear_solver_mass1 in particle_methods. Reuse? Can we put the parallelization in this solver?
10 
11 
12 module sll_m_maxwell_3d_fem_fft
13  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
14 #include "sll_assert.h"
15 #include "sll_errors.h"
16 #include "sll_memory.h"
17 #include "sll_working_precision.h"
18 
19 
20  use sll_m_gauss_legendre_integration, only: &
21  sll_f_gauss_legendre_points_and_weights
22 
23  use sll_m_linear_operator_kron, only : &
24  sll_t_linear_operator_kron
25 
26  use sll_m_linear_operator_maxwell_eb_schur, only : &
27  sll_t_linear_operator_maxwell_eb_schur
28 
29  use sll_m_linear_solver_spline_mass_fft, only : &
30  sll_t_linear_solver_spline_mass_fft
31 
32  use sll_m_low_level_bsplines, only: &
33  sll_s_uniform_bsplines_eval_basis
34 
35  use sll_m_matrix_csr, only: &
36  sll_t_matrix_csr
37 
38  use sll_m_maxwell_3d_base, only: &
39  sll_i_function_3d_real64, &
40  sll_c_maxwell_3d_base
41 
42  use sll_m_poisson_3d_fem_fft, only : &
43  sll_t_poisson_3d_fem_fft
44 
45  use sll_m_profile_functions, only: &
46  sll_t_profile_functions
47 
48  use sll_m_spline_fem_utilities, only : &
49  sll_s_spline_fem_mass_line, &
50  sll_s_spline_fem_mixedmass_line, &
51  sll_s_spline_fem_multiply_mass, &
52  sll_s_spline_fem_multiply_massmixed, &
53  sll_s_spline_fem_compute_mass_eig, &
54  sll_s_multiply_g, &
55  sll_s_multiply_gt, &
56  sll_s_multiply_c, &
57  sll_s_multiply_ct
58 
59  use sll_m_spline_fem_utilities_3d, only: &
60  sll_s_spline_fem_mass3d
61 
62  use sll_m_spline_fem_utilities_sparse, only : &
63  sll_s_spline_fem_mass1d
64 
65  implicit none
66 
67  public :: &
68  sll_t_maxwell_3d_fem_fft
69 
70  private
71 
72  type, extends(sll_c_maxwell_3d_base) :: sll_t_maxwell_3d_fem_fft
73 
74  sll_real64, allocatable :: work(:)
75  sll_real64, allocatable :: work2(:)
76  sll_real64, allocatable :: work3d(:,:,:)
77  sll_real64, allocatable :: work_d1(:)
78  sll_real64, allocatable :: work_d2_in(:)
79  sll_real64, allocatable :: work_d2_out(:)
80  sll_real64, allocatable :: work_d3_in(:)
81  sll_real64, allocatable :: work_d3_out(:)
82 
83  sll_real64, allocatable :: mass_line0_1(:)
84  sll_real64, allocatable :: mass_line0_2(:)
85  sll_real64, allocatable :: mass_line0_3(:)
86  sll_real64, allocatable :: mass_line1_1(:)
87  sll_real64, allocatable :: mass_line1_2(:)
88  sll_real64, allocatable :: mass_line1_3(:)
89  sll_real64, allocatable :: mass_line_mixed_1(:)
90  sll_real64, allocatable :: mass_line_mixed_2(:)
91  sll_real64, allocatable :: mass_line_mixed_3(:)
92  type(sll_t_linear_solver_spline_mass_fft) :: inverse_mass_0
93  type(sll_t_linear_solver_spline_mass_fft) :: inverse_mass_1(3)
94  type(sll_t_linear_solver_spline_mass_fft) :: inverse_mass_2(3)
95  type(sll_t_poisson_3d_fem_fft ) :: poisson_fft
96 
97  logical :: adiabatic_electrons = .false.
98 
99  type(sll_t_matrix_csr) :: mass0
100  type(sll_t_matrix_csr) :: mass1d(3,3)
101  type(sll_t_linear_operator_kron) :: mass1(3)
102  type(sll_t_linear_operator_kron) :: mass2(3)
103  type( sll_t_linear_operator_maxwell_eb_schur ) :: linear_op_schur_eb
104 
105  contains
106  procedure :: &
107  compute_e_from_b => sll_s_compute_e_from_b_3d_fem_fft
108  procedure :: &
109  compute_b_from_e => sll_s_compute_b_from_e_3d_fem_fft
110  procedure :: &
111  compute_curl_part => sll_s_compute_curl_part_3d_fem_fft
112  procedure :: &
113  compute_e_from_rho => sll_s_compute_e_from_rho_3d_fem_fft
114  procedure :: &
115  compute_rho_from_e => sll_s_compute_rho_from_e_3d_fem_fft
116  procedure :: &
117  compute_e_from_j => sll_s_compute_e_from_j_3d_fem_fft
118  procedure :: &
119  compute_phi_from_rho => sll_s_compute_phi_from_rho_3d_fem_fft
120  procedure :: &
121  compute_phi_from_j => sll_s_compute_phi_from_j_3d_fem_fft
122  procedure :: &
123  compute_rhs_from_function => sll_s_compute_rhs_fem_fft
124  procedure :: &
125  l2projection => l2projection_3d_fem_fft
126  procedure :: &
127  l2norm_squared => l2norm_squared_3d_fem_fft
128  procedure :: &
129  inner_product => inner_product_3d_fem_fft
130  procedure :: &
131  init => init_3d_fem_fft
132  procedure :: &
133  free => free_3d_fem_fft
134  procedure :: &
135  multiply_g
136  procedure :: &
137  multiply_gt
138  procedure :: &
139  multiply_c
140  procedure :: &
141  multiply_ct
142  procedure :: &
143  multiply_mass => multiply_mass_3d_fem_fft
144  procedure :: &
145  multiply_mass_inverse => multiply_mass_inverse_all
146 
147  end type sll_t_maxwell_3d_fem_fft
148 
149 contains
150 
151 
153  subroutine sll_s_compute_e_from_b_3d_fem_fft(self, delta_t, field_in, field_out)
154  class(sll_t_maxwell_3d_fem_fft) :: self
155  sll_real64, intent(in) :: delta_t
156  sll_real64, intent(in) :: field_in(:)
157  sll_real64, intent(inout) :: field_out(:)
158  ! local variables
159  sll_real64 :: coef
160  sll_int32 :: comp, istart, iend
161 
162  call multiply_mass_2form( self, field_in, self%work )
163 
164  call multiply_ct(self, self%work, self%work2)
165 
166  do comp=1,3
167  istart = 1+(comp-1)*self%n_total
168  iend = comp*self%n_total
169  call self%inverse_mass_1(comp)%solve( self%work2(istart:iend), self%work(istart:iend) )
170  end do
171  ! Update b from self value
172  coef = delta_t
173  field_out = field_out + coef*self%work
174  end subroutine sll_s_compute_e_from_b_3d_fem_fft
175 
176 
179  subroutine sll_s_compute_b_from_e_3d_fem_fft(self, delta_t, field_in, field_out)
180  class(sll_t_maxwell_3d_fem_fft) :: self
181  sll_real64, intent( in ) :: delta_t
182  sll_real64, intent( in ) :: field_in(:)
183  sll_real64, intent( inout ) :: field_out(:)
184 
185  call multiply_c(self, field_in, self%work)
186 
187  field_out = field_out - delta_t * self%work
188 
189  end subroutine sll_s_compute_b_from_e_3d_fem_fft
190 
191 
193  subroutine sll_s_compute_curl_part_3d_fem_fft( self, delta_t, efield, bfield, betar )
194  class(sll_t_maxwell_3d_fem_fft) :: self
195  sll_real64, intent( in ) :: delta_t
196  sll_real64, intent( inout ) :: efield(:)
197  sll_real64, intent( inout ) :: bfield(:)
198  sll_real64, optional :: betar
199  !local variables
200  sll_real64 :: factor
201 
202  if( present(betar) ) then
203  factor = betar
204  else
205  factor = 1._f64
206  end if
207 
208  ! Compute C^T M2 b
209  call multiply_mass_2form( self, bfield, self%work )
210  call self%multiply_ct( self%work, self%work2 )
211 
212  self%linear_op_schur_eb%factor = -delta_t**2*factor*0.25_f64
213  call self%linear_op_schur_eb%dot( efield, self%work )
214  self%work = self%work + delta_t*factor*self%work2
215 
216  ! Save efield dofs from previous time step for B field update
217  self%work2 = efield
218 
219  ! Invert Schur complement matrix
220  self%linear_op_schur_eb%factor = delta_t**2*factor*0.25_f64
221  call self%linear_op_schur_eb%dot_inverse( self%work, efield )
222 
223  ! Update B field
224  self%work2 = self%work2 + efield
225  call self%compute_b_from_e( delta_t*0.5_f64, self%work2, bfield )
226 
227  end subroutine sll_s_compute_curl_part_3d_fem_fft
228 
229 
231  subroutine sll_s_compute_e_from_rho_3d_fem_fft(self, field_in, field_out )
232  class(sll_t_maxwell_3d_fem_fft) :: self
233  sll_real64, intent( in ) :: field_in(:)
234  sll_real64, intent( out ) :: field_out(:)
235 
236 
237  call self%poisson_fft%compute_e_from_rho( field_in, field_out )
238 
239  end subroutine sll_s_compute_e_from_rho_3d_fem_fft
240 
241 
243  subroutine sll_s_compute_rho_from_e_3d_fem_fft(self, field_in, field_out )
244  class(sll_t_maxwell_3d_fem_fft) :: self
245  sll_real64, intent( in ) :: field_in(:)
246  sll_real64, intent( out ) :: field_out(:)
247 
248  call multiply_mass_1form( self, field_in, self%work )
249 
250  call multiply_gt( self, self%work, field_out )
251  field_out = - field_out
252  end subroutine sll_s_compute_rho_from_e_3d_fem_fft
253 
254 
256  subroutine sll_s_compute_e_from_j_3d_fem_fft(self, current, E, component)
257  class(sll_t_maxwell_3d_fem_fft) :: self
258  sll_real64, intent( in ) :: current(:)
259  sll_real64, intent( inout ) :: e(:)
260  sll_int32, intent( in ), optional :: component
261 
262  if(present(component) ) then
263  call self%inverse_mass_1(component)%solve( current, self%work(1:self%n_total) )
264  e = e - self%work(1:self%n_total)
265  else
266  call self%inverse_mass_1(1)%solve( current(1:self%n_total), self%work(1:self%n_total) )
267  call self%inverse_mass_1(2)%solve( current(self%n_total+1:2*self%n_total), self%work(self%n_total+1:2*self%n_total) )
268  call self%inverse_mass_1(3)%solve( current(2*self%n_total+1:3*self%n_total), self%work(2*self%n_total+1:3*self%n_total) )
269  e = e - self%work
270  end if
271 
272  end subroutine sll_s_compute_e_from_j_3d_fem_fft
273 
274 
276  subroutine sll_s_compute_phi_from_rho_3d_fem_fft( self, field_in, field_out, efield_dofs )
277  class(sll_t_maxwell_3d_fem_fft) :: 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%inverse_mass_0%solve( field_in, field_out )
283  call self%multiply_g( field_out, efield_dofs )
284  efield_dofs = -efield_dofs
285 
286  end subroutine sll_s_compute_phi_from_rho_3d_fem_fft
287 
288 
290  subroutine sll_s_compute_phi_from_j_3d_fem_fft( self, field_in, field_out, efield_dofs )
291  class(sll_t_maxwell_3d_fem_fft) :: self
292  sll_real64, intent( in ) :: field_in(:)
293  sll_real64, intent( inout ) :: field_out(:)
294  sll_real64, intent( out ) :: efield_dofs(:)
295 
296  call self%multiply_gt( field_in, self%work(1:self%n_total) )
297  call self%inverse_mass_0%solve( self%work(1:self%n_total), self%work2(1:self%n_total) )
298  field_out = field_out + self%work2(1:self%n_total)
299 
300  call self%multiply_g( field_out, efield_dofs )
301  efield_dofs = -efield_dofs
302 
303  end subroutine sll_s_compute_phi_from_j_3d_fem_fft
304 
305 
308  subroutine sll_s_compute_rhs_fem_fft(self, form, component, coefs_dofs, func1, func2, func3 )
309  class(sll_t_maxwell_3d_fem_fft) :: self
310  sll_int32, intent( in ) :: form
311  sll_int32, intent( in ) :: component
312  sll_real64, intent( out ) :: coefs_dofs(:)
313  procedure(sll_i_function_3d_real64) :: func1
314  procedure(sll_i_function_3d_real64), optional :: func2
315  procedure(sll_i_function_3d_real64), optional :: func3
316  ! local variables
317  sll_int32 :: i1, i2, i3,j1, j2, j3, k1, k2, k3, q(3), counter
318  sll_int32 :: degree(3)
319  sll_real64 :: c(3)
320  sll_real64 :: coef
321  sll_real64, allocatable :: xw_gauss_d1(:,:), xw_gauss_d2(:,:), xw_gauss_d3(:,:)
322  sll_real64, allocatable :: bspl_d1(:,:), bspl_d2(:,:), bspl_d3(:,:)
323  sll_real64 :: scratch(maxval(self%s_deg_0+1))
324 
325  ! Define the spline degree in the 3 dimensions, depending on form and component of the form
326  if ( form == 0 ) then
327  degree = self%s_deg_0
328  elseif (form == 1 ) then
329  degree = self%s_deg_0
330  degree(component) = self%s_deg_1(component)
331  elseif( form == 2) then
332  degree = self%s_deg_1
333  degree(component) = self%s_deg_0(component)
334  elseif( form == 3) then
335  degree = self%s_deg_1
336  else
337  print*, 'Wrong form.'
338  end if
339 
340  ! take enough Gauss points so that projection is exact for splines of degree deg
341  q = degree+1
342  ! rescale on [0,1] for compatibility with B-splines
343  allocate(xw_gauss_d1(2,q(1)))
344  allocate(bspl_d1(q(1), degree(1)+1))
345  xw_gauss_d1 = sll_f_gauss_legendre_points_and_weights(q(1), 0.0_f64, 1.0_f64)
346  ! Compute bsplines at gauss_points
347  do k1 = 1, q(1)
348  call sll_s_uniform_bsplines_eval_basis(degree(1),xw_gauss_d1(1,k1), scratch(1:degree(1)+1))
349  bspl_d1(k1,:) = scratch(1:degree(1)+1)
350  end do
351 
352  allocate(xw_gauss_d2(2,q(2)))
353  allocate(bspl_d2(q(2), degree(2)+1))
354  xw_gauss_d2 = sll_f_gauss_legendre_points_and_weights(q(2), 0.0_f64, 1.0_f64)
355  ! Compute bsplines at gauss_points
356  do k2 = 1, q(2)
357  call sll_s_uniform_bsplines_eval_basis(degree(2),xw_gauss_d2(1,k2), scratch(1:degree(2)+1))
358  bspl_d2(k2,:) = scratch(1:degree(2)+1)
359  end do
360 
361  allocate(xw_gauss_d3(2,q(3)))
362  allocate(bspl_d3(q(3), degree(3)+1))
363  xw_gauss_d3 = sll_f_gauss_legendre_points_and_weights(q(3), 0.0_f64, 1.0_f64)
364  ! Compute bsplines at gauss_points
365  do k3 = 1, q(3)
366  call sll_s_uniform_bsplines_eval_basis(degree(3),xw_gauss_d3(1,k3), scratch(1:degree(3)+1))
367  bspl_d3(k3,:) = scratch(1:degree(3)+1)
368  end do
369 
370  counter = 1
371  ! Compute coefs_dofs = int f(x)N_i(x)
372  do i3 = 1, self%n_dofs(3)
373  do i2 = 1, self%n_dofs(2)
374  do i1 = 1, self%n_dofs(1)
375  coef=0.0_f64
376  ! loop over support of B spline
377  do j3 = 1, degree(3)+1
378  do j2 = 1, degree(2)+1
379  do j1 = 1, degree(1)+1
380  ! loop over Gauss points
381  do k3 = 1, q(3)
382  c(3) = self%delta_x(3)*(xw_gauss_d3(1,k3) + real(i3 + j3 - 2,f64))
383  do k2 = 1, q(2)
384  c(2) = self%delta_x(2)*(xw_gauss_d2(1,k2) + real(i2 + j2 - 2,f64))
385  do k1 = 1, q(1)
386  c(1) = self%delta_x(1)*(xw_gauss_d1(1,k1) + real(i1 + j1 - 2,f64))
387  coef = coef + xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* &
388  xw_gauss_d3(2,k3) *&
389  func1( c ) * &
390  bspl_d1(k1,degree(1)+2-j1)*&
391  bspl_d2(k2,degree(2)+2-j2)*&
392  bspl_d3(k3,degree(3)+2-j3)
393 
394  enddo
395  enddo
396  end do
397  end do
398  end do
399  end do
400  ! rescale by cell size
401  coefs_dofs(counter) = coef*self%volume
402  counter = counter+1
403  enddo
404  end do
405  end do
406 
407  end subroutine sll_s_compute_rhs_fem_fft
408 
409 
411  subroutine l2projection_3d_fem_fft(self, form, component, coefs_dofs, func1, func2, func3 )
412  class(sll_t_maxwell_3d_fem_fft) :: self
413  sll_int32, intent( in ) :: form
414  sll_int32, intent( in ) :: component
415  sll_real64, intent( out ) :: coefs_dofs(:)
416  procedure(sll_i_function_3d_real64) :: func1
417  procedure(sll_i_function_3d_real64), optional :: func2
418  procedure(sll_i_function_3d_real64), optional :: func3
419 
420  ! Compute right-hand-side
421  call sll_s_compute_rhs_fem_fft( self, form, component, self%work(1:self%n_total), func1 )
422 
423  select case( form )
424  case( 0 )
425  call self%inverse_mass_0%solve( self%work(1:self%n_total), coefs_dofs )
426  case( 1 )
427  call self%inverse_mass_1(component)%solve( self%work(1:self%n_total), coefs_dofs )
428  case( 2 )
429  call self%inverse_mass_2(component)%solve( self%work(1:self%n_total), coefs_dofs )
430  case default
431  print*, 'L2projection for', form, '-form not implemented.'
432  end select
433 
434  end subroutine l2projection_3d_fem_fft
435 
436 
438  function l2norm_squared_3d_fem_fft(self, coefs, form, component ) result ( r )
439  class(sll_t_maxwell_3d_fem_fft) :: self
440  sll_real64 :: coefs(:)
441  sll_int32 :: form
442  sll_int32 :: component
443  sll_real64 :: r
444 
445  r = inner_product_3d_fem_fft(self, coefs, coefs, form, component)
446 
447  end function l2norm_squared_3d_fem_fft
448 
449 
451  function inner_product_3d_fem_fft(self, coefs1, coefs2, form, component) result (r)
452  class(sll_t_maxwell_3d_fem_fft) :: self
453  sll_real64 :: coefs1(:)
454  sll_real64 :: coefs2(:)
455  sll_int32 :: form
456  sll_int32, optional :: component
457  sll_real64 :: r
458 
459  if ( form == 0 ) then
460  call multiply_mass_3dkron( self, self%mass_line0_1, &
461  self%mass_line0_2, self%mass_line0_3, &
462  coefs2, self%work(1:self%n_total) )
463  elseif (form == 1 ) then
464  select case(component)
465  case (1)
466  call multiply_mass_3dkron( self, self%mass_line1_1, &
467  self%mass_line0_2, self%mass_line0_3, &
468  coefs2, self%work(1:self%n_total) )
469  case(2)
470  call multiply_mass_3dkron( self, self%mass_line0_1, &
471  self%mass_line1_2, self%mass_line0_3, &
472  coefs2, self%work(1:self%n_total) )
473  case(3)
474  call multiply_mass_3dkron( self, self%mass_line0_1, &
475  self%mass_line0_2, self%mass_line1_3, &
476  coefs2, self%work(1:self%n_total) )
477  case default
478  print*, 'wrong component.'
479  end select
480  elseif( form == 2) then
481  select case(component)
482  case (1)
483  call multiply_mass_3dkron( self, self%mass_line0_1, &
484  self%mass_line1_2, self%mass_line1_3, &
485  coefs2, self%work(1:self%n_total) )
486  case(2)
487  call multiply_mass_3dkron( self, self%mass_line1_1, &
488  self%mass_line0_2, self%mass_line1_3, &
489  coefs2, self%work(1:self%n_total) )
490  case(3)
491  call multiply_mass_3dkron( self, self%mass_line1_1, &
492  self%mass_line1_2, self%mass_line0_3, &
493  coefs2, self%work(1:self%n_total) )
494  case default
495  print*, 'wrong component.'
496  end select
497  elseif( form == 3) then
498  call multiply_mass_3dkron( self, self%mass_line1_1, &
499  self%mass_line1_2, self%mass_line1_3, &
500  coefs2, self%work(1:self%n_total) )
501  else
502  print*, 'Wrong form.'
503  end if
504 
505  r = sum(coefs1*self%work(1:self%n_total))
506  !r = r*self%volume
507 
508 
509  end function inner_product_3d_fem_fft
510 
511 
513  subroutine init_3d_fem_fft( self, domain, n_dofs, s_deg_0, adiabatic_electrons, profile )
514  class(sll_t_maxwell_3d_fem_fft), intent(out) :: self
515  sll_real64, intent(in) :: domain(3,2)
516  sll_int32, intent(in) :: n_dofs(3)
517  sll_int32, intent(in) :: s_deg_0(3)
518  logical, intent(in), optional :: adiabatic_electrons
519  type(sll_t_profile_functions), intent(in), optional :: profile
520  ! local variables
521  sll_int32 :: j
522  sll_real64, allocatable :: eig_values_mass_0_1(:)
523  sll_real64, allocatable :: eig_values_mass_0_2(:)
524  sll_real64, allocatable :: eig_values_mass_0_3(:)
525  sll_real64, allocatable :: eig_values_mass_1_1(:)
526  sll_real64, allocatable :: eig_values_mass_1_2(:)
527  sll_real64, allocatable :: eig_values_mass_1_3(:)
528  sll_real64, allocatable :: inv_eig_values_mass_0_1(:)
529  sll_real64, allocatable :: inv_eig_values_mass_0_2(:)
530  sll_real64, allocatable :: inv_eig_values_mass_0_3(:)
531  sll_real64, allocatable :: inv_eig_values_mass_1_1(:)
532  sll_real64, allocatable :: inv_eig_values_mass_1_2(:)
533  sll_real64, allocatable :: inv_eig_values_mass_1_3(:)
534  sll_real64 :: factor
535 
536  self%n_cells = n_dofs
537  self%n_dofs = n_dofs
538  self%n_total = product(n_dofs)
539  self%n_total0 = self%n_total
540  self%n_total1 = self%n_total
541 
542  self%Lx = domain(:,2) - domain(:,1)
543  self%delta_x = self%Lx / real(n_dofs,f64)
544  self%s_deg_0 = s_deg_0
545  self%s_deg_1 = s_deg_0 - 1
546 
547  self%volume = product(self%delta_x)
548 
549  if( present( adiabatic_electrons ) ) then
550  self%adiabatic_electrons = adiabatic_electrons
551  end if
552 
553  if( present( profile ) ) then
554  self%profile = profile
555  end if
556 
557  ! Allocate scratch data
558  allocate( self%work3d(n_dofs(1), n_dofs(2), n_dofs(3)) )
559  allocate( self%work(self%n_total*3) )
560  allocate( self%work2(self%n_total*3) )
561  allocate( self%work_d1( n_dofs(1) ) )
562  allocate( self%work_d2_in( n_dofs(2) ) )
563  allocate( self%work_d2_out( n_dofs(2) ) )
564  allocate( self%work_d3_in( n_dofs(3) ) )
565  allocate( self%work_d3_out( n_dofs(3) ) )
566 
567 
568  ! Sparse matrices
569  ! Assemble the mass matrices
570  ! First assemble a mass line for both degrees
571  allocate( self%mass_line0_1(s_deg_0(1)+1) )
572  allocate( self%mass_line1_1(s_deg_0(1)) )
573  allocate( self%mass_line0_2(s_deg_0(2)+1) )
574  allocate( self%mass_line1_2(s_deg_0(2)) )
575  allocate( self%mass_line0_3(s_deg_0(3)+1) )
576  allocate( self%mass_line1_3(s_deg_0(3)) )
577  allocate( self%mass_line_mixed_1(2*s_deg_0(1)) )
578  allocate( self%mass_line_mixed_2(2*s_deg_0(2)) )
579  allocate( self%mass_line_mixed_3(2*s_deg_0(3)) )
580 
581 
582  call sll_s_spline_fem_mass_line ( self%s_deg_0(1), self%mass_line0_1 )
583  call sll_s_spline_fem_mass_line ( self%s_deg_1(1), self%mass_line1_1 )
584  call sll_s_spline_fem_mixedmass_line ( self%s_deg_0(1), self%mass_line_mixed_1 )
585 
586  call sll_s_spline_fem_mass_line ( self%s_deg_0(2), self%mass_line0_2 )
587  call sll_s_spline_fem_mass_line ( self%s_deg_1(2), self%mass_line1_2 )
588  call sll_s_spline_fem_mixedmass_line ( self%s_deg_0(2), self%mass_line_mixed_2 )
589 
590  call sll_s_spline_fem_mass_line ( self%s_deg_0(3), self%mass_line0_3 )
591  call sll_s_spline_fem_mass_line ( self%s_deg_1(3), self%mass_line1_3 )
592  call sll_s_spline_fem_mixedmass_line ( self%s_deg_0(3), self%mass_line_mixed_3 )
593 
594  self%mass_line0_1 = self%delta_x(1) * self%mass_line0_1
595  self%mass_line1_1 = self%delta_x(1) * self%mass_line1_1
596  self%mass_line_mixed_1 = self%delta_x(1) * self%mass_line_mixed_1
597 
598  self%mass_line0_2 = self%delta_x(2) * self%mass_line0_2
599  self%mass_line1_2 = self%delta_x(2) * self%mass_line1_2
600  self%mass_line_mixed_2 = self%delta_x(2) * self%mass_line_mixed_2
601 
602  self%mass_line0_3 = self%delta_x(3) * self%mass_line0_3
603  self%mass_line1_3 = self%delta_x(3) * self%mass_line1_3
604  self%mass_line_mixed_3 = self%delta_x(3) * self%mass_line_mixed_3
605 
606  allocate( eig_values_mass_0_1( n_dofs(1) ) )
607  allocate( eig_values_mass_0_2( n_dofs(2) ) )
608  allocate( eig_values_mass_0_3( n_dofs(3) ) )
609  allocate( eig_values_mass_1_1( n_dofs(1) ) )
610  allocate( eig_values_mass_1_2( n_dofs(2) ) )
611  allocate( eig_values_mass_1_3( n_dofs(3) ) )
612  allocate( inv_eig_values_mass_0_1( n_dofs(1) ) )
613  allocate( inv_eig_values_mass_0_2( n_dofs(2) ) )
614  allocate( inv_eig_values_mass_0_3( n_dofs(3) ) )
615  allocate( inv_eig_values_mass_1_1( n_dofs(1) ) )
616  allocate( inv_eig_values_mass_1_2( n_dofs(2) ) )
617  allocate( inv_eig_values_mass_1_3( n_dofs(3) ) )
618 
619  call sll_s_spline_fem_compute_mass_eig( n_dofs(1), self%s_deg_0(1), self%mass_line0_1, &
620  eig_values_mass_0_1 )
621  call sll_s_spline_fem_compute_mass_eig( n_dofs(2), self%s_deg_0(2), self%mass_line0_2, &
622  eig_values_mass_0_2 )
623  call sll_s_spline_fem_compute_mass_eig( n_dofs(3), self%s_deg_0(3), self%mass_line0_3, &
624  eig_values_mass_0_3 )
625  call sll_s_spline_fem_compute_mass_eig( n_dofs(1), self%s_deg_1(1), self%mass_line1_1, &
626  eig_values_mass_1_1 )
627  call sll_s_spline_fem_compute_mass_eig( n_dofs(2), self%s_deg_1(2), self%mass_line1_2, &
628  eig_values_mass_1_2 )
629  call sll_s_spline_fem_compute_mass_eig( n_dofs(3), self%s_deg_1(3), self%mass_line1_3, &
630  eig_values_mass_1_3 )
631 
632  inv_eig_values_mass_0_1(1) = 1._f64/eig_values_mass_0_1(1)
633  inv_eig_values_mass_0_1(n_dofs(1)/2+1) = 1._f64/eig_values_mass_0_1(n_dofs(1)/2+1)
634  inv_eig_values_mass_1_1(1) = 1._f64/eig_values_mass_1_1(1)
635  inv_eig_values_mass_1_1(n_dofs(1)/2+1) = 1._f64/eig_values_mass_1_1(n_dofs(1)/2+1)
636  do j = 2, n_dofs(1)/2
637  inv_eig_values_mass_0_1(j) = 1._f64/eig_values_mass_0_1(j)
638  inv_eig_values_mass_0_1(n_dofs(1)+2-j) = 1._f64/eig_values_mass_0_1(n_dofs(1)+2-j)
639  inv_eig_values_mass_1_1(j) = 1._f64/eig_values_mass_1_1(j)
640  inv_eig_values_mass_1_1(n_dofs(1)+2-j) = 1._f64/eig_values_mass_1_1(n_dofs(1)+2-j)
641  end do
642 
643  inv_eig_values_mass_0_2(1) = 1._f64/eig_values_mass_0_2(1)
644  inv_eig_values_mass_0_2(n_dofs(2)/2+1) = 1._f64/eig_values_mass_0_2(n_dofs(2)/2+1)
645  inv_eig_values_mass_1_2(1) = 1._f64/eig_values_mass_1_2(1)
646  inv_eig_values_mass_1_2(n_dofs(2)/2+1) = 1._f64/eig_values_mass_1_2(n_dofs(2)/2+1)
647  do j = 2, n_dofs(2)/2
648  inv_eig_values_mass_0_2(j) = 1._f64/eig_values_mass_0_2(j)
649  inv_eig_values_mass_0_2(n_dofs(2)+2-j) = 1._f64/eig_values_mass_0_2(n_dofs(2)+2-j)
650  inv_eig_values_mass_1_2(j) = 1._f64/eig_values_mass_1_2(j)
651  inv_eig_values_mass_1_2(n_dofs(2)+2-j) = 1._f64/eig_values_mass_1_2(n_dofs(2)+2-j)
652  end do
653 
654  inv_eig_values_mass_0_3(1) = 1._f64/eig_values_mass_0_3(1)
655  inv_eig_values_mass_0_3(n_dofs(3)/2+1) = 1._f64/eig_values_mass_0_3(n_dofs(3)/2+1)
656  inv_eig_values_mass_1_3(1) = 1._f64/eig_values_mass_1_3(1)
657  inv_eig_values_mass_1_3(n_dofs(3)/2+1) = 1._f64/eig_values_mass_1_3(n_dofs(3)/2+1)
658  do j = 2, n_dofs(3)/2
659  inv_eig_values_mass_0_3(j) = 1._f64/eig_values_mass_0_3(j)
660  inv_eig_values_mass_0_3(n_dofs(3)+2-j) = 1._f64/eig_values_mass_0_3(n_dofs(3)+2-j)
661  inv_eig_values_mass_1_3(j) = 1._f64/eig_values_mass_1_3(j)
662  inv_eig_values_mass_1_3(n_dofs(3)+2-j) = 1._f64/eig_values_mass_1_3(n_dofs(3)+2-j)
663  end do
664 
665  if(self%adiabatic_electrons) then
666  factor = self%profile%T_e( 0._f64 )/self%profile%rho_0( 0._f64 )
667  call self%inverse_mass_0%create( n_dofs, inv_eig_values_mass_0_1*factor, inv_eig_values_mass_0_2*factor, inv_eig_values_mass_0_3*factor )
668  else
669  call self%inverse_mass_0%create( n_dofs, inv_eig_values_mass_0_1, inv_eig_values_mass_0_2, inv_eig_values_mass_0_3 )
670  end if
671 
672  call self%inverse_mass_1(1)%create( n_dofs, inv_eig_values_mass_1_1, inv_eig_values_mass_0_2, inv_eig_values_mass_0_3 )
673  call self%inverse_mass_1(2)%create( n_dofs, inv_eig_values_mass_0_1, inv_eig_values_mass_1_2, inv_eig_values_mass_0_3 )
674  call self%inverse_mass_1(3)%create( n_dofs, inv_eig_values_mass_0_1, inv_eig_values_mass_0_2, inv_eig_values_mass_1_3 )
675 
676  call self%inverse_mass_2(1)%create( n_dofs, inv_eig_values_mass_0_1, inv_eig_values_mass_1_2, inv_eig_values_mass_1_3 )
677  call self%inverse_mass_2(2)%create( n_dofs, inv_eig_values_mass_1_1, inv_eig_values_mass_0_2, inv_eig_values_mass_1_3 )
678  call self%inverse_mass_2(3)%create( n_dofs, inv_eig_values_mass_1_1, inv_eig_values_mass_1_2, inv_eig_values_mass_0_3 )
679 
680 
681  ! Poisson solver based on fft inversion
682  call self%poisson_fft%init( self%n_dofs, self%s_deg_0, self%delta_x )
683 
684  ! Next put together the 1d parts of the 3d Kronecker product
685  call sll_s_spline_fem_mass1d( self%n_dofs(1), self%s_deg_0(1), self%mass_line0_1, self%mass1d(1,1) )
686  call sll_s_spline_fem_mass1d( self%n_dofs(2), self%s_deg_0(2), self%mass_line0_2, self%mass1d(1,2) )
687  call sll_s_spline_fem_mass1d( self%n_dofs(3), self%s_deg_0(3), self%mass_line0_3, self%mass1d(1,3) )
688 
689  call sll_s_spline_fem_mass1d( self%n_dofs(1), self%s_deg_1(1), self%mass_line1_1, self%mass1d(2,1) )
690  call sll_s_spline_fem_mass1d( self%n_dofs(2), self%s_deg_1(2), self%mass_line1_2, self%mass1d(2,2) )
691  call sll_s_spline_fem_mass1d( self%n_dofs(3), self%s_deg_1(3), self%mass_line1_3, self%mass1d(2,3) )
692 
693  ! Only for Schur complement eb solver
694  call self%mass1(1)%create( linop_a=self%mass1d(1,3), &
695  linop_b=self%mass1d(1,2), &
696  linop_c=self%mass1d(2,1) )
697  call self%mass1(2)%create( linop_a=self%mass1d(1,3), &
698  linop_b=self%mass1d(2,2), &
699  linop_c=self%mass1d(1,1) )
700  call self%mass1(3)%create( linop_a=self%mass1d(2,3), &
701  linop_b=self%mass1d(1,2), &
702  linop_c=self%mass1d(1,1))
703 
704  call self%mass2(1)%create( linop_a=self%mass1d(2,3), &
705  linop_b=self%mass1d(2,2), &
706  linop_c=self%mass1d(1,1) )
707  call self%mass2(2)%create( linop_a=self%mass1d(2,3), &
708  linop_b=self%mass1d(1,2), &
709  linop_c=self%mass1d(2,1) )
710  call self%mass2(3)%create( linop_a=self%mass1d(1,3), &
711  linop_b=self%mass1d(2,2), &
712  linop_c=self%mass1d(2,1))
713 
714 
715  call self%linear_op_schur_eb%create( eig_values_mass_0_1, eig_values_mass_0_2, eig_values_mass_0_3, eig_values_mass_1_1, eig_values_mass_1_2, eig_values_mass_1_3, self%n_dofs, self%delta_x )
716 
717  end subroutine init_3d_fem_fft
718 
719 
721  subroutine free_3d_fem_fft(self)
722  class(sll_t_maxwell_3d_fem_fft) :: self
723  !local variable
724  sll_int32 :: j
725 
726  call self%inverse_mass_0%free()
727  do j=1, 3
728  call self%inverse_mass_1(j)%free()
729  call self%inverse_mass_2(j)%free()
730  end do
731  call self%poisson_fft%free()
732  deallocate( self%work )
733  deallocate( self%work3d )
734  deallocate( self%work2 )
735  deallocate( self%work_d1 )
736  deallocate( self%work_d2_in )
737  deallocate( self%work_d2_out )
738  deallocate( self%work_d3_in )
739  deallocate( self%work_d3_out )
740 
741  deallocate( self%mass_line0_1 )
742  deallocate( self%mass_line1_1 )
743  deallocate( self%mass_line0_2 )
744  deallocate( self%mass_line1_2 )
745  deallocate( self%mass_line0_3 )
746  deallocate( self%mass_line1_3 )
747  deallocate( self%mass_line_mixed_1 )
748  deallocate( self%mass_line_mixed_2 )
749  deallocate( self%mass_line_mixed_3 )
750 
751  end subroutine free_3d_fem_fft
752 
753 
755  subroutine multiply_g( self, field_in, field_out )
756  class(sll_t_maxwell_3d_fem_fft) :: self
757  sll_real64, intent( in ) :: field_in(:)
758  sll_real64, intent( out ) :: field_out(:)
759 
760  call sll_s_multiply_g(self%n_dofs, self%delta_x, field_in, field_out)
761 
762  end subroutine multiply_g
763 
764 
766  subroutine multiply_gt(self, field_in, field_out)
767  class(sll_t_maxwell_3d_fem_fft) :: self
768  sll_real64, intent( in ) :: field_in(:)
769  sll_real64, intent( out ) :: field_out(:)
770 
771  call sll_s_multiply_gt(self%n_dofs, self%delta_x, field_in, field_out)
772 
773  end subroutine multiply_gt
774 
775 
777  subroutine multiply_c(self, field_in, field_out)
778  class(sll_t_maxwell_3d_fem_fft) :: self
779  sll_real64, intent( in ) :: field_in(:)
780  sll_real64, intent( out ) :: field_out(:)
781 
782  call sll_s_multiply_c(self%n_dofs, self%delta_x, field_in, field_out)
783 
784  end subroutine multiply_c
785 
786 
788  subroutine multiply_ct(self, field_in, field_out)
789  class(sll_t_maxwell_3d_fem_fft) :: self
790  sll_real64, intent( in ) :: field_in(:)
791  sll_real64, intent( out ) :: field_out(:)
792 
793  call sll_s_multiply_ct(self%n_dofs, self%delta_x, field_in, field_out)
794 
795  end subroutine multiply_ct
796 
797 
799  subroutine multiply_mass_3d_fem_fft( self, deg, coefs_in, coefs_out )
800  class(sll_t_maxwell_3d_fem_fft) :: self
801  sll_int32, intent( in ) :: deg(:)
802  sll_real64, intent( in ) :: coefs_in(:)
803  sll_real64, intent( out ) :: coefs_out(:)
804 
805  if( size(deg) ==1 )then
806  select case(deg(1))
807  case(0)
808  call multiply_mass_3dkron( self, self%mass_line0_1, &
809  self%mass_line0_2, self%mass_line0_3, &
810  coefs_in, coefs_out )
811  case(1)
812  call multiply_mass_1form( self, coefs_in, coefs_out )
813  case(2)
814  call multiply_mass_2form( self, coefs_in, coefs_out )
815  case default
816  sll_error('maxwell_3d_fem_fft','multiply mass for other form not yet implemented')
817  end select
818  else if( size(deg) == 3 ) then
819  call multiply_mass_all( self, deg, coefs_in, coefs_out )
820  end if
821 
822  end subroutine multiply_mass_3d_fem_fft
823 
824 
826  subroutine multiply_mass_1form( self, coefs_in, coefs_out )
827  class(sll_t_maxwell_3d_fem_fft), intent( inout ) :: self
828  sll_real64, intent( in ) :: coefs_in(:)
829  sll_real64, intent( out ) :: coefs_out(:)
830  !local variables
831  sll_int32:: iend, istart
832 
833  istart = 1
834  iend = self%n_total
835  call multiply_mass_3dkron( self, self%mass_line1_1, &
836  self%mass_line0_2, self%mass_line0_3, &
837  coefs_in(istart:iend), coefs_out(istart:iend) )
838  istart = iend+1
839  iend = iend + self%n_total
840  call multiply_mass_3dkron( self, self%mass_line0_1, &
841  self%mass_line1_2, self%mass_line0_3, &
842  coefs_in(istart:iend), coefs_out(istart:iend) )
843  istart = iend+1
844  iend = iend + self%n_total
845  call multiply_mass_3dkron( self, self%mass_line0_1, &
846  self%mass_line0_2, self%mass_line1_3, &
847  coefs_in(istart:iend), coefs_out(istart:iend) )
848 
849  end subroutine multiply_mass_1form
850 
851 
853  subroutine multiply_mass_2form( self, coefs_in, coefs_out )
854  class(sll_t_maxwell_3d_fem_fft), intent( inout ) :: self
855  sll_real64, intent(in) :: coefs_in(:)
856  sll_real64, intent(out) :: coefs_out(:)
857  !local variables
858  sll_int32:: iend, istart
859 
860  istart = 1
861  iend = self%n_total
862  call multiply_mass_3dkron( self, self%mass_line0_1, &
863  self%mass_line1_2, self%mass_line1_3, &
864  coefs_in(istart:iend), coefs_out(istart:iend) )
865  istart = iend+1
866  iend = iend + self%n_total
867  call multiply_mass_3dkron( self, self%mass_line1_1, &
868  self%mass_line0_2, self%mass_line1_3, &
869  coefs_in(istart:iend), coefs_out(istart:iend) )
870  istart = iend+1
871  iend = iend + self%n_total
872  call multiply_mass_3dkron( self, self%mass_line1_1, &
873  self%mass_line1_2, self%mass_line0_3, &
874  coefs_in(istart:iend), coefs_out(istart:iend) )
875 
876  end subroutine multiply_mass_2form
877 
878 
880  subroutine multiply_mass_3dkron( self, mass_line_1, mass_line_2, mass_line_3, coefs_in, coefs_out )
881  class(sll_t_maxwell_3d_fem_fft), intent( inout ) :: self
882  !sll_int32, intent( in ) :: deg(3) !< \a deg(i) specifies the degree of the 1d mass matrix in dimension \a i (Note: 1 for 0-form, 2 for 1-form, 3 for 0-1-form mix)
883  sll_real64, intent(in) :: mass_line_1(:)
884  sll_real64, intent(in) :: mass_line_2(:)
885  sll_real64, intent(in) :: mass_line_3(:)
886  sll_real64, intent(in) :: coefs_in(:)
887  sll_real64, intent(out) :: coefs_out(:)
888  ! Local variables
889  sll_int32 :: i,j,k,istart,iend
890  sll_int32 :: deg(3)
891 
892  deg(1) = size(mass_line_1)-1
893  deg(2) = size(mass_line_2)-1
894  deg(3) = size(mass_line_3)-1
895 
896  istart = 1
897  iend = self%n_dofs(1)
898  do k=1,self%n_dofs(3)
899  do j=1,self%n_dofs(2)
900  !print*, coefs_in(istart:iend)
901  call sll_s_spline_fem_multiply_mass ( self%n_dofs(1), deg(1), &
902  mass_line_1, coefs_in(istart:iend), self%work_d1 )
903  !print*, self%mass_line_1(:,1)
904  !print*, self%work_d1
905  !stop
906  self%work3d(:,j,k) = self%work_d1
907  istart = iend+1
908  iend = iend + self%n_dofs(1)
909  end do
910  end do
911 
912  do k=1,self%n_dofs(3)
913  do i =1,self%n_dofs(1)
914  self%work_d2_in = self%work3d(i,:,k)
915  call sll_s_spline_fem_multiply_mass ( self%n_dofs(2), deg(2), &
916  mass_line_2, self%work_d2_in, self%work_d2_out )
917  self%work3d(i,:,k) = self%work_d2_out
918  end do
919  end do
920 
921  istart = 1
922  do j=1,self%n_dofs(2)
923  do i =1,self%n_dofs(1)
924  self%work_d3_in = self%work3d(i,j,:)
925  call sll_s_spline_fem_multiply_mass ( self%n_dofs(3), deg(3), &
926  mass_line_3, self%work_d3_in, self%work_d3_out )
927 
928  do k=1,self%n_dofs(3)
929  coefs_out(istart+(k-1)*self%n_dofs(1)*self%n_dofs(2)) = self%work_d3_out(k)
930  end do
931  istart = istart +1
932  end do
933  end do
934 
935  end subroutine multiply_mass_3dkron
936 
937 
939  subroutine multiply_mass_all( self, deg, coefs_in, coefs_out )
940  class(sll_t_maxwell_3d_fem_fft) :: self
941  sll_int32, intent( in ) :: deg(:)
942  sll_real64, intent( in ) :: coefs_in(:)
943  sll_real64, intent( out ) :: coefs_out(:)
944  ! Local variables
945  sll_int32 :: i,j,k,istart,iend
946 
947  istart = 1
948  iend = self%n_dofs(1)
949  do k=1,self%n_dofs(3)
950  do j=1,self%n_dofs(2)
951  select case ( deg(1) )
952  case( 1 )
953  call sll_s_spline_fem_multiply_mass ( self%n_dofs(1), self%s_deg_0(1), &
954  self%mass_line0_1, coefs_in(istart:iend), self%work_d1 )
955  case( 2 )
956  call sll_s_spline_fem_multiply_mass ( self%n_dofs(1), self%s_deg_1(1), &
957  self%mass_line1_1, coefs_in(istart:iend), self%work_d1 )
958  case ( 3 )
959  call sll_s_spline_fem_multiply_massmixed ( self%n_dofs(1), self%s_deg_0(1), &
960  self%mass_line_mixed_1, coefs_in(istart:iend), self%work_d1 )
961  case default
962  print*, 'multiply_mass is not implemented for that spline degree.'
963  end select
964  self%work3d(:,j,k) = self%work_d1
965  istart = iend+1
966  iend = iend + self%n_dofs(1)
967  end do
968  end do
969 
970  do k=1,self%n_dofs(3)
971  do i =1,self%n_dofs(1)
972  self%work_d2_in = self%work3d(i,:,k)
973  select case ( deg(2) )
974  case( 1 )
975  call sll_s_spline_fem_multiply_mass ( self%n_dofs(2), self%s_deg_0(2), &
976  self%mass_line0_2, self%work_d2_in, self%work_d2_out )
977  case( 2 )
978  call sll_s_spline_fem_multiply_mass ( self%n_dofs(2), self%s_deg_1(2), &
979  self%mass_line1_2, self%work_d2_in, self%work_d2_out )
980  case ( 3 )
981  call sll_s_spline_fem_multiply_massmixed ( self%n_dofs(2), self%s_deg_0(2), &
982  self%mass_line_mixed_2, self%work_d2_in, self%work_d2_out )
983  case default
984  print*, 'multiply_mass is not implemented for that spline degree.'
985  end select
986  self%work3d(i,:,k) = self%work_d2_out
987  end do
988  end do
989 
990  istart = 1
991  do j=1,self%n_dofs(2)
992  do i =1,self%n_dofs(1)
993  self%work_d3_in = self%work3d(i,j,:)
994  select case ( deg(3) )
995  case( 1 )
996  call sll_s_spline_fem_multiply_mass ( self%n_dofs(3), self%s_deg_0(3), &
997  self%mass_line0_3, self%work_d3_in, self%work_d3_out )
998  case( 2 )
999  call sll_s_spline_fem_multiply_mass ( self%n_dofs(3), self%s_deg_1(3), &
1000  self%mass_line1_3, self%work_d3_in, self%work_d3_out )
1001  case ( 3 )
1002  call sll_s_spline_fem_multiply_massmixed ( self%n_dofs(3), self%s_deg_0(3), &
1003  self%mass_line_mixed_3, self%work_d3_in, self%work_d3_out )
1004  case default
1005  print*, 'multiply_mass is not implemented for that spline degree.'
1006  end select
1007 
1008  do k=1,self%n_dofs(3)
1009  coefs_out(istart+(k-1)*self%n_dofs(1)*self%n_dofs(2)) = self%work_d3_out(k)
1010  end do
1011  istart = istart +1
1012  end do
1013  end do
1014 
1015 
1016  end subroutine multiply_mass_all
1017 
1018 
1020  subroutine multiply_mass_inverse_all( self, form, coefs_in, coefs_out )
1021  class(sll_t_maxwell_3d_fem_fft) :: self
1022  sll_int32, intent( in ) :: form
1023  sll_real64, intent( in ) :: coefs_in(:)
1024  sll_real64, intent( out ) :: coefs_out(:)
1025 
1026  select case( form )
1027  case( 0 )
1028  call self%inverse_mass_0%solve( coefs_in, coefs_out )
1029  case( 1 )
1030  call self%inverse_mass_1(1)%solve( coefs_in(1:self%n_total), coefs_out(1:self%n_total) )
1031  call self%inverse_mass_1(2)%solve( coefs_in(self%n_total+1:2*self%n_total), coefs_out(self%n_total+1:2*self%n_total) )
1032  call self%inverse_mass_1(3)%solve( coefs_in(2*self%n_total+1:3*self%n_total), coefs_out(2*self%n_total+1:3*self%n_total) )
1033  case( 2 )
1034  call self%inverse_mass_2(1)%solve( coefs_in(1:self%n_total), coefs_out(1:self%n_total) )
1035  call self%inverse_mass_2(2)%solve( coefs_in(self%n_total+1:2*self%n_total), coefs_out(self%n_total+1:2*self%n_total) )
1036  call self%inverse_mass_2(3)%solve( coefs_in(2*self%n_total+1:3*self%n_total), coefs_out(2*self%n_total+1:3*self%n_total) )
1037  case default
1038  print*, 'inverse_mass for', form, '-form not implemented.'
1039  end select
1040 
1041 
1042 
1043  end subroutine multiply_mass_inverse_all
1044 
1045 
1046 end module sll_m_maxwell_3d_fem_fft
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_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_kron
module for a linear operator of kronecker solver
Definition: sll_m_linear_operator_kron.F90:15
sll_m_linear_operator_maxwell_eb_schur
This linear operator implements the compatible spline FEM operator for the curl part of Maxwell's equ...
Definition: sll_m_linear_operator_maxwell_eb_schur.F90:5
sll_m_linear_solver_spline_mass_fft
Invert a circulant matrix based on diagonalization in Fourier space (3d version)
Definition: sll_m_linear_solver_spline_mass_fft.F90:10
sll_m_low_level_bsplines
Low level arbitrary degree splines.
Definition: sll_m_low_level_bsplines.F90:27
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_fem_fft
Module interface to solve Maxwell's equations in 3D The linear systems are solved based on FFT diagno...
Definition: sll_m_maxwell_3d_fem_fft.F90:12
sll_m_maxwell_3d_fem_fft::sll_s_compute_e_from_rho_3d_fem_fft
subroutine sll_s_compute_e_from_rho_3d_fem_fft(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_fem_fft.F90:232
sll_m_maxwell_3d_fem_fft::sll_s_compute_b_from_e_3d_fem_fft
subroutine sll_s_compute_b_from_e_3d_fem_fft(self, delta_t, field_in, field_out)
Compute Bz from Ey using strong 1D Faraday equation for spline coefficients $B_z^{new}(x_j) = B_z^{ol...
Definition: sll_m_maxwell_3d_fem_fft.F90:180
sll_m_maxwell_3d_fem_fft::sll_s_compute_e_from_b_3d_fem_fft
subroutine sll_s_compute_e_from_b_3d_fem_fft(self, delta_t, field_in, field_out)
compute Ey from Bz using weak Ampere formulation
Definition: sll_m_maxwell_3d_fem_fft.F90:154
sll_m_maxwell_3d_fem_fft::multiply_mass_2form
subroutine multiply_mass_2form(self, coefs_in, coefs_out)
Helper function for multiply_mass.
Definition: sll_m_maxwell_3d_fem_fft.F90:854
sll_m_maxwell_3d_fem_fft::free_3d_fem_fft
subroutine free_3d_fem_fft(self)
Finalization.
Definition: sll_m_maxwell_3d_fem_fft.F90:722
sll_m_maxwell_3d_fem_fft::multiply_ct
subroutine multiply_ct(self, field_in, field_out)
Multiply by transpose of discrete curl matrix.
Definition: sll_m_maxwell_3d_fem_fft.F90:789
sll_m_maxwell_3d_fem_fft::sll_s_compute_rhs_fem_fft
subroutine sll_s_compute_rhs_fem_fft(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_fem_fft.F90:309
sll_m_maxwell_3d_fem_fft::sll_s_compute_phi_from_rho_3d_fem_fft
subroutine sll_s_compute_phi_from_rho_3d_fem_fft(self, field_in, field_out, efield_dofs)
Compute phi from rho_i integrated over the time interval.
Definition: sll_m_maxwell_3d_fem_fft.F90:277
sll_m_maxwell_3d_fem_fft::sll_s_compute_rho_from_e_3d_fem_fft
subroutine sll_s_compute_rho_from_e_3d_fem_fft(self, field_in, field_out)
compute rho from e using weak Gauss law ( rho = G^T M_1 e )
Definition: sll_m_maxwell_3d_fem_fft.F90:244
sll_m_maxwell_3d_fem_fft::multiply_g
subroutine multiply_g(self, field_in, field_out)
Multiply by dicrete gradient matrix.
Definition: sll_m_maxwell_3d_fem_fft.F90:756
sll_m_maxwell_3d_fem_fft::sll_s_compute_e_from_j_3d_fem_fft
subroutine sll_s_compute_e_from_j_3d_fem_fft(self, current, E, component)
Compute E_i from j_i integrated over the time interval using weak Ampere formulation.
Definition: sll_m_maxwell_3d_fem_fft.F90:257
sll_m_maxwell_3d_fem_fft::l2projection_3d_fem_fft
subroutine l2projection_3d_fem_fft(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_fem_fft.F90:412
sll_m_maxwell_3d_fem_fft::multiply_c
subroutine multiply_c(self, field_in, field_out)
Multiply by discrete curl matrix.
Definition: sll_m_maxwell_3d_fem_fft.F90:778
sll_m_maxwell_3d_fem_fft::init_3d_fem_fft
subroutine init_3d_fem_fft(self, domain, n_dofs, s_deg_0, adiabatic_electrons, profile)
Initialization.
Definition: sll_m_maxwell_3d_fem_fft.F90:514
sll_m_maxwell_3d_fem_fft::multiply_gt
subroutine multiply_gt(self, field_in, field_out)
Multiply by transpose of dicrete gradient matrix.
Definition: sll_m_maxwell_3d_fem_fft.F90:767
sll_m_maxwell_3d_fem_fft::sll_s_compute_curl_part_3d_fem_fft
subroutine sll_s_compute_curl_part_3d_fem_fft(self, delta_t, efield, bfield, betar)
Solve curl part of Maxwell's equations.
Definition: sll_m_maxwell_3d_fem_fft.F90:194
sll_m_maxwell_3d_fem_fft::multiply_mass_inverse_all
subroutine multiply_mass_inverse_all(self, form, coefs_in, coefs_out)
Multiply by the inverse mass matrix.
Definition: sll_m_maxwell_3d_fem_fft.F90:1021
sll_m_maxwell_3d_fem_fft::sll_s_compute_phi_from_j_3d_fem_fft
subroutine sll_s_compute_phi_from_j_3d_fem_fft(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_fem_fft.F90:291
sll_m_maxwell_3d_fem_fft::inner_product_3d_fem_fft
real(kind=f64) function inner_product_3d_fem_fft(self, coefs1, coefs2, form, component)
Compute inner product.
Definition: sll_m_maxwell_3d_fem_fft.F90:452
sll_m_maxwell_3d_fem_fft::multiply_mass_3d_fem_fft
subroutine multiply_mass_3d_fem_fft(self, deg, coefs_in, coefs_out)
Multiply by the mass matrix.
Definition: sll_m_maxwell_3d_fem_fft.F90:800
sll_m_maxwell_3d_fem_fft::multiply_mass_all
subroutine multiply_mass_all(self, deg, coefs_in, coefs_out)
Multiply by the mass matrix.
Definition: sll_m_maxwell_3d_fem_fft.F90:940
sll_m_maxwell_3d_fem_fft::multiply_mass_1form
subroutine multiply_mass_1form(self, coefs_in, coefs_out)
Helper function for multiply_mass.
Definition: sll_m_maxwell_3d_fem_fft.F90:827
sll_m_maxwell_3d_fem_fft::l2norm_squared_3d_fem_fft
real(kind=f64) function l2norm_squared_3d_fem_fft(self, coefs, form, component)
Compute square of the L2norm.
Definition: sll_m_maxwell_3d_fem_fft.F90:439
sll_m_maxwell_3d_fem_fft::multiply_mass_3dkron
subroutine multiply_mass_3dkron(self, mass_line_1, mass_line_2, mass_line_3, coefs_in, coefs_out)
Multiply by the mass matrix.
Definition: sll_m_maxwell_3d_fem_fft.F90:881
sll_m_poisson_3d_fem_fft
Definition: sll_m_poisson_3d_fem_fft.F90:1
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_sparse
Utilites for Maxwell solver's with spline finite elements using sparse matrix linear solvers.
Definition: sll_m_spline_fem_utilities_sparse.F90:10
sll_m_spline_fem_utilities_sparse::sll_s_spline_fem_mass1d
subroutine, public sll_s_spline_fem_mass1d(n_cells, s_deg, mass_line, matrix)
Set up 1d mass matrix (sparsity pattern and values)
Definition: sll_m_spline_fem_utilities_sparse.F90:191
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_spline_fem_multiply_massmixed
subroutine, public sll_s_spline_fem_multiply_massmixed(n_cells, degree, mass, invec, outvec)
Multiplication of the mix mass matrix given by a mass line mass.
Definition: sll_m_spline_fem_utilities.F90:644
sll_m_spline_fem_utilities::sll_s_spline_fem_multiply_mass
subroutine, public sll_s_spline_fem_multiply_mass(n_cells, degree, mass, invec, outvec)
Multiply the vector invec with the spline FEM mass matrix with mass line mass.
Definition: sll_m_spline_fem_utilities.F90:594
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
sll_m_spline_fem_utilities::sll_s_spline_fem_mixedmass_line
subroutine, public sll_s_spline_fem_mixedmass_line(deg, mass_line)
Helper function to find line of mass matrix ( N_i^p N_j^{p-1}))
Definition: sll_m_spline_fem_utilities.F90:526
sll_m_spline_fem_utilities::sll_s_spline_fem_mass_line
subroutine, public sll_s_spline_fem_mass_line(degree, mass_line)
Computes the mass line for a mass matrix with degree splines.
Definition: sll_m_spline_fem_utilities.F90:65
sll_m_spline_fem_utilities::sll_s_spline_fem_compute_mass_eig
subroutine, public sll_s_spline_fem_compute_mass_eig(n_cells, degree, mass_line, eig_values_mass)
Compute eigenvalues of mass matrix with given line mass_line.
Definition: sll_m_spline_fem_utilities.F90:572
sll_m_linear_operator_kron::sll_t_linear_operator_kron
class for a linear operator
Definition: sll_m_linear_operator_kron.F90:34
sll_m_linear_operator_maxwell_eb_schur::sll_t_linear_operator_maxwell_eb_schur
Definition: sll_m_linear_operator_maxwell_eb_schur.F90:22
sll_m_linear_solver_spline_mass_fft::sll_t_linear_solver_spline_mass_fft
Linear solver for FFT-based inversion of 3d tensor product of circulant matrices (extending the abstr...
Definition: sll_m_linear_solver_spline_mass_fft.F90:26
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_fem_fft::sll_t_maxwell_3d_fem_fft
Definition: sll_m_maxwell_3d_fem_fft.F90:72
sll_m_poisson_3d_fem_fft::sll_t_poisson_3d_fem_fft
Definition: sll_m_poisson_3d_fem_fft.F90:30
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