Report Typos and Errors    
Semi-Lagrangian Library
Modular library for kinetic and gyrokinetic simulations of plasmas in fusion energy devices.
sll_m_preconditioner_curl_solver_fft.F90
Go to the documentation of this file.
1 
9 
11 #include "sll_working_precision.h"
12 #include "sll_errors.h"
13 
14  use sll_m_linear_solver_abstract, only : &
16 
17  use sll_m_constants
18 
19  use sll_m_fft
20 
22 
23  implicit none
24  private
25 
27 
28 
31  sll_int32 :: n_total
32  sll_int32 :: n_dofs(3)
33  sll_real64 :: factor = 1._f64
34 
35  ! For Fourier variant
36  sll_comp64, allocatable :: eig_values_d1(:)
37  sll_comp64, allocatable :: eig_values_d1t(:)
38  sll_comp64, allocatable :: eig_values_d2(:)
39  sll_comp64, allocatable :: eig_values_d2t(:)
40  sll_comp64, allocatable :: eig_values_d3(:)
41  sll_comp64, allocatable :: eig_values_d3t(:)
42  sll_real64, allocatable :: eig_values_mass_0_1(:)
43  sll_real64, allocatable :: eig_values_mass_0_2(:)
44  sll_real64, allocatable :: eig_values_mass_0_3(:)
45  sll_real64, allocatable :: eig_values_mass_1_1(:)
46  sll_real64, allocatable :: eig_values_mass_1_2(:)
47  sll_real64, allocatable :: eig_values_mass_1_3(:)
48 
49 
50  type(sll_t_fft) :: fft1
51  type(sll_t_fft) :: fft2
52  type(sll_t_fft) :: fft3
53  type(sll_t_fft) :: ifft1
54  type(sll_t_fft) :: ifft2
55  type(sll_t_fft) :: ifft3
56 
57  contains
58 
59  procedure :: read_from_file => read_from_file_preconditioner
60  procedure :: set_verbose => set_verbose_preconditioner
61  procedure :: solve_real => solve_real_preconditioner
62  procedure :: print_info => print_info_preconditioner
63  procedure :: create => create_preconditioner
64  procedure :: free => free_preconditioner
65 
67 
68 contains
69 
70  subroutine create_preconditioner( self, n_dofs, delta_x, degree )
71  class( sll_t_preconditioner_curl_solver_fft), intent( inout ) :: self
72  sll_int32 :: n_dofs(3)
73  sll_real64 :: delta_x(3)
74  sll_int32 :: degree(3)
75  !local variables
76  sll_real64 :: angle
77  sll_int32 :: j
78  sll_comp64 :: array1d_x(n_dofs(1))
79  sll_comp64 :: array1d_y(n_dofs(2))
80  sll_comp64 :: array1d_z(n_dofs(3))
81  sll_real64, allocatable :: mass_line_0(:)
82  sll_real64, allocatable :: mass_line_1(:)
83 
84 
85  allocate( self%eig_values_mass_0_1(n_dofs(1)) )
86  allocate( self%eig_values_mass_0_2(n_dofs(2)) )
87  allocate( self%eig_values_mass_0_3(n_dofs(3)) )
88  allocate( self%eig_values_mass_1_1(n_dofs(1)) )
89  allocate( self%eig_values_mass_1_2(n_dofs(2)) )
90  allocate( self%eig_values_mass_1_3(n_dofs(3)) )
91 
92 
93  self%n_dofs = n_dofs
94  self%n_total = product(self%n_dofs)
95 
96  self%n_rows = self%n_total*3
97  self%n_cols = self%n_total*3
98 
99  self%n_global_rows = self%n_rows
100  self%n_global_cols = self%n_cols
101 
102  ! Fourier product
103 
104  ! Eigenvalues of derivative matrices
105  allocate( self%eig_values_d1(1:self%n_dofs(1)) )
106  allocate( self%eig_values_d1t(1:self%n_dofs(1) ) )
107 
108  self%eig_values_d1(1) = cmplx(0.0_f64, 0.0_f64, f64)
109  self%eig_values_d1t(1) = cmplx(0.0_f64, 0.0_f64, f64)
110  do j=2,self%n_dofs(1)
111  angle = sll_p_twopi*real(j-1,f64)/real(self%n_dofs(1), f64)
112 
113  self%eig_values_d1(j) = cmplx((1.0_f64 - cos(angle))/delta_x(1),sin(angle)/delta_x(1), f64 )
114  self%eig_values_d1t(j) = cmplx((1.0_f64 - cos(angle))/delta_x(1),-sin(angle)/delta_x(1), f64 )
115  end do
116 
117 
118  allocate( self%eig_values_d2(1:self%n_dofs(2)) )
119  allocate( self%eig_values_d2t(1:self%n_dofs(2) ) )
120 
121  self%eig_values_d2(1) = cmplx(0.0_f64, 0.0_f64, f64)
122  self%eig_values_d2t(1) = cmplx(0.0_f64, 0.0_f64, f64)
123  do j=2,self%n_dofs(2)
124  angle = sll_p_twopi*real(j-1,f64)/real(self%n_dofs(2), f64)
125 
126  self%eig_values_d2(j) = cmplx((1.0_f64 - cos(angle))/delta_x(2),sin(angle)/delta_x(2), f64 )
127  self%eig_values_d2t(j) = cmplx((1.0_f64 - cos(angle))/delta_x(2),-sin(angle)/delta_x(2), f64 )
128  end do
129 
130  allocate( self%eig_values_d3(1:self%n_dofs(3)) )
131  allocate( self%eig_values_d3t(1:self%n_dofs(3) ) )
132 
133  self%eig_values_d3(1) = cmplx(0.0_f64, 0.0_f64, f64)
134  self%eig_values_d3t(1) = cmplx(0.0_f64, 0.0_f64, f64)
135  do j=2,self%n_dofs(3)
136  angle = sll_p_twopi*real(j-1,f64)/real(self%n_dofs(3), f64)
137 
138  self%eig_values_d3(j) = cmplx((1.0_f64 - cos(angle))/delta_x(3),sin(angle)/delta_x(3), f64 )
139  self%eig_values_d3t(j) = cmplx((1.0_f64 - cos(angle))/delta_x(3),-sin(angle)/delta_x(3), f64 )
140  end do
141 
142  allocate( mass_line_0(degree(1)+1) )
143  allocate( mass_line_1(degree(1)) )
144  ! Eigenvalues of mass matrices
145  call sll_s_spline_fem_mass_line( degree(1), mass_line_0 )
146  call sll_s_spline_fem_mass_line( degree(1)-1, mass_line_1 )
147  call sll_s_spline_fem_compute_mass_eig( self%n_dofs(1), degree(1), mass_line_0*delta_x(1), self%eig_values_mass_0_1 )
148  call sll_s_spline_fem_compute_mass_eig( self%n_dofs(1), degree(1)-1, mass_line_1*delta_x(1), self%eig_values_mass_1_1 )
149  deallocate( mass_line_0 )
150  deallocate( mass_line_1 )
151 
152 
153  allocate( mass_line_0(degree(2)+1) )
154  allocate( mass_line_1(degree(2)) )
155  ! Eigenvalues of mass matrices
156  call sll_s_spline_fem_mass_line( degree(2), mass_line_0 )
157  call sll_s_spline_fem_mass_line( degree(2)-1, mass_line_1 )
158  call sll_s_spline_fem_compute_mass_eig( self%n_dofs(2), degree(2), mass_line_0*delta_x(2), self%eig_values_mass_0_2 )
159  call sll_s_spline_fem_compute_mass_eig( self%n_dofs(2), degree(2)-1, mass_line_1*delta_x(2), self%eig_values_mass_1_2 )
160  deallocate( mass_line_0 )
161  deallocate( mass_line_1 )
162 
163  allocate( mass_line_0(degree(3)+1) )
164  allocate( mass_line_1(degree(3)) )
165  ! Eigenvalues of mass matrices
166  call sll_s_spline_fem_mass_line( degree(3), mass_line_0 )
167  call sll_s_spline_fem_mass_line( degree(3)-1, mass_line_1 )
168  call sll_s_spline_fem_compute_mass_eig( self%n_dofs(3), degree(3), mass_line_0*delta_x(3), self%eig_values_mass_0_3 )
169  call sll_s_spline_fem_compute_mass_eig( self%n_dofs(3), degree(3)-1, mass_line_1*delta_x(3), self%eig_values_mass_1_3 )
170  deallocate( mass_line_0 )
171  deallocate( mass_line_1 )
172 
173  ! Initialize fft
174  call sll_s_fft_init_c2c_1d( self%fft1, self%n_dofs(1), array1d_x, array1d_x, sll_p_fft_forward)
175  call sll_s_fft_init_c2c_1d( self%fft2, self%n_dofs(2), array1d_y, array1d_y, sll_p_fft_forward)
176  call sll_s_fft_init_c2c_1d( self%fft3, self%n_dofs(3), array1d_z, array1d_z, sll_p_fft_forward)
177  call sll_s_fft_init_c2c_1d( self%ifft1, self%n_dofs(1), array1d_x, array1d_x, &
178  sll_p_fft_backward, normalized=.true.)
179  call sll_s_fft_init_c2c_1d( self%ifft2, self%n_dofs(2), array1d_y, array1d_y, &
180  sll_p_fft_backward, normalized=.true.)
181  call sll_s_fft_init_c2c_1d( self%ifft3, self%n_dofs(3), array1d_z, array1d_z, &
182  sll_p_fft_backward, normalized=.true.)
183 
184 
185  end subroutine create_preconditioner
186 
187  subroutine solve_real_preconditioner(self, rhs, unknown)
188  class( sll_t_preconditioner_curl_solver_fft), intent( inout ) :: self
189  real(kind=f64), intent(in ) :: rhs(:)
190  real(kind=f64), intent( out ) :: unknown(:)
191  !local variables
192  sll_comp64 :: scratch1(self%n_dofs(1),self%n_dofs(2), self%n_dofs(3),3)
193  sll_comp64 :: scratch2(self%n_dofs(1),self%n_dofs(2), self%n_dofs(3),3)
194  sll_comp64 :: mat(3,3), imat(3,3)
195  sll_comp64 :: array1d_x(self%n_dofs(1))
196  sll_comp64 :: array1d_y(self%n_dofs(2))
197  sll_comp64 :: array1d_z(self%n_dofs(3))
198  sll_int32 :: i,j,k
199 
200 
201 
202  ! Compute Fourier transform of input
203  call fft3d( self, array1d_x, array1d_y, array1d_z, 1, rhs(1:self%n_total), scratch1 )
204  call fft3d( self, array1d_x, array1d_y, array1d_z, 2, rhs(1+self%n_total:2*self%n_total), scratch1 )
205  call fft3d( self, array1d_x, array1d_y, array1d_z, 3, rhs(1+2*self%n_total:3*self%n_total), scratch1 )
206 
207 
208  ! Apply eigenvalues to Fourier coefficient
209  do k=1,self%n_dofs(3)
210  do j=1,self%n_dofs(2)
211  do i=1,self%n_dofs(1)
212  if ( i == 1 .or. j==1 .or. k == 1 ) then
213  scratch2(i,j,k,1) = cmplx(0.0_f64, 0.0_f64, f64)
214  scratch2(i,j,k,2) = cmplx(0.0_f64, 0.0_f64, f64)
215  scratch2(i,j,k,3) = cmplx(0.0_f64, 0.0_f64, f64)
216  else
217  mat(1,1) = (self%eig_values_d3t(k)*self%eig_values_d3(k)*&
218  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_0_2(j)*&
219  self%eig_values_mass_1_3(k),kind=f64) + &
220  self%eig_values_d2t(j)*self%eig_values_d2(j)*&
221  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_1_2(j)*&
222  self%eig_values_mass_0_3(k),kind=f64)) +&
223  cmplx(self%factor,kind=f64)* &
224  self%eig_values_d1(k)*self%eig_values_d1t(k)*&
225  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_0_2(j)*&
226  self%eig_values_mass_0_3(k),kind=f64)*&
227  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_0_2(j)*&
228  self%eig_values_mass_0_3(k),kind=f64)
229  mat(1,2) = - ( self%eig_values_d2t(j)*self%eig_values_d1(i)* &
230  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_1_2(j)*&
231  self%eig_values_mass_0_3(k),kind=f64))+&
232  cmplx(self%factor,kind=f64)* &
233  self%eig_values_d1(k)*self%eig_values_d2t(k)*&
234  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_0_2(j)*&
235  self%eig_values_mass_0_3(k),kind=f64)*&
236  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_1_2(j)*&
237  self%eig_values_mass_0_3(k),kind=f64)
238  mat(1,3) = - ( self%eig_values_d3t(k)*self%eig_values_d1(i)* &
239  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_0_2(j)*&
240  self%eig_values_mass_1_3(k),kind=f64))+&
241  cmplx(self%factor,kind=f64)* &
242  self%eig_values_d1(k)*self%eig_values_d3t(k)*&
243  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_0_2(j)*&
244  self%eig_values_mass_0_3(k),kind=f64)*&
245  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_0_2(j)*&
246  self%eig_values_mass_1_3(k),kind=f64)
247 
248  mat(2,2) = (self%eig_values_d3t(k)*self%eig_values_d3(k)*&
249  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_1_2(j)*&
250  self%eig_values_mass_1_3(k),kind=f64) + &
251  self%eig_values_d1t(i)*self%eig_values_d1(i)*&
252  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_1_2(j)*&
253  self%eig_values_mass_0_3(k),kind=f64))+&
254  cmplx(self%factor,kind=f64)* &
255  self%eig_values_d2(k)*self%eig_values_d2t(k)*&
256  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_1_2(j)*&
257  self%eig_values_mass_0_3(k),kind=f64)*&
258  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_1_2(j)*&
259  self%eig_values_mass_0_3(k),kind=f64)
260  mat(2,1) = - ( self%eig_values_d1t(i)*self%eig_values_d2(j)* &
261  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_1_2(j)*&
262  self%eig_values_mass_0_3(k),kind=f64))+&
263  cmplx(self%factor,kind=f64)* &
264  self%eig_values_d2(k)*self%eig_values_d1t(k)*&
265  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_1_2(j)*&
266  self%eig_values_mass_0_3(k),kind=f64)*&
267  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_0_2(j)*&
268  self%eig_values_mass_0_3(k),kind=f64)
269  mat(2,3) = - ( self%eig_values_d3t(k)*self%eig_values_d2(j)* &
270  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_1_2(j)*&
271  self%eig_values_mass_1_3(k),kind=f64))+&
272  cmplx(self%factor,kind=f64)* &
273  self%eig_values_d2(k)*self%eig_values_d3t(k)*&
274  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_1_2(j)*&
275  self%eig_values_mass_0_3(k),kind=f64)*&
276  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_0_2(j)*&
277  self%eig_values_mass_1_3(k),kind=f64)
278 
279  mat(3,3) = (self%eig_values_d1t(i)*self%eig_values_d1(i)*&
280  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_0_2(j)*&
281  self%eig_values_mass_1_3(k),kind=f64) + &
282  self%eig_values_d2t(j)*self%eig_values_d2(j)*&
283  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_1_2(j)*&
284  self%eig_values_mass_1_3(k),kind=f64))+&
285  cmplx(self%factor,kind=f64)* &
286  self%eig_values_d3(k)*self%eig_values_d3t(k)*&
287  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_0_2(j)*&
288  self%eig_values_mass_1_3(k),kind=f64)*&
289  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_0_2(j)*&
290  self%eig_values_mass_1_3(k),kind=f64)
291  mat(3,2) = - ( self%eig_values_d2t(j)*self%eig_values_d3(k)* &
292  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_1_2(j)*&
293  self%eig_values_mass_1_3(k),kind=f64))+&
294  cmplx(self%factor,kind=f64)* &
295  self%eig_values_d3(k)*self%eig_values_d2t(k)*&
296  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_0_2(j)*&
297  self%eig_values_mass_1_3(k),kind=f64)*&
298  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_1_2(j)*&
299  self%eig_values_mass_0_3(k),kind=f64)
300  mat(3,1) = - ( self%eig_values_d1t(i)*self%eig_values_d3(k)* &
301  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_0_2(j)*&
302  self%eig_values_mass_1_3(k),kind=f64))+&
303  cmplx(self%factor,kind=f64)* &
304  self%eig_values_d3(k)*self%eig_values_d1t(k)*&
305  cmplx(self%eig_values_mass_0_1(i)*self%eig_values_mass_0_2(j)*&
306  self%eig_values_mass_1_3(k),kind=f64)*&
307  cmplx(self%eig_values_mass_1_1(i)*self%eig_values_mass_0_2(j)*&
308  self%eig_values_mass_0_3(k),kind=f64)
309 
310  call invert3d( mat, imat )
311 
312  scratch2(i,j,k,1) = imat(1,1)*scratch1(i,j,k,1) + &
313  imat(1,2)*scratch1(i,j,k,2) + &
314  imat(1,3)*scratch1(i,j,k,3)
315  scratch2(i,j,k,2) = imat(2,1)*scratch1(i,j,k,1) + &
316  imat(2,2)*scratch1(i,j,k,2) + &
317  imat(2,3)*scratch1(i,j,k,3)
318  scratch2(i,j,k,3) = imat(3,1)*scratch1(i,j,k,1) + &
319  imat(3,2)*scratch1(i,j,k,2) + &
320  imat(3,3)*scratch1(i,j,k,3)
321  end if
322  end do
323  end do
324  end do
325  !scratch2 = scratch1
326 
327  ! Compute inverse Fourier transform of result
328  call ifft3d( self, array1d_x, array1d_y, array1d_z, 1, scratch2, unknown(1:self%n_total) )
329  call ifft3d( self, array1d_x, array1d_y, array1d_z, 2, scratch2, unknown(1+self%n_total:2*self%n_total) )
330  call ifft3d( self, array1d_x, array1d_y, array1d_z, 3, scratch2, unknown(1+2*self%n_total:3*self%n_total) )
331 
332  end subroutine solve_real_preconditioner
333 
334 
336  subroutine fft3d( self, array1d_x, array1d_y, array1d_z, inde, x, scratch1 )
337  class( sll_t_preconditioner_curl_solver_fft), intent( inout ) :: self
338  sll_comp64, intent( inout ) :: array1d_x(:)
339  sll_comp64, intent( inout ) :: array1d_y(:)
340  sll_comp64, intent( inout ) :: array1d_z(:)
341  sll_int32, intent( in ) :: inde
342  sll_real64, intent( in ) :: x(:)
343  sll_comp64, intent( out ) :: scratch1(:,:,:,:)
344  !local variables
345  sll_int32 :: ind, i,j,k
346 
347 
348  ind=0
349  do k=1,self%n_dofs(3)
350  do j=1,self%n_dofs(2)
351  do i=1,self%n_dofs(1)
352  ind = ind+1
353  array1d_x(i) = cmplx( x(ind), 0_f64, kind=f64)
354  end do
355  call sll_s_fft_exec_c2c_1d( self%fft1, array1d_x, array1d_x)
356  do i=1,self%n_dofs(1)
357  scratch1(i,j,k,inde) = array1d_x(i)
358  end do
359  end do
360  end do
361  do k=1,self%n_dofs(3)
362  do i=1,self%n_dofs(1)
363  do j=1,self%n_dofs(2)
364  array1d_y(j) = scratch1(i,j,k,inde)
365  end do
366  call sll_s_fft_exec_c2c_1d( self%fft2, array1d_y, array1d_y)
367  do j=1,self%n_dofs(2)
368  scratch1(i,j,k,inde) = array1d_y(j)
369  end do
370  end do
371  end do
372  do j=1,self%n_dofs(2)
373  do i=1,self%n_dofs(1)
374  do k=1,self%n_dofs(3)
375  array1d_z(k) = scratch1(i,j,k,inde)
376  end do
377  call sll_s_fft_exec_c2c_1d( self%fft3, array1d_z, array1d_z)
378  do k=1,self%n_dofs(3)
379  scratch1(i,j,k,inde) = array1d_z(k)
380  end do
381  end do
382  end do
383 
384  end subroutine fft3d
385 
387  subroutine ifft3d( self, array1d_x, array1d_y, array1d_z, inde, scratch, y )
388  class( sll_t_preconditioner_curl_solver_fft), intent( inout ) :: self
389  sll_comp64, intent( inout ) :: array1d_x(:)
390  sll_comp64, intent( inout ) :: array1d_y(:)
391  sll_comp64, intent( inout ) :: array1d_z(:)
392  sll_int32, intent( in ) :: inde
393  sll_real64, intent( out ) :: y(:)
394  sll_comp64, intent( inout ) :: scratch(:,:,:,:)
395  !local variables
396  sll_int32 :: ind, i,j,k
397 
398 
399  do j=1,self%n_dofs(2)
400  do i=1,self%n_dofs(1)
401  do k=1,self%n_dofs(3)
402  array1d_z(k) = scratch(i,j,k,inde)
403  end do
404  call sll_s_fft_exec_c2c_1d( self%ifft3, array1d_z, array1d_z)
405  do k=1,self%n_dofs(3)
406  scratch(i,j,k,inde) = array1d_z(k)
407  end do
408  end do
409  end do
410 
411  do k=1,self%n_dofs(3)
412  do i=1,self%n_dofs(1)
413  do j=1,self%n_dofs(2)
414  array1d_y(j) = scratch(i,j,k,inde)
415  end do
416  call sll_s_fft_exec_c2c_1d( self%ifft2, array1d_y, array1d_y)
417  do j=1,self%n_dofs(2)
418  scratch(i,j,k,inde) = array1d_y(j)
419  end do
420  end do
421  end do
422 
423  ind=0
424  do k=1,self%n_dofs(3)
425  do j=1,self%n_dofs(2)
426  do i=1,self%n_dofs(1)
427  array1d_x(i) = scratch(i,j,k,inde)
428  end do
429  call sll_s_fft_exec_c2c_1d( self%ifft1, array1d_x, array1d_x)
430 
431  do i=1,self%n_dofs(1)
432  ind = ind+1
433  y(ind) = real( array1d_x(i), kind=f64 )
434  end do
435  end do
436  end do
437 
438  end subroutine ifft3d
439 
441  subroutine invert3d( mat, mat_inv )
442  sll_comp64, intent( in ) :: mat(3,3)
443  sll_comp64, intent( out ) :: mat_inv(3,3)
444 
445  sll_comp64 :: det
446 
447  det = mat(1,1)*mat(2,2)*mat(3,3) + mat(1,2)*mat(2,3)*mat(3,1) + mat(1,3)*mat(2,1)*mat(3,2) - mat(1,3)*mat(2,2)*mat(3,1) - mat(3,3)*mat(1,2)*mat(2,1) - mat(1,1)*mat(2,3)*mat(3,2)
448 
449  mat_inv(1,1) = mat(2,2)*mat(3,3) - mat(2,3)*mat(3,2)
450  mat_inv(1,2) = mat(1,3)*mat(3,2) - mat(1,2)*mat(3,3)
451  mat_inv(1,3) = mat(1,2)*mat(2,3) - mat(1,3)*mat(2,2)
452 
453  mat_inv(2,1) = mat(2,3)*mat(3,1) - mat(2,1)*mat(3,3)
454  mat_inv(2,2) = mat(1,1)*mat(3,3) - mat(1,3)*mat(3,1)
455  mat_inv(2,3) = mat(1,3)*mat(2,1) - mat(1,1)*mat(2,3)
456 
457  mat_inv(3,1) = mat(2,1)*mat(3,2) - mat(2,2)*mat(3,1)
458  mat_inv(3,2) = mat(1,2)*mat(3,1) - mat(1,1)*mat(3,2)
459  mat_inv(3,3) = mat(1,1)*mat(2,2) - mat(1,2)*mat(2,1)
460 
461  mat_inv = mat_inv/det
462 
463  end subroutine invert3d
464 
465  subroutine read_from_file_preconditioner(self, filename)
466  class( sll_t_preconditioner_curl_solver_fft), intent( inout ) :: self
467  character(len=*), intent( in ) :: filename
468 
469  end subroutine read_from_file_preconditioner
470 
471  subroutine print_info_preconditioner(self)
472  class( sll_t_preconditioner_curl_solver_fft), intent( in ) :: self
473 
474  end subroutine print_info_preconditioner
475 
476  subroutine set_verbose_preconditioner( self, verbose )
477  class( sll_t_preconditioner_curl_solver_fft), intent( inout ) :: self
478  logical, intent( in ) :: verbose
479 
480  self%verbose = verbose
481 
482  end subroutine set_verbose_preconditioner
483 
484 
485  subroutine free_preconditioner(self)
486  class( sll_t_preconditioner_curl_solver_fft), intent( inout ) :: self
487 
488  end subroutine free_preconditioner
489 
Fortran module where set some physical and mathematical constants.
real(kind=f64), parameter, public sll_p_twopi
FFT interface for FFTW.
integer, parameter, public sll_p_fft_forward
Flag to specify forward FFT (negative sign) [value -1].
subroutine, public sll_s_fft_exec_c2c_1d(plan, array_in, array_out)
Compute fast Fourier transform in complex to complex mode.
subroutine, public sll_s_fft_init_c2c_1d(plan, nx, array_in, array_out, direction, normalized, aligned, optimization)
Create new 1d complex to complex plan.
integer, parameter, public sll_p_fft_backward
Flag to specify backward FFT (positive sign) [value 1].
module for abstract linear solver
Invert a circulant matrix based on diagonalization in Fourier space (3d version)
subroutine create_preconditioner(self, n_dofs, delta_x, degree)
subroutine ifft3d(self, array1d_x, array1d_y, array1d_z, inde, scratch, y)
Helper function.
subroutine fft3d(self, array1d_x, array1d_y, array1d_z, inde, x, scratch1)
Helper function.
subroutine invert3d(mat, mat_inv)
Helper function to invert 3x3 matrix.
Utilites for Maxwell solver's with spline finite elements.
subroutine, public sll_s_spline_fem_mass_line(degree, mass_line)
Computes the mass line for a mass matrix with degree splines.
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.
Type for FFT plan in SeLaLib.
Linear solver for FFT-based inversion of 3d tensor product of circulant matrices (extending the abstr...
    Report Typos and Errors