11 #include "sll_working_precision.h"
12 #include "sll_errors.h"
32 sll_int32 :: n_dofs(3)
33 sll_real64 :: factor = 1._f64
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(:)
72 sll_int32 :: n_dofs(3)
73 sll_real64 :: delta_x(3)
74 sll_int32 :: degree(3)
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(:)
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)) )
94 self%n_total = product(self%n_dofs)
96 self%n_rows = self%n_total*3
97 self%n_cols = self%n_total*3
99 self%n_global_rows = self%n_rows
100 self%n_global_cols = self%n_cols
105 allocate( self%eig_values_d1(1:self%n_dofs(1)) )
106 allocate( self%eig_values_d1t(1:self%n_dofs(1) ) )
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)
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 )
118 allocate( self%eig_values_d2(1:self%n_dofs(2)) )
119 allocate( self%eig_values_d2t(1:self%n_dofs(2) ) )
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)
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 )
130 allocate( self%eig_values_d3(1:self%n_dofs(3)) )
131 allocate( self%eig_values_d3t(1:self%n_dofs(3) ) )
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)
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 )
142 allocate( mass_line_0(degree(1)+1) )
143 allocate( mass_line_1(degree(1)) )
149 deallocate( mass_line_0 )
150 deallocate( mass_line_1 )
153 allocate( mass_line_0(degree(2)+1) )
154 allocate( mass_line_1(degree(2)) )
160 deallocate( mass_line_0 )
161 deallocate( mass_line_1 )
163 allocate( mass_line_0(degree(3)+1) )
164 allocate( mass_line_1(degree(3)) )
170 deallocate( mass_line_0 )
171 deallocate( mass_line_1 )
189 real(kind=f64),
intent(in ) :: rhs(:)
190 real(kind=f64),
intent( out ) :: unknown(:)
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))
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 )
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)
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)
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)
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)
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)
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) )
336 subroutine fft3d( self, array1d_x, array1d_y, array1d_z, inde, x, scratch1 )
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(:,:,:,:)
345 sll_int32 :: ind, i,j,k
349 do k=1,self%n_dofs(3)
350 do j=1,self%n_dofs(2)
351 do i=1,self%n_dofs(1)
353 array1d_x(i) = cmplx( x(ind), 0_f64, kind=f64)
356 do i=1,self%n_dofs(1)
357 scratch1(i,j,k,inde) = array1d_x(i)
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)
367 do j=1,self%n_dofs(2)
368 scratch1(i,j,k,inde) = array1d_y(j)
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)
378 do k=1,self%n_dofs(3)
379 scratch1(i,j,k,inde) = array1d_z(k)
387 subroutine ifft3d( self, array1d_x, array1d_y, array1d_z, inde, scratch, y )
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(:,:,:,:)
396 sll_int32 :: ind, i,j,k
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)
405 do k=1,self%n_dofs(3)
406 scratch(i,j,k,inde) = array1d_z(k)
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)
417 do j=1,self%n_dofs(2)
418 scratch(i,j,k,inde) = array1d_y(j)
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)
431 do i=1,self%n_dofs(1)
433 y(ind) = real( array1d_x(i), kind=f64 )
442 sll_comp64,
intent( in ) :: mat(3,3)
443 sll_comp64,
intent( out ) :: mat_inv(3,3)
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)
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)
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)
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)
461 mat_inv = mat_inv/det
467 character(len=*),
intent( in ) :: filename
478 logical,
intent( in ) :: verbose
480 self%verbose = verbose
Fortran module where set some physical and mathematical constants.
real(kind=f64), parameter, public sll_p_twopi
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 set_verbose_preconditioner(self, verbose)
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 read_from_file_preconditioner(self, filename)
subroutine print_info_preconditioner(self)
subroutine invert3d(mat, mat_inv)
Helper function to invert 3x3 matrix.
subroutine free_preconditioner(self)
subroutine solve_real_preconditioner(self, rhs, unknown)
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.
class for abstract linear solver
Linear solver for FFT-based inversion of 3d tensor product of circulant matrices (extending the abstr...