Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_maxwell_clamped_3d_trafo.F90
Go to the documentation of this file.
1 
8 
9 module sll_m_maxwell_clamped_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_poisson_clamped_3d, only : &
26  sll_t_linear_operator_poisson_clamped_3d
27 
28  use sll_m_linear_operator_schur_eb_cl_3d, only : &
29  sll_t_linear_operator_schur_eb_cl_3d
30 
31  use sll_m_linear_solver_cg, only : &
32  sll_t_linear_solver_cg
33 
34  use sll_m_low_level_bsplines, only: &
35  sll_s_uniform_bsplines_eval_basis
36 
37  use sll_m_mapping_3d, only: &
38  sll_t_mapping_3d
39 
40  use sll_m_matrix_csr, only: &
41  sll_t_matrix_csr
42 
43  use sll_m_maxwell_3d_base, only: &
44  sll_i_function_3d_real64, &
45  sll_c_maxwell_3d_base
46 
47  use sll_m_preconditioner_fft, only : &
48  sll_t_preconditioner_fft
49 
50  use sll_m_preconditioner_singular, only : &
51  sll_t_preconditioner_singular
52 
53  use sll_m_profile_functions, only: &
54  sll_t_profile_functions
55 
56  use sll_m_spline_fem_utilities_3d_clamped, only: &
57  sll_s_spline_fem_mass3d_clamped, &
58  sll_s_spline_fem_mixedmass3d_clamped, &
59  sll_s_multiply_g_clamped, &
60  sll_s_multiply_gt_clamped, &
61  sll_s_multiply_c_clamped, &
62  sll_s_multiply_ct_clamped
63 
64  use sll_m_spline_fem_utilities_3d, only: &
65  sll_s_spline_fem_mass3d, &
66  sll_s_spline_fem_mixedmass3d
67 
68  use sll_m_splines_pp, only: &
69  sll_t_spline_pp_1d, &
70  sll_s_spline_pp_init_1d, &
71  sll_s_spline_pp_free_1d, &
72  sll_f_spline_pp_horner_1d
73 
74  use sll_m_preconditioner_poisson_fft, only : &
75  sll_t_preconditioner_poisson_fft
76 
77  implicit none
78 
79  public :: &
80  sll_t_maxwell_clamped_3d_trafo
81 
82  private
83  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
84 
85  type, extends(sll_c_maxwell_3d_base) :: sll_t_maxwell_clamped_3d_trafo
86 
87  type(sll_t_spline_pp_1d), pointer :: spline0_pp
88  type(sll_t_spline_pp_1d), pointer :: spline1_pp
89 
90  sll_real64, allocatable :: work0(:)
91  sll_real64, allocatable :: work01(:)
92  sll_real64, allocatable :: work02(:)
93  sll_real64, allocatable :: work1(:)
94  sll_real64, allocatable :: work12(:)
95  sll_real64, allocatable :: work2(:)
96  sll_real64, allocatable :: work22(:)
97  type(sll_t_matrix_csr) :: mass0
98  type(sll_t_matrix_csr) :: mass1(3,3)
99  type(sll_t_matrix_csr) :: mass2(3,3)
100  type(sll_t_linear_operator_block) :: mass1_operator
101  type(sll_t_linear_operator_block) :: mass2_operator
102  type(sll_t_linear_solver_cg) :: mass0_solver
103  type(sll_t_linear_solver_cg) :: mass1_solver
104  type(sll_t_linear_solver_cg) :: mass2_solver
105  type(sll_t_linear_operator_poisson_clamped_3d) :: poisson_matrix
106  type(sll_t_linear_solver_cg) :: poisson_solver
107  type( sll_t_linear_operator_schur_eb_cl_3d ) :: linear_op_schur_eb
108  type( sll_t_linear_solver_cg ) :: linear_solver_schur_eb
109  type(sll_t_mapping_3d), pointer :: map
110  type(sll_t_preconditioner_fft) :: preconditioner_fft
111  type(sll_t_preconditioner_singular) :: preconditioner1
112  type(sll_t_preconditioner_singular) :: preconditioner2
113 
114  type(sll_t_preconditioner_poisson_fft ) :: poisson_preconditioner
115  logical :: adiabatic_electrons = .false.
116 
117  contains
118 
119  procedure :: &
120  compute_e_from_b => sll_s_compute_e_from_b_3d_trafo
121  procedure :: &
122  compute_b_from_e => sll_s_compute_b_from_e_3d_trafo
123  procedure :: &
124  compute_curl_part => sll_s_compute_curl_part_3d_trafo
125  procedure :: &
126  compute_e_from_rho => sll_s_compute_e_from_rho_3d_trafo
127  procedure :: &
128  compute_rho_from_e => sll_s_compute_rho_from_e_3d_trafo
129  procedure :: &
130  compute_e_from_j => sll_s_compute_e_from_j_3d_trafo
131  procedure :: &
132  compute_phi_from_rho => sll_s_compute_phi_from_rho_3d_trafo
133  procedure :: &
134  compute_phi_from_j => sll_s_compute_phi_from_j_3d_trafo
135  procedure :: &
136  compute_rhs_from_function => sll_s_compute_rhs_trafo
137  procedure :: &
138  l2projection => l2projection_3d_trafo
139  procedure :: &
140  l2norm_squared => l2norm_squared_3d_trafo
141  procedure :: &
142  inner_product => inner_product_3d_trafo
143  procedure :: &
144  init => init_3d_trafo
145  procedure :: &
146  init_from_file => init_from_file_3d_trafo
147  procedure :: &
148  free => free_3d_trafo
149  procedure :: &
150  multiply_g
151  procedure :: &
152  multiply_gt
153  procedure :: &
154  multiply_c
155  procedure :: &
156  multiply_ct
157  procedure :: &
158  multiply_mass => multiply_mass_3d_trafo
159  procedure :: &
160  multiply_mass_inverse => multiply_mass_inverse_3d_trafo
161 
162  end type sll_t_maxwell_clamped_3d_trafo
163 
164 contains
165 
166 
168  subroutine sll_s_compute_e_from_b_3d_trafo( self, delta_t, field_in, field_out )
169  class(sll_t_maxwell_clamped_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 self%multiply_mass( [2], field_in, self%work2 )
175  call self%multiply_ct( self%work2, self%work1 )
176 
177 
178  call self%multiply_mass_inverse( 1, self%work1, self%work12)
179  ! Update e from self value
180  field_out = field_out + delta_t*self%work12
181 
182 
183  end subroutine sll_s_compute_e_from_b_3d_trafo
184 
185 
187  subroutine sll_s_compute_b_from_e_3d_trafo( self, delta_t, field_in, field_out )
188  class(sll_t_maxwell_clamped_3d_trafo) :: self
189  sll_real64, intent( in ) :: delta_t
190  sll_real64, intent( in ) :: field_in(:)
191  sll_real64, intent( inout ) :: field_out(:)
192 
193  call self%multiply_c( field_in, self%work2 )
194  ! Update b from self value
195  field_out = field_out - delta_t * self%work2
196 
197  end subroutine sll_s_compute_b_from_e_3d_trafo
198 
199 
201  subroutine sll_s_compute_curl_part_3d_trafo( self, delta_t, efield, bfield, betar )
202  class(sll_t_maxwell_clamped_3d_trafo) :: self
203  sll_real64, intent(in) :: delta_t
204  sll_real64, intent(inout) :: efield(:)
205  sll_real64, intent(inout) :: bfield(:)
206  sll_real64, optional :: betar
207  !local variables
208  sll_real64 :: factor
209 
210  if( present(betar) ) then
211  factor = betar
212  else
213  factor = 1._f64
214  end if
215 
216  self%work0 = 0._f64
217  self%work1 = 0._f64
218  self%work12 = 0._f64
219  self%work2 = 0._f64
220 
221  ! Compute C^T M2 b
222  call self%multiply_mass( [2], bfield, self%work2 )
223  call self%multiply_ct( self%work2, self%work1 )
224 
225  self%linear_op_schur_eb%sign = -delta_t**2*0.25_f64*factor
226  call self%linear_op_schur_eb%dot( efield, self%work12 )
227  self%work12 = self%work12 + delta_t*factor*self%work1
228 
229  ! Save efield dofs from previous time step for B field update
230  self%work1 = efield
231 
232  ! Invert Schur complement matrix
233  self%linear_op_schur_eb%sign = delta_t**2*0.25_f64
234  call self%linear_solver_schur_eb%set_guess( efield )
235  call self%linear_solver_schur_eb%solve( self%work12, efield)
236 
237  ! Update B field
238  self%work1 = self%work1 + efield
239  call self%compute_B_from_E( delta_t*0.5_f64, self%work1, bfield)
240 
241  end subroutine sll_s_compute_curl_part_3d_trafo
242 
243 
245  subroutine sll_s_compute_e_from_rho_3d_trafo( self, field_in, field_out )
246  class(sll_t_maxwell_clamped_3d_trafo) :: self
247  sll_real64, intent( in ) :: field_in(:)
248  sll_real64, intent( out ) :: field_out(:)
249 
250  call self%poisson_solver%solve( field_in, self%work0 )
251 
252  call self%multiply_g( self%work0, field_out )
253  field_out = -field_out
254 
255 
256  end subroutine sll_s_compute_e_from_rho_3d_trafo
257 
258 
260  subroutine sll_s_compute_rho_from_e_3d_trafo( self, field_in, field_out )
261  class(sll_t_maxwell_clamped_3d_trafo) :: self
262  sll_real64, intent( in ) :: field_in(:)
263  sll_real64, intent( out ) :: field_out(:)
264 
265  call self%multiply_mass( [1], field_in, self%work1 )
266  call self%multiply_gt( self%work1, field_out )
267  field_out = - field_out
268 
269  end subroutine sll_s_compute_rho_from_e_3d_trafo
270 
271 
273  subroutine sll_s_compute_e_from_j_3d_trafo( self, current, E, component )
274  class(sll_t_maxwell_clamped_3d_trafo) :: self
275  sll_real64, intent( in ) :: current(:)
276  sll_real64, intent( inout ) :: e(:)
277  sll_int32, intent( in ), optional :: component
278 
279  if( present(component) )then
280  print*, 'compute_e_from_j_3d_trafo not implemented for single component'
281  else
282  call self%multiply_mass_inverse( 1, current, self%work1 )
283  e = e - self%work1
284  end if
285 
286  end subroutine sll_s_compute_e_from_j_3d_trafo
287 
288 
290  subroutine sll_s_compute_phi_from_rho_3d_trafo( self, field_in, field_out, efield_dofs )
291  class(sll_t_maxwell_clamped_3d_trafo) :: 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_mass_inverse( 0, field_in, field_out )
297  call self%multiply_g( field_out, efield_dofs )
298  efield_dofs = -efield_dofs
299 
300  end subroutine sll_s_compute_phi_from_rho_3d_trafo
301 
302 
304  subroutine sll_s_compute_phi_from_j_3d_trafo( self, field_in, field_out, efield_dofs )
305  class(sll_t_maxwell_clamped_3d_trafo) :: self
306  sll_real64, intent( in ) :: field_in(:)
307  sll_real64, intent( inout ) :: field_out(:)
308  sll_real64, intent( out ) :: efield_dofs(:)
309 
310  call self%multiply_gt( field_in, self%work0 )
311  call self%multiply_mass_inverse( 0, self%work0, self%work01 )
312  field_out = field_out + self%work01
313 
314  call self%multiply_g( field_out, efield_dofs )
315  efield_dofs = -efield_dofs
316 
317  end subroutine sll_s_compute_phi_from_j_3d_trafo
318 
319 
322  subroutine sll_s_compute_rhs_trafo( self, form, component, coefs_dofs, func1, func2, func3 )
323  class(sll_t_maxwell_clamped_3d_trafo) :: self
324  sll_int32, intent( in ) :: form
325  sll_int32, intent( in ) :: component
326  sll_real64, intent( out ) :: coefs_dofs(:)
327  procedure(sll_i_function_3d_real64) :: func1
328  procedure(sll_i_function_3d_real64), optional :: func2
329  procedure(sll_i_function_3d_real64), optional :: func3
330  ! local variables
331  sll_int32 :: degree(3)
332 
333  ! Define the spline degree in the 3 dimensions, depending on form and component of the form
334  if ( form == 0 ) then
335  call rhs_zeroform( self, self%s_deg_0, func1, coefs_dofs )
336  elseif (form == 1 ) then
337  degree = self%s_deg_0
338  degree(component) = self%s_deg_1(component)
339  call rhs_oneform( self, degree, func1, func2, func3, component, coefs_dofs )
340  elseif( form == 2) then
341  degree = self%s_deg_1
342  degree(component) = self%s_deg_0(component)
343  call rhs_twoform( self, degree, func1, func2, func3, component, coefs_dofs )
344  else
345  print*, 'Wrong form.'
346  end if
347 
348  end subroutine sll_s_compute_rhs_trafo
349 
350 
352  subroutine rhs_zeroform( self, deg, func, coefs_dofs )
353  class(sll_t_maxwell_clamped_3d_trafo) :: self
354  sll_int32, intent( in ) :: deg(3)
355  procedure(sll_i_function_3d_real64) :: func
356  sll_real64, intent( out ) :: coefs_dofs(:)
357  ! local variables
358  sll_int32 :: i1, i2, i3,j1, j2, j3, k1, k2, k3, q(3), index1d
359  sll_real64 :: jacobian, c(3)
360  sll_real64, allocatable :: xw_gauss_d1(:,:), xw_gauss_d2(:,:), xw_gauss_d3(:,:)
361  sll_real64, allocatable :: bspl_d1(:,:), bspl_d2(:,:), bspl_d3(:,:)
362  sll_real64 :: scratch(maxval(deg+1))
363  type(sll_t_spline_pp_1d), pointer :: spline_pp
364 
365  q = deg+1
366  allocate(xw_gauss_d2(2, q(2)))
367  allocate(bspl_d2(q(2), deg(2)+1))
368  xw_gauss_d2 = sll_f_gauss_legendre_points_and_weights(q(2), 0.0_f64, 1.0_f64)
369  ! Compute bsplines at gauss_points
370  do k2 = 1, q(2)
371  call sll_s_uniform_bsplines_eval_basis(deg(2),xw_gauss_d2(1,k2), scratch(1:deg(2)+1))
372  bspl_d2(k2,:) = scratch(1:deg(2)+1)
373  end do
374 
375  allocate(xw_gauss_d3(2,q(3)))
376  allocate(bspl_d3(q(3), deg(3)+1 ))
377  xw_gauss_d3 = sll_f_gauss_legendre_points_and_weights(q(3), 0.0_f64, 1.0_f64)
378  ! Compute bsplines at gauss_points
379  do k3 = 1, q(3)
380  call sll_s_uniform_bsplines_eval_basis(deg(3),xw_gauss_d3(1,k3), scratch(1:deg(3)+1) )
381  bspl_d3(k3,:) = scratch(1:deg(3)+1)
382  end do
383 
384  ! rescale on [0,1] for compatibility with B-splines
385  allocate(xw_gauss_d1(2, q(1)))
386  allocate(bspl_d1(q(1), deg(1)+1 ))
387  xw_gauss_d1 = sll_f_gauss_legendre_points_and_weights(q(1), 0.0_f64, 1.0_f64)
388 
389  if( deg(1) == self%s_deg_0(1) ) then
390  spline_pp => self%spline0_pp
391  else if (deg(1) == self%s_deg_1(1))then
392  spline_pp => self%spline1_pp
393  else
394  print*, "error in compute rhs"
395  end if
396  ! Compute bsplines at gauss_points
397  bspl_d1 = 0._f64
398  do k1 = 1, q(1)
399  do j1 = 1, deg(1)+1
400  bspl_d1(k1,j1) = sll_f_spline_pp_horner_1d( deg(1), spline_pp%poly_coeffs, xw_gauss_d1(1,k1), j1)
401  end do
402  end do
403  coefs_dofs = 0._f64
404  ! Compute coefs_dofs = int f(x)N_i(x)
405  !loop over cells
406  do i3 = 1, self%n_cells(3)
407  do i2 = 1, self%n_cells(2)
408  !! divide in i1=1, deg-1, deg,n_cells-deg+1,
409  do i1 = 1, deg(1)-1
410  ! loop over support of B spline
411  do j3 = 1, deg(3)+1
412  do j2 = 1, deg(2)+1
413  do j1 = 1, deg(1)+1
414  index1d = i1+j1-1+modulo(i2-deg(2)+j2-2,self%n_cells(2))*(self%n_cells(1)+deg(1))+modulo(i3-deg(3)+j3-2,self%n_cells(3))*(self%n_cells(1)+deg(1))*self%n_cells(2)
415  ! loop over Gauss points
416  do k3 = 1, q(3)
417  c(3) = self%delta_x(3)*(xw_gauss_d3(1,k3) + real(i3-1,f64))
418  do k2 = 1, q(2)
419  c(2) = self%delta_x(2)*(xw_gauss_d2(1,k2) + real(i2-1,f64))
420  do k1 = 1, q(1)
421  c(1) = self%delta_x(1)*(xw_gauss_d1(1,k1)+ real(i1 - 1,f64))
422  jacobian = self%map%jacobian( c )
423  coefs_dofs(index1d) = coefs_dofs(index1d) + &
424  self%volume * xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* &
425  xw_gauss_d3(2,k3) *&
426  func(self%map%get_x(c) ) * &
427  sll_f_spline_pp_horner_1d( deg(1), spline_pp%poly_coeffs_boundary_left(:,:,i1), xw_gauss_d1(1,k1), j1)*&
428  bspl_d2(k2,j2)*&
429  bspl_d3(k3,j3) * abs(jacobian)
430 
431  enddo
432  enddo
433  end do
434  end do
435  end do
436  end do
437  enddo
438  do i1 = max(deg(1), 1), min(self%n_cells(1)+1-deg(1), self%n_cells(1))
439  ! loop over support of B spline
440  do j3 = 1, deg(3)+1
441  do j2 = 1, deg(2)+1
442  do j1 = 1, deg(1)+1
443  index1d = i1+j1-1+modulo(i2-deg(2)+j2-2,self%n_cells(2))*(self%n_cells(1)+deg(1))+modulo(i3-deg(3)+j3-2,self%n_cells(3))*(self%n_cells(1)+deg(1))*self%n_cells(2)
444  ! loop over Gauss points
445  do k3 = 1, q(3)
446  c(3) = self%delta_x(3)*(xw_gauss_d3(1,k3) + real(i3-1,f64))
447  do k2 = 1, q(2)
448  c(2) = self%delta_x(2)*(xw_gauss_d2(1,k2) + real(i2-1,f64))
449  do k1 = 1, q(1)
450  c(1) = self%delta_x(1)*(xw_gauss_d1(1,k1)+ real(i1 - 1,f64))
451  jacobian = self%map%jacobian( c )
452  coefs_dofs(index1d) = coefs_dofs(index1d) + &
453  self%volume * xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* &
454  xw_gauss_d3(2,k3) *&
455  func(self%map%get_x(c) ) * &
456  bspl_d1(k1,j1)*&
457  bspl_d2(k2,j2)*&
458  bspl_d3(k3,j3) * abs(jacobian)
459 
460  enddo
461  enddo
462  end do
463  end do
464  end do
465  end do
466  enddo
467  do i1 = self%n_cells(1)-deg(1)+2, self%n_cells(1)
468  ! loop over support of B spline
469  do j3 = 1, deg(3)+1
470  do j2 = 1, deg(2)+1
471  do j1 = 1, deg(1)+1
472  index1d = i1+j1-1+modulo(i2-deg(2)+j2-2,self%n_cells(2))*(self%n_cells(1)+deg(1))+modulo(i3-deg(3)+j3-2,self%n_cells(3))*(self%n_cells(1)+deg(1))*self%n_cells(2)
473  ! loop over Gauss points
474  do k3 = 1, q(3)
475  c(3) = self%delta_x(3)*(xw_gauss_d3(1,k3) + real(i3-1,f64))
476  do k2 = 1, q(2)
477  c(2) = self%delta_x(2)*(xw_gauss_d2(1,k2) + real(i2-1,f64))
478  do k1 = 1, q(1)
479  c(1) = self%delta_x(1)*(xw_gauss_d1(1,k1)+ real(i1 - 1,f64))
480  jacobian = self%map%jacobian( c )
481  coefs_dofs(index1d) = coefs_dofs(index1d) + &
482  self%volume * xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* &
483  xw_gauss_d3(2,k3) *&
484  func( self%map%get_x(c) ) * &
485  sll_f_spline_pp_horner_1d( deg(1), spline_pp%poly_coeffs_boundary_right(:,:,i1-self%n_cells(1)+deg(1)-1), xw_gauss_d1(1,k1), j1)*&
486  bspl_d2(k2,j2)*&
487  bspl_d3(k3,j3) * abs(jacobian)
488 
489  enddo
490  enddo
491  end do
492  end do
493  end do
494  end do
495  enddo
496  end do
497  end do
498 
499 
500  end subroutine rhs_zeroform
501 
502 
503 
505  subroutine rhs_oneform( self, deg, func1, func2, func3, component, coefs_dofs )
506  class(sll_t_maxwell_clamped_3d_trafo) :: self
507  sll_int32, intent( in ) :: deg(3)
508  procedure(sll_i_function_3d_real64) :: func1
509  procedure(sll_i_function_3d_real64) :: func2
510  procedure(sll_i_function_3d_real64) :: func3
511  sll_int32, intent( in ) :: component
512  sll_real64, intent( out ) :: coefs_dofs(:)
513  ! local variables
514  sll_int32 :: i1, i2, i3,j1, j2, j3, k1, k2, k3, q(3),index1d
515  sll_real64 :: jacobian, jmatrix(3,3), c(3)
516  sll_real64, allocatable :: xw_gauss_d1(:,:), xw_gauss_d2(:,:), xw_gauss_d3(:,:)
517  sll_real64, allocatable :: bspl_d1(:,:), bspl_d2(:,:), bspl_d3(:,:)
518  sll_real64 :: scratch(maxval(deg+1))
519  type(sll_t_spline_pp_1d), pointer :: spline_pp
520 
521  q = deg+1
522  allocate(xw_gauss_d2(2, q(2)))
523  allocate(bspl_d2(q(2), deg(2)+1))
524  xw_gauss_d2 = sll_f_gauss_legendre_points_and_weights(q(2), 0.0_f64, 1.0_f64)
525  ! Compute bsplines at gauss_points
526  do k2 = 1, q(2)
527  call sll_s_uniform_bsplines_eval_basis(deg(2),xw_gauss_d2(1,k2), scratch(1:deg(2)+1))
528  bspl_d2(k2,:) = scratch(1:deg(2)+1)
529  end do
530 
531  allocate(xw_gauss_d3(2,q(3)))
532  allocate(bspl_d3(q(3), deg(3)+1 ))
533  xw_gauss_d3 = sll_f_gauss_legendre_points_and_weights(q(3), 0.0_f64, 1.0_f64)
534  ! Compute bsplines at gauss_points
535  do k3 = 1, q(3)
536  call sll_s_uniform_bsplines_eval_basis(deg(3),xw_gauss_d3(1,k3), scratch(1:deg(3)+1) )
537  bspl_d3(k3,:) = scratch(1:deg(3)+1)
538  end do
539 
540  ! rescale on [0,1] for compatibility with B-splines
541  allocate(xw_gauss_d1(2,q(1)))
542  allocate(bspl_d1(q(1), deg(1)+1 ))
543  xw_gauss_d1 = sll_f_gauss_legendre_points_and_weights(q(1), 0.0_f64, 1.0_f64)
544 
545  if( deg(1) == self%s_deg_0(1) ) then
546  spline_pp => self%spline0_pp
547  else if (deg(1) == self%s_deg_1(1))then
548  spline_pp => self%spline1_pp
549  else
550  print*, "error in compute rhs"
551  end if
552  ! Compute bsplines at gauss_points
553  bspl_d1 = 0._f64
554  do k1 = 1, q(1)
555  do j1 = 1, deg(1)+1
556  bspl_d1(k1,j1) = sll_f_spline_pp_horner_1d( deg(1), spline_pp%poly_coeffs, xw_gauss_d1(1,k1), j1)
557  end do
558  end do
559  coefs_dofs = 0._f64
560  ! Compute coefs_dofs = int f(x)N_i(x)
561  !loop over cells
562  do i3 = 1, self%n_cells(3)
563  do i2 = 1, self%n_cells(2)
564  !! divide in i1=1, deg-1, deg,n_cells-deg+1,
565  do i1 = 1, deg(1)-1
566  ! loop over support of B spline
567  do j3 = 1, deg(3)+1
568  do j2 = 1, deg(2)+1
569  do j1 = 1, deg(1)+1
570  index1d = i1+j1-1+modulo(i2-deg(2)+j2-2,self%n_cells(2))*(self%n_cells(1)+deg(1))+modulo(i3-deg(3)+j3-2,self%n_cells(3))*(self%n_cells(1)+deg(1))*self%n_cells(2)
571  ! loop over Gauss points
572  do k3 = 1, q(3)
573  c(3) = self%delta_x(3)*(xw_gauss_d3(1,k3) + real(i3-1,f64))
574  do k2 = 1, q(2)
575  c(2) = self%delta_x(2)*(xw_gauss_d2(1,k2) + real(i2-1,f64))
576  do k1 = 1, q(1)
577  c(1) = self%delta_x(1)*(xw_gauss_d1(1,k1)+ real(i1-1,f64))
578  jacobian = self%map%jacobian( c )
579  jmatrix = self%map%jacobian_matrix_inverse_transposed( c )
580  coefs_dofs(index1d) = coefs_dofs(index1d) + &
581  self%volume * xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* &
582  xw_gauss_d3(2,k3) *&
583  ( func1( self%map%get_x(c) )*jmatrix(1,component) + &
584  func2( self%map%get_x(c) )*jmatrix(2,component) + &
585  func3( self%map%get_x(c) )*jmatrix(3,component)) * &
586  sll_f_spline_pp_horner_1d( deg(1), spline_pp%poly_coeffs_boundary_left(:,:,i1), xw_gauss_d1(1,k1), j1)*&
587  bspl_d2(k2,j2)*&
588  bspl_d3(k3,j3)* abs(jacobian)
589  enddo
590  enddo
591  end do
592  end do
593  end do
594  end do
595  enddo
596  do i1 = max(deg(1), 1), min(self%n_cells(1)+1-deg(1), self%n_cells(1))
597  ! loop over support of B spline
598  do j3 = 1, deg(3)+1
599  do j2 = 1, deg(2)+1
600  do j1 = 1, deg(1)+1
601  index1d = i1+j1-1+modulo(i2-deg(2)+j2-2,self%n_cells(2))*(self%n_cells(1)+deg(1))+modulo(i3-deg(3)+j3-2,self%n_cells(3))*(self%n_cells(1)+deg(1))*self%n_cells(2)
602  ! loop over Gauss points
603  do k3 = 1, q(3)
604  c(3) = self%delta_x(3)*(xw_gauss_d3(1,k3) + real(i3-1,f64))
605  do k2 = 1, q(2)
606  c(2) = self%delta_x(2)*(xw_gauss_d2(1,k2) + real(i2-1,f64))
607  do k1 = 1, q(1)
608  c(1) = self%delta_x(1)*(xw_gauss_d1(1,k1)+ real(i1-1,f64))
609  jacobian = self%map%jacobian( c )
610  jmatrix = self%map%jacobian_matrix_inverse_transposed( c )
611  coefs_dofs(index1d) = coefs_dofs(index1d) + &
612  self%volume * xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* &
613  xw_gauss_d3(2,k3) *&
614  ( func1( self%map%get_x(c) )*jmatrix(1,component) + &
615  func2( self%map%get_x(c) )*jmatrix(2,component) + &
616  func3( self%map%get_x(c) )*jmatrix(3,component)) * &
617  bspl_d1(k1,j1)*&
618  bspl_d2(k2,j2)*&
619  bspl_d3(k3,j3) * abs(jacobian)
620  enddo
621  enddo
622  end do
623  end do
624  end do
625  end do
626  enddo
627  do i1 = self%n_cells(1)-deg(1)+2, self%n_cells(1)
628  ! loop over support of B spline
629  do j3 = 1, deg(3)+1
630  do j2 = 1, deg(2)+1
631  do j1 = 1, deg(1)+1
632  index1d = i1+j1-1+modulo(i2-deg(2)+j2-2,self%n_cells(2))*(self%n_cells(1)+deg(1))+modulo(i3-deg(3)+j3-2,self%n_cells(3))*(self%n_cells(1)+deg(1))*self%n_cells(2)
633  ! loop over Gauss points
634  do k3 = 1, q(3)
635  c(3) = self%delta_x(3)*(xw_gauss_d3(1,k3) + real(i3-1,f64))
636  do k2 = 1, q(2)
637  c(2) = self%delta_x(2)*(xw_gauss_d2(1,k2) + real(i2-1,f64))
638  do k1 = 1, q(1)
639  c(1) = self%delta_x(1)*(xw_gauss_d1(1,k1)+ real(i1-1,f64))
640  jacobian = self%map%jacobian( c )
641  jmatrix = self%map%jacobian_matrix_inverse_transposed( c )
642  coefs_dofs(index1d) = coefs_dofs(index1d) + &
643  self%volume * xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* &
644  xw_gauss_d3(2,k3) *&
645  ( func1( self%map%get_x(c) )*jmatrix(1,component) + &
646  func2( self%map%get_x(c) )*jmatrix(2,component) + &
647  func3( self%map%get_x(c) )*jmatrix(3,component)) * &
648  sll_f_spline_pp_horner_1d( deg(1), spline_pp%poly_coeffs_boundary_right(:,:,i1-self%n_cells(1)+deg(1)-1), xw_gauss_d1(1,k1), j1)*&
649  bspl_d2(k2,j2)*&
650  bspl_d3(k3,j3) *abs(jacobian)
651  enddo
652  enddo
653  end do
654  end do
655  end do
656  end do
657  enddo
658  end do
659  end do
660 
661 
662  end subroutine rhs_oneform
663 
664 
665 
667  subroutine rhs_twoform( self, deg, func1, func2, func3, component, coefs_dofs )
668  class(sll_t_maxwell_clamped_3d_trafo) :: self
669  sll_int32, intent( in ) :: deg(3)
670  procedure(sll_i_function_3d_real64) :: func1
671  procedure(sll_i_function_3d_real64) :: func2
672  procedure(sll_i_function_3d_real64) :: func3
673  sll_int32, intent( in ) :: component
674  sll_real64, intent( out ) :: coefs_dofs(:)
675  ! local variables
676  sll_int32 :: i1, i2, i3,j1, j2, j3, k1, k2, k3, q(3), index1d
677  sll_real64 :: jmatrix(3,3), c(3)
678  sll_real64, allocatable :: xw_gauss_d1(:,:), xw_gauss_d2(:,:), xw_gauss_d3(:,:)
679  sll_real64, allocatable :: bspl_d1(:,:), bspl_d2(:,:), bspl_d3(:,:)
680  sll_real64 :: scratch(maxval(deg+1))
681  type(sll_t_spline_pp_1d), pointer :: spline_pp
682 
683  q = deg+1
684  allocate(xw_gauss_d2(2, q(2)))
685  allocate(bspl_d2(q(2), deg(2)+1))
686  xw_gauss_d2 = sll_f_gauss_legendre_points_and_weights(q(2), 0.0_f64, 1.0_f64)
687  ! Compute bsplines at gauss_points
688  do k2 = 1, q(2)
689  call sll_s_uniform_bsplines_eval_basis(deg(2),xw_gauss_d2(1,k2), scratch(1:deg(2)+1))
690  bspl_d2(k2,:) = scratch(1:deg(2)+1)
691  end do
692 
693  allocate(xw_gauss_d3(2,q(3)))
694  allocate(bspl_d3(q(3), deg(3)+1 ))
695  xw_gauss_d3 = sll_f_gauss_legendre_points_and_weights(q(3), 0.0_f64, 1.0_f64)
696  ! Compute bsplines at gauss_points
697  do k3 = 1, q(3)
698  call sll_s_uniform_bsplines_eval_basis(deg(3),xw_gauss_d3(1,k3), scratch(1:deg(3)+1) )
699  bspl_d3(k3,:) = scratch(1:deg(3)+1)
700  end do
701 
702  ! rescale on [0,1] for compatibility with B-splines
703  allocate(xw_gauss_d1(2,q(1)))
704  allocate(bspl_d1(q(1), deg(1)+1 ))
705  xw_gauss_d1 = sll_f_gauss_legendre_points_and_weights(q(1), 0.0_f64, 1.0_f64)
706 
707  if( deg(1) == self%s_deg_0(1) ) then
708  spline_pp => self%spline0_pp
709  else if (deg(1) == self%s_deg_1(1))then
710  spline_pp => self%spline1_pp
711  else
712  print*, "error in compute rhs"
713  end if
714  ! Compute bsplines at gauss_points
715  bspl_d1 = 0._f64
716  do k1 = 1, q(1)
717  do j1 = 1, deg(1)+1
718  bspl_d1(k1,j1) = sll_f_spline_pp_horner_1d( deg(1), spline_pp%poly_coeffs, xw_gauss_d1(1,k1), j1)
719  end do
720  end do
721  coefs_dofs = 0._f64
722  ! Compute coefs_dofs = int f(x)N_i(x)
723  !loop over cells
724  do i3 = 1, self%n_cells(3)
725  do i2 = 1, self%n_cells(2)
726  !! divide in i1=1, deg-1, deg,n_cells-deg+1,
727  do i1 = 1, deg(1)-1
728  ! loop over support of B spline
729  do j3 = 1, deg(3)+1
730  do j2 = 1, deg(2)+1
731  do j1 = 1, deg(1)+1
732  index1d=i1+j1-1+modulo(i2-deg(2)+j2-2,self%n_cells(2))*(self%n_cells(1)+deg(1))+modulo(i3-deg(3)+j3-2,self%n_cells(3))*(self%n_cells(1)+deg(1))*self%n_cells(2)
733  ! loop over Gauss points
734  do k3=1, q(3)
735  c(3) = self%delta_x(3)*(xw_gauss_d3(1,k3) + real(i3-1,f64))
736  do k2=1, q(2)
737  c(2) = self%delta_x(2)*(xw_gauss_d2(1,k2) + real(i2-1,f64))
738  do k1=1, q(1)
739  c(1) = self%delta_x(1)*(xw_gauss_d1(1,k1)+ real(i1 - 1,f64))
740  jmatrix=self%map%jacobian_matrix( c )!* sign( 1._f64, self%map%jacobian( c ) )
741  coefs_dofs(index1d) = coefs_dofs(index1d) + &
742  self%volume * xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* &
743  xw_gauss_d3(2,k3) *&
744  ( func1( self%map%get_x(c) )*jmatrix(1,component) + &
745  func2( self%map%get_x(c) )*jmatrix(2,component) + &
746  func3( self%map%get_x(c) )*jmatrix(3,component)) * &
747  sll_f_spline_pp_horner_1d( deg(1), spline_pp%poly_coeffs_boundary_left(:,:,i1), xw_gauss_d1(1,k1), j1)*&
748  bspl_d2(k2,j2)*&
749  bspl_d3(k3,j3)
750 
751  enddo
752  enddo
753  end do
754  end do
755  end do
756  end do
757  enddo
758  do i1 = max(deg(1), 1), min(self%n_cells(1)+1-deg(1), self%n_cells(1))
759  ! loop over support of B spline
760  do j3 = 1, deg(3)+1
761  do j2 = 1, deg(2)+1
762  do j1 = 1, deg(1)+1
763  index1d=i1+j1-1+modulo(i2-deg(2)+j2-2,self%n_cells(2))*(self%n_cells(1)+deg(1))+modulo(i3-deg(3)+j3-2,self%n_cells(3))*(self%n_cells(1)+deg(1))*self%n_cells(2)
764  ! loop over Gauss points
765  do k3=1, q(3)
766  c(3) = self%delta_x(3)*(xw_gauss_d3(1,k3) + real(i3-1,f64))
767  do k2=1, q(2)
768  c(2) = self%delta_x(2)*(xw_gauss_d2(1,k2) + real(i2-1,f64))
769  do k1=1, q(1)
770  c(1) = self%delta_x(1)*(xw_gauss_d1(1,k1)+ real(i1 - 1,f64))
771  jmatrix=self%map%jacobian_matrix( c )!* sign( 1._f64, self%map%jacobian( c ) )
772  coefs_dofs(index1d) = coefs_dofs(index1d) + &
773  self%volume * xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* &
774  xw_gauss_d3(2,k3) *&
775  ( func1( self%map%get_x(c) )*jmatrix(1,component) + &
776  func2( self%map%get_x(c) )*jmatrix(2,component) + &
777  func3( self%map%get_x(c) )*jmatrix(3,component)) * &
778  bspl_d1(k1,j1)*&
779  bspl_d2(k2,j2)*&
780  bspl_d3(k3,j3)
781 
782  enddo
783  enddo
784  end do
785  end do
786  end do
787  end do
788  enddo
789  do i1 = self%n_cells(1)-deg(1)+2, self%n_cells(1)
790  ! loop over support of B spline
791  do j3 = 1, deg(3)+1
792  do j2 = 1, deg(2)+1
793  do j1 = 1, deg(1)+1
794  index1d=i1+j1-1+modulo(i2-deg(2)+j2-2,self%n_cells(2))*(self%n_cells(1)+deg(1))+modulo(i3-deg(3)+j3-2,self%n_cells(3))*(self%n_cells(1)+deg(1))*self%n_cells(2)
795  ! loop over Gauss points
796  do k3=1, q(3)
797  c(3) = self%delta_x(3)*(xw_gauss_d3(1,k3) + real(i3-1,f64))
798  do k2=1, q(2)
799  c(2) = self%delta_x(2)*(xw_gauss_d2(1,k2) + real(i2-1,f64))
800  do k1=1, q(1)
801  c(1) = self%delta_x(1)*(xw_gauss_d1(1,k1)+ real(i1 - 1,f64))
802  jmatrix=self%map%jacobian_matrix( c )!* sign( 1._f64, self%map%jacobian( c ) )
803  coefs_dofs(index1d) = coefs_dofs(index1d) + &
804  self%volume * xw_gauss_d1(2,k1)* xw_gauss_d2(2,k2)* &
805  xw_gauss_d3(2,k3) *&
806  ( func1( self%map%get_x(c) )*jmatrix(1,component) + &
807  func2( self%map%get_x(c) )*jmatrix(2,component) + &
808  func3( self%map%get_x(c) )*jmatrix(3,component)) * &
809  sll_f_spline_pp_horner_1d( deg(1), spline_pp%poly_coeffs_boundary_right(:,:,i1-self%n_cells(1)+deg(1)-1), xw_gauss_d1(1,k1), j1)*&
810  bspl_d2(k2,j2)*&
811  bspl_d3(k3,j3)
812  enddo
813  enddo
814  end do
815  end do
816  end do
817  end do
818  enddo
819  end do
820  end do
821 
822  end subroutine rhs_twoform
823 
824 
826  subroutine l2projection_3d_trafo( self, form, component, coefs_dofs, func1, func2, func3 )
827  class(sll_t_maxwell_clamped_3d_trafo) :: self
828  sll_int32, intent( in ) :: form
829  sll_int32, intent( in ) :: component
830  sll_real64, intent( out ) :: coefs_dofs(:)
831  procedure(sll_i_function_3d_real64) :: func1
832  procedure(sll_i_function_3d_real64), optional :: func2
833  procedure(sll_i_function_3d_real64), optional :: func3
834 
835  if( present(func2) .and. present(func3) ) then
836  select case( form )
837  case( 1 )
838  ! Compute right-hand-side
839  call sll_s_compute_rhs_trafo( self, 1, 1, self%work1(1:self%n_total1), func1, func2, func3 )
840  call sll_s_compute_rhs_trafo( self, 1, 2, self%work1(1+self%n_total1:self%n_total1+self%n_total0), func1, func2, func3 )
841  call sll_s_compute_rhs_trafo( self, 1, 3, self%work1(1+self%n_total1+self%n_total0:self%n_total1+self%n_total0*2), func1, func2, func3 )
842  call self%multiply_mass_inverse( 1, self%work1, coefs_dofs )
843  case( 2 )
844  ! Compute right-hand-side
845  call sll_s_compute_rhs_trafo( self, 2, 1, self%work2(1:self%n_total0), func1, func2, func3 )
846  call sll_s_compute_rhs_trafo( self, 2, 2, self%work2(1+self%n_total0:self%n_total0+self%n_total1), func1, func2, func3 )
847  call sll_s_compute_rhs_trafo( self, 2, 3, self%work2(1+self%n_total0+self%n_total1:self%n_total0+self%n_total1*2), func1, func2, func3 )
848  call self%multiply_mass_inverse( 2, self%work2, coefs_dofs )
849  case default
850  print*, 'L2projection for', form, '-form not implemented.'
851  end select
852  else
853  select case( form )
854  case(0)
855  call sll_s_compute_rhs_trafo( self, form, 0, self%work0, func1 )
856  call self%multiply_mass_inverse( 0, self%work0, coefs_dofs )
857  case default
858  print*,'l2 projection not for single function or component implemented'
859  end select
860  end if
861 
862  end subroutine l2projection_3d_trafo
863 
864 
866  function l2norm_squared_3d_trafo( self, coefs, form, component ) result ( r )
867  class(sll_t_maxwell_clamped_3d_trafo) :: self
868  sll_real64 :: coefs(:)
869  sll_int32 :: form
870  sll_int32 :: component
871  sll_real64 :: r
872 
873  r = inner_product_3d_trafo( self, coefs, coefs, form, component )
874 
875  end function l2norm_squared_3d_trafo
876 
877 
879  function inner_product_3d_trafo( self, coefs1, coefs2, form, component ) result ( r )
880  class(sll_t_maxwell_clamped_3d_trafo) :: self
881  sll_real64 :: coefs1(:)
882  sll_real64 :: coefs2(:)
883  sll_int32 :: form
884  sll_int32, optional :: component
885  sll_real64 :: r
886 
887  !local variables
888  sll_int32 :: istart, iend
889 
890  if ( form == 0 ) then
891  call self%multiply_mass( [0], coefs2, self%work0 )
892  r = sum(coefs1*self%work0)
893  elseif ( form == 1 ) then
894  if (present(component)) then
895  if( component == 1)then
896  istart = 1
897  iend = self%n_total1
898  else if( component == 2)then
899  istart = 1+self%n_total1
900  iend = self%n_total1+self%n_total0
901  else if( component == 3)then
902  istart = 1+self%n_total1+self%n_total0
903  iend = self%n_total1+self%n_total0*2
904  end if
905  else
906  istart=1
907  iend=self%n_total1+2*self%n_total0
908  end if
909  call self%multiply_mass( [1], coefs2, self%work1 )
910  r = sum(coefs1(istart:iend)*self%work1(istart:iend))
911  elseif( form == 2) then
912  if (present(component)) then
913  if( component == 1)then
914  istart = 1
915  iend = self%n_total0
916  else if( component == 2)then
917  istart = 1+self%n_total0
918  iend = self%n_total0+self%n_total1
919  else if( component == 3)then
920  istart = 1+self%n_total0+self%n_total1
921  iend = self%n_total0+self%n_total1*2
922  end if
923  else
924  istart=1
925  iend=self%n_total0+2*self%n_total1
926  end if
927  call self%multiply_mass( [2], coefs2, self%work2 )
928  r = sum(coefs1(istart:iend)*self%work2(istart:iend))
929  else
930  print*, 'Wrong form.'
931  end if
932 
933  end function inner_product_3d_trafo
934 
935 
937  subroutine init_3d_trafo( self, domain, n_cells, s_deg_0, boundary, map, mass_tolerance, poisson_tolerance, solver_tolerance, adiabatic_electrons, profile )
938  class(sll_t_maxwell_clamped_3d_trafo), intent( inout ) :: self
939  sll_real64, intent(in) :: domain(3,2)
940  sll_int32, intent( in ) :: n_cells(3)
941  sll_int32, intent( in ) :: s_deg_0(3)
942  sll_int32, intent( in ) :: boundary(3)
943  type(sll_t_mapping_3d), target, intent( inout ) :: map
944  sll_real64, intent(in), optional :: mass_tolerance
945  sll_real64, intent(in), optional :: poisson_tolerance
946  sll_real64, intent(in), optional :: solver_tolerance
947  logical, intent(in), optional :: adiabatic_electrons
948  type(sll_t_profile_functions), intent(in), optional :: profile
949  ! local variables
950  sll_int32 :: j, deg1(3), deg2(3)
951 
952 
953  self%Lx = map%Lx
954 
955  if (present( mass_tolerance) ) then
956  self%mass_solver_tolerance = mass_tolerance
957  else
958  self%mass_solver_tolerance = 1d-12
959  end if
960  if (present( poisson_tolerance) ) then
961  self%poisson_solver_tolerance = poisson_tolerance
962  else
963  self%poisson_solver_tolerance = 1d-12
964  end if
965  if (present( solver_tolerance) ) then
966  self%solver_tolerance = solver_tolerance
967  else
968  self%solver_tolerance = 1d-12
969  end if
970 
971  if( present( adiabatic_electrons ) ) then
972  self%adiabatic_electrons = adiabatic_electrons
973  end if
974 
975  if( present( profile ) ) then
976  self%profile = profile
977  end if
978 
979  self%n_cells = n_cells
980  self%n_dofs = n_cells
981  self%n_dofs(1) = n_cells(1)+s_deg_0(1)
982  self%n_total = product(n_cells)
983  self%n_total0 = product(self%n_dofs)
984  self%n_total1 = (self%n_dofs(1)-1)*n_cells(2)*n_cells(3)
985  self%delta_x = 1._f64 / real(n_cells,f64)
986  self%s_deg_0 = s_deg_0
987  self%s_deg_1 = s_deg_0 - 1
988  self%volume = product(self%delta_x)
989  self%map => map
990 
991  allocate( self%spline0_pp )
992  allocate( self%spline1_pp )
993  call sll_s_spline_pp_init_1d( self%spline0_pp, s_deg_0(1), self%n_cells(1), boundary(1))
994  call sll_s_spline_pp_init_1d( self%spline1_pp, s_deg_0(1)-1, self%n_cells(1), boundary(1))
995 
996  ! Allocate scratch data
997  allocate(self%work0(1:self%n_total0))
998  allocate(self%work01(1:self%n_total0))
999  allocate(self%work02(1:self%n_total0))
1000  allocate(self%work1(1:self%n_total1+2*self%n_total0))
1001  allocate(self%work12(1:self%n_total1+2*self%n_total0))
1002  allocate(self%work2(1:self%n_total0+2*self%n_total1))
1003  allocate(self%work22(1:self%n_total0+2*self%n_total1))
1004 
1005 
1006 !!!!! Assemble the sparse diagonal mass matrices
1007  !0-form
1008  deg1 = self%s_deg_0
1009  if(self%adiabatic_electrons) then
1010  call sll_s_spline_fem_mass3d_clamped( self%n_cells, deg1, -1, self%mass0, profile_m0, self%spline0_pp )
1011  else
1012  call sll_s_spline_fem_mass3d_clamped( self%n_cells, deg1, 0, self%mass0, profile_0, self%spline0_pp )
1013  end if
1014  ! first diagonal entry
1015  deg1(1) = self%s_deg_1(1)
1016  if(deg1(1) == 0 ) then
1017  call sll_s_spline_fem_mass3d( self%n_cells, deg1, 1, self%mass1(1,1), profile_1)
1018  else
1019  call sll_s_spline_fem_mass3d_clamped( self%n_cells, deg1, 1, self%mass1(1,1),profile_1, self%spline1_pp )
1020  end if
1021  deg1 = self%s_deg_1
1022  deg1(1) = self%s_deg_0(1)
1023  call sll_s_spline_fem_mass3d_clamped( self%n_cells, deg1, 1, self%mass2(1,1), profile_2, self%spline0_pp )
1024  !second and third diagonal entry
1025  do j = 2, 3
1026  deg1 = self%s_deg_0
1027  deg1(j) = self%s_deg_1(j)
1028  call sll_s_spline_fem_mass3d_clamped( self%n_cells, deg1, j, self%mass1(j,j), profile_1, self%spline0_pp )
1029 
1030  deg1 = self%s_deg_1
1031  deg1(j) = self%s_deg_0(j)
1032  if(deg1(1) == 0 ) then
1033  call sll_s_spline_fem_mass3d( self%n_cells, deg1, j, self%mass2(j,j), profile_2 )
1034  else
1035  call sll_s_spline_fem_mass3d_clamped( self%n_cells, deg1, j, self%mass2(j,j), profile_2, self%spline1_pp )
1036  end if
1037  end do
1038 !!!!assemble mixed mass for one differential form( 1 or 2 )
1039  if(self%map%flag2d )then
1040  !off diagonal entries for the upper 2x2 submatrix
1041  deg1 = self%s_deg_0
1042  deg1(1) = self%s_deg_1(1)
1043  deg2 = self%s_deg_0
1044  deg2(2) = self%s_deg_1(2)
1045  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg1, deg2, [1,2], self%mass1(1,2), profile_m1, self%spline1_pp, self%spline0_pp )
1046  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg2, deg1, [2,1], self%mass1(2,1), profile_m1, self%spline0_pp, self%spline1_pp )
1047  deg1 = self%s_deg_1
1048  deg1(1) = self%s_deg_0(1)
1049  deg2 = self%s_deg_1
1050  deg2(2) = self%s_deg_0(2)
1051  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg1, deg2, [1,2], self%mass2(1,2), profile_m2, self%spline0_pp, self%spline1_pp )
1052  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg2, deg1, [2,1], self%mass2(2,1), profile_m2, self%spline1_pp, self%spline0_pp )
1053  if(self%map%flag3d)then
1054  ! off diagonal entries for the full 3d matrix
1055  deg1=self%s_deg_0
1056  deg1(2)=self%s_deg_1(2)
1057  deg2=self%s_deg_0
1058  deg2(3) = self%s_deg_1(3)
1059  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg1, deg2, [2,3], self%mass1(2,3), profile_m1, self%spline0_pp, self%spline0_pp )
1060  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg2, deg1, [3,2], self%mass1(3,2), profile_m1, self%spline0_pp, self%spline0_pp )
1061  deg1 = self%s_deg_1
1062  deg1(2) = self%s_deg_0(2)
1063  deg2 = self%s_deg_1
1064  deg2(3) = self%s_deg_0(3)
1065  if(deg1(1) == 0 .and. deg2(1) == 0) then
1066  call sll_s_spline_fem_mixedmass3d( self%n_cells, deg1, deg2, [2,3], self%mass2(2,3), profile_m2 )
1067  call sll_s_spline_fem_mixedmass3d( self%n_cells, deg2, deg1, [3,2], self%mass2(3,2), profile_m2 )
1068  else
1069  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg1, deg2, [2,3], self%mass2(2,3), profile_m2, self%spline1_pp, self%spline1_pp )
1070  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg2, deg1, [3,2], self%mass2(3,2), profile_m2, self%spline1_pp, self%spline1_pp )
1071  end if
1072  deg1=self%s_deg_0
1073  deg1(3)=self%s_deg_1(3)
1074  deg2=self%s_deg_0
1075  deg2(1) = self%s_deg_1(1)
1076  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg1, deg2, [3,1], self%mass1(3,1), profile_m1, self%spline0_pp, self%spline1_pp )
1077  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg2, deg1, [1,3], self%mass1(1,3), profile_m1, self%spline1_pp, self%spline0_pp )
1078  deg1 = self%s_deg_1
1079  deg1(3) = self%s_deg_0(3)
1080  deg2 = self%s_deg_1
1081  deg2(1) = self%s_deg_0(1)
1082  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg1, deg2, [3,1], self%mass2(3,1), profile_m2, self%spline1_pp, self%spline0_pp )
1083  call sll_s_spline_fem_mixedmass3d_clamped( self%n_cells, deg2, deg1, [1,3], self%mass2(1,3), profile_m2, self%spline0_pp, self%spline1_pp )
1084  end if
1085  end if
1086 
1087  call self%mass1_operator%create( 3, 3 )
1088  call self%mass2_operator%create( 3, 3 )
1089 
1090  do j=1,3
1091  call self%mass1_operator%set( j, j, self%mass1(j,j) )
1092  call self%mass2_operator%set( j, j, self%mass2(j,j) )
1093  end do
1094  if(self%map%flag2d)then
1095  call self%mass1_operator%set( 1, 2, self%mass1(1,2) )
1096  call self%mass1_operator%set( 2, 1, self%mass1(2,1) )
1097  call self%mass2_operator%set( 1, 2, self%mass2(1,2) )
1098  call self%mass2_operator%set( 2, 1, self%mass2(2,1) )
1099  if(self%map%flag3d)then
1100  do j = 2,3
1101  call self%mass1_operator%set( j, modulo(j,3)+1, self%mass1(j, modulo(j,3)+1) )
1102  call self%mass1_operator%set( modulo(j,3)+1, j, self%mass1(modulo(j,3)+1, j) )
1103  call self%mass2_operator%set( j, modulo(j,3)+1, self%mass2(j, modulo(j,3)+1) )
1104  call self%mass2_operator%set( modulo(j,3)+1, j, self%mass2(modulo(j,3)+1, j) )
1105  end do
1106  end if
1107  end if
1108 
1109  call self%preconditioner_fft%init( self%Lx, n_cells, s_deg_0, .true. )
1110 
1111  self%work12 = 1._f64
1112  call self%mass1_operator%dot(self%work12, self%work1)
1113  do j = 1, self%n_dofs(2)*self%n_dofs(3)
1114  self%work1(self%n_total1+1+(j-1)*self%n_dofs(1)) = 1._f64
1115  self%work1(self%n_total1+j*self%n_dofs(1)) = 1._f64
1116  self%work1(self%n_total1+self%n_total0+1+(j-1)*self%n_dofs(1)) = 1._f64
1117  self%work1(self%n_total1+self%n_total0+j*self%n_dofs(1)) = 1._f64
1118  end do
1119  self%work1 = 1._f64/sqrt(self%work1)
1120 
1121  self%work22 = 1._f64
1122  call self%mass2_operator%dot(self%work22, self%work2)
1123  do j = 1, self%n_dofs(2)*self%n_dofs(3)
1124  self%work2(1+(j-1)*self%n_dofs(1)) = 1._f64
1125  self%work2(j*self%n_dofs(1)) = 1._f64
1126  end do
1127  self%work2 = 1._f64/sqrt(self%work2)
1128 
1129  call self%preconditioner1%create( self%preconditioner_fft%inverse_mass1_3d, self%work1, 2*self%n_total0+self%n_total1 )
1130  call self%preconditioner2%create( self%preconditioner_fft%inverse_mass2_3d, self%work2, 2*self%n_total1+self%n_total0 )
1131 
1132 
1133  self%work1 = 0._f64
1134  self%work12 = 0._f64
1135  self%work2 = 0._f64
1136  self%work22 = 0._f64
1137 
1138 
1139  call self%mass0_solver%create( self%mass0 )
1140  call self%mass1_solver%create( self%mass1_operator, self%preconditioner1)
1141  call self%mass2_solver%create( self%mass2_operator, self%preconditioner2)
1142  self%mass0_solver%atol = self%mass_solver_tolerance
1143  self%mass1_solver%atol = self%mass_solver_tolerance/maxval(self%Lx)
1144  self%mass2_solver%atol = self%mass_solver_tolerance/maxval(self%Lx)**2
1145  !self%mass0_solver%verbose = .true.
1146  !self%mass1_solver%verbose = .true.
1147  !self%mass2_solver%verbose = .true.
1148 
1149 
1150 
1151  call self%poisson_matrix%create( self%mass1_operator, self%s_deg_0, self%n_dofs, self%delta_x)
1152  ! Preconditioner for Poisson solver based on FFT inversion for periodic case
1153  do j= 1, self%n_total0
1154  self%work01 = 0._f64
1155  self%work01(j) = 1._f64
1156  call self%poisson_matrix%dot( self%work01, self%work02 )
1157  if (abs(self%work02(j))< 1d-13) then
1158  self%work0(j) = 1._f64
1159  else
1160  self%work0(j) = 1._f64/sqrt(self%work02(j))
1161  end if
1162  end do
1163  self%work0 = 1._f64/sqrt(self%work0)
1164  call self%poisson_preconditioner%create( self%n_dofs, self%s_deg_0, self%delta_x, self%work0 )
1165  call self%poisson_solver%create( self%poisson_matrix, self%poisson_preconditioner )
1166  self%poisson_solver%atol = self%poisson_solver_tolerance
1167  !self%poisson_solver%verbose = .true.
1168 
1169  ! Only for Schur complement eb solver
1170  call self%linear_op_schur_eb%create( self%mass1_operator, self%mass2_operator, self%n_dofs, self%delta_x, self%s_deg_0 )
1171  call self%linear_solver_schur_eb%create( self%linear_op_schur_eb, self%preconditioner1 )
1172  self%linear_solver_schur_eb%atol = self%solver_tolerance
1173  !self%linear_solver_schur_eb%verbose = .true.
1174  !self%linear_solver_schur_eb%n_maxiter = 2000
1175 
1176  contains
1177  function profile_m0( x, component)
1178  sll_real64 :: profile_m0
1179  sll_real64, intent(in) :: x(3)
1180  sll_int32, optional, intent(in) :: component(:)
1181 
1182  profile_m0 = self%map%jacobian( x ) * self%profile%rho_0( x(1) )/ self%profile%T_e( x(1) )
1183 
1184  end function profile_m0
1185 
1186  function profile_0(x, component)
1187  sll_real64 :: profile_0
1188  sll_real64, intent(in) :: x(3)
1189  sll_int32, optional, intent(in) :: component(:)
1190 
1191  profile_0 = self%map%jacobian( x )
1192 
1193  end function profile_0
1194 
1195  function profile_1(x, component)
1196  sll_real64 :: profile_1
1197  sll_real64, intent(in) :: x(3)
1198  sll_int32, optional, intent(in) :: component(:)
1199 
1200  profile_1 = self%map%metric_inverse_single_jacobian( x, component(1), component(1) )
1201 
1202  end function profile_1
1203 
1204  function profile_2(x, component)
1205  sll_real64 :: profile_2
1206  sll_real64, intent(in) :: x(3)
1207  sll_int32, optional, intent(in) :: component(:)
1208 
1209  profile_2 = self%map%metric_single_jacobian( x, component(1), component(1) )
1210 
1211  end function profile_2
1212 
1213  function profile_m1(x, component)
1214  sll_real64 :: profile_m1
1215  sll_real64, intent(in) :: x(3)
1216  sll_int32, optional, intent(in) :: component(:)
1217 
1218  profile_m1 = self%map%metric_inverse_single_jacobian( x, component(1), component(2) )
1219  end function profile_m1
1220 
1221  function profile_m2(x, component)
1222  sll_real64 :: profile_m2
1223  sll_real64, intent(in) :: x(3)
1224  sll_int32, optional, intent(in) :: component(:)
1225 
1226  profile_m2 = self%map%metric_single_jacobian( x, component(1), component(2) )
1227 
1228  end function profile_m2
1229 
1230  end subroutine init_3d_trafo
1231 
1232 
1234  subroutine init_from_file_3d_trafo( self, domain, n_cells, s_deg_0, boundary, map, nml_file, adiabatic_electrons, profile )
1235  class(sll_t_maxwell_clamped_3d_trafo), intent( inout ) :: self
1236  sll_real64, intent(in) :: domain(3,2)
1237  sll_int32, intent( in ) :: n_cells(3)
1238  sll_int32, intent( in ) :: s_deg_0(3)
1239  sll_int32, intent( in ) :: boundary(3)
1240  type(sll_t_mapping_3d), target, intent( inout ) :: map
1241  character(len=*), intent(in) :: nml_file
1242  logical, intent(in), optional :: adiabatic_electrons
1243  type(sll_t_profile_functions), intent(in), optional :: profile
1244  ! local variables
1245  character(len=256) :: file_prefix
1246  sll_int32 :: input_file
1247  sll_int32 :: io_stat, io_stat0, io_stat1, rank, file_id
1248  sll_real64 :: mass_tolerance
1249  sll_real64 :: poisson_tolerance
1250  sll_real64 :: maxwell_tolerance
1251 
1252  namelist /output/ file_prefix
1253  namelist /maxwell_solver/ mass_tolerance, poisson_tolerance
1254  namelist /time_solver/ maxwell_tolerance
1255 
1256  rank = sll_f_get_collective_rank(sll_v_world_collective)
1257 
1258 
1259  if( present( adiabatic_electrons) ) then
1260  self%adiabatic_electrons = adiabatic_electrons
1261  end if
1262 
1263  if( present( profile ) ) then
1264  self%profile = profile
1265  end if
1266 
1267  ! Read in solver tolerance
1268  open(newunit = input_file, file=nml_file, status='old',iostat=io_stat)
1269  if (io_stat /= 0) then
1270  if (rank == 0 ) then
1271  print*, 'sll_m_maxwell_3d_trafo: Input file does not exist. Set default tolerance.'
1272  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1273  write(file_id, *) 'mass solver tolerance:', 1d-12
1274  write(file_id, *) 'poisson solver tolerance:', 1d-12
1275  close(file_id)
1276  end if
1277  call self%init( domain, n_cells, s_deg_0, boundary, map )
1278  else
1279  read(input_file, output,iostat=io_stat0)
1280  read(input_file, maxwell_solver,iostat=io_stat)
1281  read(input_file, time_solver,iostat=io_stat1)
1282  if (io_stat /= 0 .and. io_stat1 /= 0) then
1283  if (rank == 0 ) then
1284  print*, 'sll_m_maxwell_3d_trafo: Input parameter does not exist. Set default tolerance.'
1285  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1286  write(file_id, *) 'mass solver tolerance:', 1d-12
1287  write(file_id, *) 'poisson solver tolerance:', 1d-12
1288  close(file_id)
1289  end if
1290  call self%init( domain, n_cells, s_deg_0, boundary, map )
1291  else if (io_stat /= 0 .and. io_stat1 == 0) then
1292  call self%init( domain, n_cells, s_deg_0, boundary, map, solver_tolerance=maxwell_tolerance )
1293  else if (io_stat == 0 .and. io_stat1 /= 0) then
1294  if (rank == 0 ) then
1295  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1296  write(file_id, *) 'mass solver tolerance:', mass_tolerance
1297  write(file_id, *) 'poisson solver tolerance:', poisson_tolerance
1298  close(file_id)
1299  end if
1300  call self%init( domain, n_cells, s_deg_0, boundary, map, mass_tolerance, poisson_tolerance )
1301  else
1302  if (rank == 0 ) then
1303  open(newunit=file_id, file=trim(file_prefix)//'_parameters_used.dat', position = 'append', status='old', action='write', iostat=io_stat)
1304  write(file_id, *) 'mass solver tolerance:', mass_tolerance
1305  write(file_id, *) 'poisson solver tolerance:', poisson_tolerance
1306  close(file_id)
1307  end if
1308  call self%init( domain, n_cells, s_deg_0, boundary, map, mass_tolerance, poisson_tolerance, maxwell_tolerance )
1309  end if
1310  close(input_file)
1311  end if
1312 
1313  end subroutine init_from_file_3d_trafo
1314 
1315 
1317  subroutine free_3d_trafo( self )
1318  class(sll_t_maxwell_clamped_3d_trafo) :: self
1319 
1320  self%map => null()
1321  call self%mass0_solver%free()
1322  call self%mass1_solver%free()
1323  call self%mass2_solver%free()
1324  call self%mass1_operator%free()
1325  call self%mass2_operator%free()
1326  call self%poisson_solver%free()
1327  call self%poisson_matrix%free()
1328  call self%linear_solver_schur_eb%free()
1329  call self%linear_op_schur_eb%free()
1330 
1331  call sll_s_spline_pp_free_1d( self%spline0_pp )
1332  call sll_s_spline_pp_free_1d( self%spline1_pp )
1333 
1334  deallocate( self%work0 )
1335  deallocate( self%work01 )
1336  deallocate( self%work1 )
1337  deallocate( self%work12 )
1338  deallocate( self%work2 )
1339  deallocate( self%work22 )
1340 
1341 
1342  end subroutine free_3d_trafo
1343 
1344 
1346  subroutine multiply_g( self, field_in, field_out )
1347  class(sll_t_maxwell_clamped_3d_trafo) :: self
1348  sll_real64, intent( in ) :: field_in(:)
1349  sll_real64, intent( out ) :: field_out(:)
1350 
1351  call sll_s_multiply_g_clamped(self%n_dofs, self%delta_x, self%s_deg_0, field_in, field_out)
1352 
1353  end subroutine multiply_g
1354 
1355 
1357  subroutine multiply_gt(self, field_in, field_out)
1358  class(sll_t_maxwell_clamped_3d_trafo) :: self
1359  sll_real64, intent( in ) :: field_in(:)
1360  sll_real64, intent( out ) :: field_out(:)
1361 
1362  call sll_s_multiply_gt_clamped(self%n_dofs, self%delta_x, self%s_deg_0, field_in, field_out)
1363 
1364  end subroutine multiply_gt
1365 
1366 
1368  subroutine multiply_c(self, field_in, field_out)
1369  class(sll_t_maxwell_clamped_3d_trafo) :: self
1370  sll_real64, intent( in ) :: field_in(:)
1371  sll_real64, intent( out ) :: field_out(:)
1372 
1373  call sll_s_multiply_c_clamped(self%n_dofs, self%delta_x, self%s_deg_0, field_in, field_out)
1374 
1375  end subroutine multiply_c
1376 
1377 
1379  subroutine multiply_ct(self, field_in, field_out)
1380  class(sll_t_maxwell_clamped_3d_trafo) :: self
1381  sll_real64, intent( in ) :: field_in(:)
1382  sll_real64, intent( out ) :: field_out(:)
1383 
1384  call sll_s_multiply_ct_clamped(self%n_dofs, self%delta_x, self%s_deg_0, field_in, field_out)
1385 
1386  end subroutine multiply_ct
1387 
1388 
1390  subroutine multiply_mass_3d_trafo( self, deg, coefs_in, coefs_out )
1391  class(sll_t_maxwell_clamped_3d_trafo) :: self
1392  sll_int32, intent( in ) :: deg(:)
1393  sll_real64, intent( in ) :: coefs_in(:)
1394  sll_real64, intent( out ) :: coefs_out(:)
1395 
1396  if( size(deg) == 1 )then
1397  sll_assert(deg(1)>=0 .and. deg(1)<=2)
1398  select case(deg(1))
1399  case(0)
1400  call self%mass0%dot( coefs_in, coefs_out )
1401  case(1)
1402  call self%mass1_operator%dot( coefs_in, coefs_out )
1403  case(2)
1404  call self%mass2_operator%dot( coefs_in, coefs_out )
1405  case default
1406  print*, 'multiply mass for other form not yet implemented'
1407  stop
1408  end select
1409  else if( size(deg) == 3 ) then
1410  coefs_out = 0._f64
1411  end if
1412 
1413  end subroutine multiply_mass_3d_trafo
1414 
1415 
1417  subroutine multiply_mass_inverse_3d_trafo( self, form, coefs_in, coefs_out )
1418  class(sll_t_maxwell_clamped_3d_trafo) :: self
1419  sll_int32, intent( in ) :: form
1420  sll_real64, intent( in ) :: coefs_in(:)
1421  sll_real64, intent( out ) :: coefs_out(:)
1422  !local variable
1423  sll_int32 :: j
1424 
1425 
1426  sll_assert(form>=0 .and. form<=2)
1427  select case(form)
1428  case(0)
1429  call self%mass0_solver%solve( coefs_in, coefs_out )
1430  do j = 1, self%n_dofs(2)*self%n_dofs(3)
1431  coefs_out(1+(j-1)*self%n_dofs(1)) = 0._f64
1432  coefs_out(j*self%n_dofs(1)) = 0._f64
1433  end do
1434  case(1)
1435  call self%mass1_solver%solve( coefs_in, coefs_out )
1436  do j = 1, self%n_dofs(2)*self%n_dofs(3)
1437  coefs_out(self%n_total1+1+(j-1)*self%n_dofs(1)) = 0._f64
1438  coefs_out(self%n_total1+j*self%n_dofs(1)) = 0._f64
1439  coefs_out(self%n_total1+self%n_total0+1+(j-1)*self%n_dofs(1)) = 0._f64
1440  coefs_out(self%n_total1+self%n_total0+j*self%n_dofs(1)) = 0._f64
1441  end do
1442  case(2)
1443  call self%mass2_solver%solve( coefs_in, coefs_out )
1444  do j = 1, self%n_dofs(2)*self%n_dofs(3)
1445  coefs_out(1+(j-1)*self%n_dofs(1)) = 0._f64
1446  coefs_out(j*self%n_dofs(1)) = 0._f64
1447  end do
1448  case default
1449  print*, 'multiply inverse mass for other form not yet implemented'
1450  stop
1451  end select
1452  end subroutine multiply_mass_inverse_3d_trafo
1453 
1454 
1455 
1456 end module sll_m_maxwell_clamped_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_poisson_clamped_3d
Definition: sll_m_linear_operator_poisson_clamped_3d.F90:1
sll_m_linear_operator_schur_eb_cl_3d
Definition: sll_m_linear_operator_schur_eb_cl_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_clamped_3d_trafo
Module interface to solve Maxwell's equations with coordinate transformation in 3D.
Definition: sll_m_maxwell_clamped_3d_trafo.F90:9
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:169
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:880
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:323
sll_m_maxwell_clamped_3d_trafo::multiply_mass_3d_trafo
subroutine multiply_mass_3d_trafo(self, deg, coefs_in, coefs_out)
Multiply by the mass matrix.
Definition: sll_m_maxwell_clamped_3d_trafo.F90:1391
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:305
sll_m_maxwell_clamped_3d_trafo::multiply_mass_inverse_3d_trafo
subroutine multiply_mass_inverse_3d_trafo(self, form, coefs_in, coefs_out)
Multiply by the inverse mass matrix.
Definition: sll_m_maxwell_clamped_3d_trafo.F90:1418
sll_m_maxwell_clamped_3d_trafo::multiply_gt
subroutine multiply_gt(self, field_in, field_out)
Multiply by transpose of dicrete gradient matrix.
Definition: sll_m_maxwell_clamped_3d_trafo.F90:1358
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:202
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:274
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:668
sll_m_maxwell_clamped_3d_trafo::free_3d_trafo
subroutine free_3d_trafo(self)
Finalization.
Definition: sll_m_maxwell_clamped_3d_trafo.F90:1318
sll_m_maxwell_clamped_3d_trafo::multiply_g
subroutine multiply_g(self, field_in, field_out)
Multiply by dicrete gradient matrix.
Definition: sll_m_maxwell_clamped_3d_trafo.F90:1347
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:506
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:188
sll_m_maxwell_clamped_3d_trafo::init_from_file_3d_trafo
subroutine init_from_file_3d_trafo(self, domain, n_cells, s_deg_0, boundary, map, nml_file, adiabatic_electrons, profile)
Initialization from nml file.
Definition: sll_m_maxwell_clamped_3d_trafo.F90:1235
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:291
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:261
sll_m_maxwell_clamped_3d_trafo::init_3d_trafo
subroutine init_3d_trafo(self, domain, n_cells, s_deg_0, boundary, map, mass_tolerance, poisson_tolerance, solver_tolerance, adiabatic_electrons, profile)
Initialization.
Definition: sll_m_maxwell_clamped_3d_trafo.F90:938
sll_m_maxwell_clamped_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 ( rho = G^T M_...
Definition: sll_m_maxwell_clamped_3d_trafo.F90:246
sll_m_maxwell_clamped_3d_trafo::multiply_c
subroutine multiply_c(self, field_in, field_out)
Multiply by discrete curl matrix.
Definition: sll_m_maxwell_clamped_3d_trafo.F90:1369
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:827
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:353
sll_m_maxwell_clamped_3d_trafo::multiply_ct
subroutine multiply_ct(self, field_in, field_out)
Multiply by transpose of discrete curl matrix.
Definition: sll_m_maxwell_clamped_3d_trafo.F90:1380
sll_m_maxwell_clamped_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_clamped_3d_trafo.F90:867
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_preconditioner_poisson_fft
This module is a wrapper around the spline FEM Poisson solver for the uniform grid with periodic boun...
Definition: sll_m_preconditioner_poisson_fft.F90:8
sll_m_preconditioner_singular
Module interface to solve Maxwell's equations with coordinate transformation in 3D The linear systems...
Definition: sll_m_preconditioner_singular.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_clamped
Utilites for Maxwell solver's with spline finite elements using sparse matrices.
Definition: sll_m_spline_fem_utilities_3d_clamped.F90:9
sll_m_spline_fem_utilities_3d_clamped::sll_s_multiply_ct_clamped
subroutine, public sll_s_multiply_ct_clamped(n_dofs, delta_x, s_deg_0, field_in, field_out)
Multiply by transpose of discrete curl matrix.
Definition: sll_m_spline_fem_utilities_3d_clamped.F90:1003
sll_m_spline_fem_utilities_3d_clamped::sll_s_spline_fem_mixedmass3d_clamped
subroutine, public sll_s_spline_fem_mixedmass3d_clamped(n_cells, deg1, deg2, component, matrix, profile, spline1, spline2, n_cells_min, n_cells_max)
Set up 3d clamped mixed mass matrix for specific spline degree and profile function.
Definition: sll_m_spline_fem_utilities_3d_clamped.F90:272
sll_m_spline_fem_utilities_3d_clamped::sll_s_spline_fem_mass3d_clamped
subroutine, public sll_s_spline_fem_mass3d_clamped(n_cells, deg, component, matrix, profile, spline, n_cells_min, n_cells_max)
Set up 3d clamped mass matrix for specific spline degrees and profile function.
Definition: sll_m_spline_fem_utilities_3d_clamped.F90:54
sll_m_spline_fem_utilities_3d_clamped::sll_s_multiply_g_clamped
subroutine, public sll_s_multiply_g_clamped(n_dofs, delta_x, s_deg_0, field_in, field_out)
Multiply by dicrete gradient matrix.
Definition: sll_m_spline_fem_utilities_3d_clamped.F90:539
sll_m_spline_fem_utilities_3d_clamped::sll_s_multiply_c_clamped
subroutine, public sll_s_multiply_c_clamped(n_dofs, delta_x, s_deg_0, field_in, field_out)
Multiply by discrete curl matrix.
Definition: sll_m_spline_fem_utilities_3d_clamped.F90:777
sll_m_spline_fem_utilities_3d_clamped::sll_s_multiply_gt_clamped
subroutine, public sll_s_multiply_gt_clamped(n_dofs, delta_x, s_deg_0, field_in, field_out)
Multiply by transpose of dicrete gradient matrix.
Definition: sll_m_spline_fem_utilities_3d_clamped.F90:653
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_splines_pp
Splines in pp form.
Definition: sll_m_splines_pp.F90:8
sll_m_splines_pp::sll_s_spline_pp_init_1d
subroutine, public sll_s_spline_pp_init_1d(self, degree, n_cells, boundary)
Initialize sll_t_spline_pp_1d object (Set poly_coeffs depending on degree)
Definition: sll_m_splines_pp.F90:1387
sll_m_splines_pp::sll_s_spline_pp_free_1d
subroutine, public sll_s_spline_pp_free_1d(self)
Destructor 1d.
Definition: sll_m_splines_pp.F90:1685
sll_m_splines_pp::sll_f_spline_pp_horner_1d
real(kind=f64) function, public sll_f_spline_pp_horner_1d(degree, pp_coeffs, x, index)
Perform a 1d hornerschema on the pp_coeffs at index.
Definition: sll_m_splines_pp.F90:1175
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_poisson_clamped_3d::sll_t_linear_operator_poisson_clamped_3d
Definition: sll_m_linear_operator_poisson_clamped_3d.F90:20
sll_m_linear_operator_schur_eb_cl_3d::sll_t_linear_operator_schur_eb_cl_3d
Definition: sll_m_linear_operator_schur_eb_cl_3d.F90:21
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_clamped_3d_trafo::sll_t_maxwell_clamped_3d_trafo
Definition: sll_m_maxwell_clamped_3d_trafo.F90:85
sll_m_preconditioner_fft::sll_t_preconditioner_fft
Definition: sll_m_preconditioner_fft.F90:35
sll_m_preconditioner_poisson_fft::sll_t_preconditioner_poisson_fft
Definition: sll_m_preconditioner_poisson_fft.F90:28
sll_m_preconditioner_singular::sll_t_preconditioner_singular
Definition: sll_m_preconditioner_singular.F90:30
sll_m_profile_functions::sll_t_profile_functions
Definition: sll_m_profile_functions.F90:19
sll_m_splines_pp::sll_t_spline_pp_1d
arbitrary degree 1d spline
Definition: sll_m_splines_pp.F90:93
    Report Typos and Errors
Generated on Mon Oct 23 2023 19:15:40 for Semi-Lagrangian Library by   doxygen 1.9.1