module WaveKernelClass use SparseClass use FEMDomainClass use FEMSolverClass use DEMDomainClass use SpectreAnalysisClass implicit none type :: WaveKernel_ type(CRS_) :: OmegaSqMatrix real(Real64), allocatable :: Mmatrix_diag(:) real(Real64), allocatable :: DampingRatio(:) real(real64) :: tol = dble(1.0e-25) real(real64), allocatable :: v_in1(:) real(real64), allocatable :: v_in2(:) real(real64), allocatable :: v_out1(:) real(real64), allocatable :: v_out2(:) real(real64), allocatable :: u_in1(:) real(real64), allocatable :: u_in2(:) real(real64), allocatable :: u_out1(:) real(real64), allocatable :: u_out2(:) real(real64), allocatable :: hanning_buffer(:, :) integer(int32) :: itrmax = 100 ! EbO-Scheme integer(int32) :: EbO_Algorithm = FEMDomain_Overset_GPP real(real64) :: EbO_relative_penalty = 1000.0d0 ! custom low-pass filter real(real64), allocatable, private :: weights(:) contains procedure, pass :: initWaveKernel procedure, pass :: init_by_femdomain_pointers_WK !procedure, pass :: initWaveKernel_DEM generic :: init => initWaveKernel, init_by_femdomain_pointers_WK !,initWaveKernel_DEM procedure :: setLPF => setLPF_WaveKernel procedure :: bandpass => bandpass_WaveKernel procedure :: lowpass => lowpass_WaveKernel procedure :: movingAverage => movingAverage_WaveKernel ! >>> RECCOMENDED >>> procedure, pass :: getDisplacement_and_Velocity_WaveKernel procedure, pass :: getDisplacement_and_Velocity_MPI_WaveKernel procedure, pass :: getDisplacement_and_Velocity_complex64_WaveKernel generic :: getDisplacement_and_Velocity & => getDisplacement_and_Velocity_WaveKernel, & getDisplacement_and_Velocity_complex64_WaveKernel, & getDisplacement_and_Velocity_MPI_WaveKernel ! <<< RECCOMENDED <<< procedure :: getDisplacement => getDisplacementWaveKernel procedure :: getVelocity => getVelocityWaveKernel end type contains ! ############################################################## subroutine initWaveKernel(this, FEMDomain, DOF, YoungModulus, PoissonRatio, & DampingRatio, Density) class(WaveKernel_), intent(inout) :: this integer(int32), intent(in) :: DOF type(FEMDomain_), intent(inout) :: FEMDomain real(real64), intent(in) :: YoungModulus(:), Density(:) real(real64), optional, intent(in) :: PoissonRatio(:), DampingRatio(:) type(CRS_) :: Imatrix, Mmatrix, Cmatrix !Mmatrix = FEMDomain%MassMatrix(DOF=DOF, Density=Density) call FEMDomain%setMassMatrix(DOF=DOF, Density=Density,MassMatrix=Mmatrix) this%Mmatrix_diag = Mmatrix%diag(cell_centered=.true.) call Mmatrix%remove() if (DOF == 1) then !this%OmegaSqMatrix = FEMDomain%StiffnessMatrix( & ! YoungModulus=YoungModulus) call FEMDomain%setStiffnessMatrix( & YoungModulus=YoungModulus,StiffnessMatrix=this%OmegaSqMatrix) else if (.not. present(PoissonRatio)) then print *, "[initWaveKernel] Please input PoissonRatio" stop end if !this%OmegaSqMatrix = FEMDomain%StiffnessMatrix( & ! YoungModulus=YoungModulus, PoissonRatio=PoissonRatio) call FEMDomain%setStiffnessMatrix( & YoungModulus=YoungModulus,& PoissonRatio=PoissonRatio,& StiffnessMatrix=this%OmegaSqMatrix) end if !this%OmegaSqMatrix = this%OmegaSqMatrix%divide_by(this%Mmatrix_diag) call this%OmegaSqMatrix%divide_by_vector(this%Mmatrix_diag) this%DampingRatio = zeros(this%OmegaSqMatrix%size()) if (present(DampingRatio)) then !call Imatrix%eyes(this%OmegaSqMatrix%size() ) Cmatrix = to_diag(DampingRatio*DampingRatio) ! vector to diagonal matrix this%OmegaSqMatrix = this%OmegaSqMatrix - Cmatrix this%DampingRatio = DampingRatio end if end subroutine initWaveKernel ! ############################################################## !! ############################################################## !subroutine initWaveKernel_DEM(this, DEMDomain, DOF, YoungModulus, PoissonRatio, & ! DampingRatio, Density) !class(WaveKernel_), intent(inout) :: this !integer(int32), intent(in) :: DOF !type(DEMDomain_), intent(inout) :: DEMDomain !real(real64), intent(in) :: YoungModulus(:), Density(:) !real(real64), optional, intent(in) :: PoissonRatio(:), DampingRatio(:) !type(CRS_) :: Imatrix, Mmatrix, Cmatrix ! !Mmatrix = DEMDomain%MassMatrix(DOF=DOF, Density=Density) ! !if (DOF == 1) then !this%OmegaSqMatrix = DEMDomain%StiffnessMatrix( & ! YoungModulus=YoungModulus) ! !else !if (.not. present(PoissonRatio)) then !print *, "[initWaveKernel] Please input PoissonRatio" !stop !end if !this%OmegaSqMatrix = DEMDomain%StiffnessMatrix( & ! YoungModulus=YoungModulus, PoissonRatio=PoissonRatio) ! !end if ! !this%Mmatrix_diag = Mmatrix%diag(cell_centered=.true.) !this%OmegaSqMatrix = this%OmegaSqMatrix%divide_by(this%Mmatrix_diag) ! !this%DampingRatio = zeros(this%OmegaSqMatrix%size()) ! !if (present(DampingRatio)) then !!call Imatrix%eyes(this%OmegaSqMatrix%size() ) !Cmatrix = to_diag(DampingRatio*DampingRatio) ! vector to diagonal matrix !this%OmegaSqMatrix = this%OmegaSqMatrix - Cmatrix !this%DampingRatio = DampingRatio !end if ! !end subroutine initWaveKernel ! ############################################################## ! ############################################################## subroutine init_by_femdomain_pointers_WK(this, FEMDomainPointers, DOF, YoungModulus, PoissonRatio, & DampingRatio, Density) class(WaveKernel_), intent(inout) :: this integer(int32), intent(in) :: DOF type(FEMDomainp_), intent(inout) :: FEMDomainPointers(:) type(FEMSolver_) :: solver real(real64), intent(in) :: YoungModulus(:), Density(:) real(real64), optional, intent(in) :: PoissonRatio(:), DampingRatio(:) integer(int32) :: DomainID, ElementID type(CRS_) :: Imatrix, Mmatrix, Cmatrix, Kmatrix ! initialize wave kernel function by given FEMDomain-pointers call solver%init(NumDomain=size(FEMDomainPointers)) call solver%setDomain(FEMDomainPointers=FEMDomainPointers) call solver%setCRS(DOF=DOF) !$OMP parallel do private(ElementID) do DomainID = 1, size(FEMDomainPointers) do ElementID = 1, FEMDomainPointers(DomainID)%femdomainp%ne() call solver%setMatrix(DomainID=DomainID, ElementID=ElementID, DOF=DOF, & Matrix=FEMDomainPointers(DomainID)%femdomainp%StiffnessMatrix( & ElementID=ElementID, & E=YoungModulus(get_element_idx(FEMDomainPointers, DomainID, ElementID)), & v=PoissonRatio(get_element_idx(FEMDomainPointers, DomainID, ElementID)))) end do end do !$OMP end parallel do call solver%setEbOM(penalty=this%EbO_relative_penalty, DOF=3) Kmatrix = solver%getCRS() call solver%zeros() ! mass matrix !$OMP parallel do private(ElementID) do DomainID = 1, size(FEMDomainPointers) do ElementID = 1, FEMDomainPointers(DomainID)%femdomainp%ne() call solver%setMatrix(DomainID=DomainID, ElementID=ElementID, DOF=DOF, & Matrix=FEMDomainPointers(DomainID)%femdomainp%MassMatrix( & ElementID=ElementID, & Density=Density(get_element_idx(FEMDomainPointers, DomainID, ElementID)), & DOF=3)) end do end do !$OMP end parallel do Mmatrix = solver%getCRS() this%OmegaSqMatrix = Kmatrix this%Mmatrix_diag = Mmatrix%diag(cell_centered=.true.) this%OmegaSqMatrix = this%OmegaSqMatrix%divide_by(this%Mmatrix_diag) this%DampingRatio = zeros(this%OmegaSqMatrix%size()) if (present(DampingRatio)) then !call Imatrix%eyes(this%OmegaSqMatrix%size() ) Cmatrix = to_diag(DampingRatio*DampingRatio) ! vector to diagonal matrix this%OmegaSqMatrix = this%OmegaSqMatrix - Cmatrix this%DampingRatio = DampingRatio end if end subroutine init_by_femdomain_pointers_WK ! ############################################################## ! ############################################################## subroutine bandpass_WaveKernel(this, u_n, v_n, freq_range, dt, fix_idx) class(WaveKernel_), intent(inout) :: this real(real64), intent(inout) :: u_n(:), v_n(:) real(real64), intent(in) :: freq_range(1:2) integer(int32), optional, intent(in) :: fix_idx(:) real(real64), allocatable :: u(:) real(real64), intent(in) :: dt integer(int32) :: j real(real64) :: h if (.not. allocated(this%u_in1)) then this%u_in1 = u_n this%u_in2 = u_n this%u_out1 = u_n this%u_out2 = u_n end if if (.not. allocated(this%v_in1)) then this%v_in1 = v_n this%v_in2 = v_n this%v_out1 = v_n this%v_out2 = v_n end if u_n(fix_idx) = 0.0d0 !$OMP parallel do do j = 1, size(u_n) u_n(j) = bandpass_filter(x=u_n(j), & freq_range=freq_range, & sampling_Hz=1.0d0/dt, & in1=this%u_in1(j), & in2=this%u_in2(j), & out1=this%u_out1(j), & out2=this%u_out2(j) & ) end do !$OMP end parallel do v_n(fix_idx) = 0.0d0 !$OMP parallel do do j = 1, size(v_n) v_n(j) = bandpass_filter(x=v_n(j), & freq_range=freq_range, & sampling_Hz=1.0d0/dt, & in1=this%v_in1(j), & in2=this%v_in2(j), & out1=this%v_out1(j), & out2=this%v_out2(j) & ) end do !$OMP end parallel do end subroutine ! ############################################################## ! ############################################################## subroutine lowpass_WaveKernel(this, u_n, v_n, cutoff_frequency, dt, fix_idx) class(WaveKernel_), intent(inout) :: this real(real64), intent(inout) :: u_n(:), v_n(:) real(real64), intent(in) :: cutoff_frequency integer(int32), optional, intent(in) :: fix_idx(:) real(real64), allocatable :: u(:) real(real64), intent(in) :: dt integer(int32) :: j real(real64) :: h if (.not. allocated(this%u_in1)) then this%u_in1 = u_n end if if (.not. allocated(this%v_in1)) then this%v_in1 = v_n end if u_n(fix_idx) = 0.0d0 !$OMP parallel do do j = 1, size(u_n) u_n(j) = lowpass_filter(x_n=u_n(j), & fc=cutoff_frequency, & sampling_Hz=1.0d0/dt, & buf=this%u_in1(j) & ) end do !$OMP end parallel do v_n(fix_idx) = 0.0d0 !$OMP parallel do do j = 1, size(v_n) v_n(j) = lowpass_filter(x_n=v_n(j), & fc=cutoff_frequency, & sampling_Hz=1.0d0/dt, & buf=this%v_in1(j) & ) end do !$OMP end parallel do end subroutine ! ############################################################## ! ############################################################## subroutine movingAverage_WaveKernel(this, u_n, v_n, fix_idx) class(WaveKernel_), intent(inout) :: this real(real64), intent(inout) :: u_n(:), v_n(:) integer(int32), optional, intent(in) :: fix_idx(:) real(real64), allocatable :: u(:) integer(int32) :: j real(real64) :: h if (.not. allocated(this%u_in1)) then this%u_in1 = u_n this%u_in2 = u_n this%u_out1 = u_n this%u_out2 = u_n end if if (.not. allocated(this%v_in1)) then this%v_in1 = v_n this%v_in2 = v_n this%v_out1 = v_n this%v_out2 = v_n end if u_n(fix_idx) = 0.0d0 !$OMP parallel do do j = 1, size(u_n) u_n(j) = movingAverage_filter(x=u_n(j), & in1=this%u_in1(j), & in2=this%u_in2(j) & ) end do !$OMP end parallel do v_n(fix_idx) = 0.0d0 !$OMP parallel do do j = 1, size(v_n) v_n(j) = movingAverage_filter(x=v_n(j), & in1=this%v_in1(j), & in2=this%v_in2(j) & ) end do !$OMP end parallel do end subroutine ! ############################################################## ! ############################################################## subroutine getDisplacement_and_Velocity_WaveKernel(this, u_n, v_n, dt, & fix_idx, cutoff_frequency, debug_mode, u, v, RHS) class(WaveKernel_), intent(inout) :: this real(real64), intent(in) :: u_n(:), v_n(:) real(real64), optional, intent(in) :: RHS(:) real(real64), allocatable :: du(:), dv(:) real(real64), allocatable :: u(:), v(:) integer(int32), optional, allocatable, intent(in) :: fix_idx(:) real(real64), intent(in) :: dt logical, optional, intent(in) :: debug_mode real(real64), optional, intent(in) :: cutoff_frequency integer(int32) :: j, k ! [CAUTION!!] only undamped is implemented. du = u_n if (present(fix_idx)) then if (allocated(fix_idx)) then du(fix_idx) = 0.0d0 end if end if ! U_n から U? u = u_n v = 0.0d0*v_n do k = 1, this%itrmax if (present(fix_idx)) then if (allocated(fix_idx)) then du(fix_idx) = 0.0d0 end if end if du = this%OmegaSqMatrix%matmul(du) if (present(fix_idx)) then if (allocated(fix_idx)) then du(fix_idx) = 0.0d0 end if end if if (norm(LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*du) & < this%tol) exit u = u + LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*du if (k == 1) then v = -dt*du else v = v - LPF_t_sinc_sqrt_taylor_coefficient(k=k - 1, t=dt, f_c=cutoff_frequency)*du end if end do ! V_n から V? dv = v_n if (present(fix_idx)) then if (allocated(fix_idx)) then dv(fix_idx) = 0.0d0 end if end if ! k=1 u = u + dt*dv v = v + dv if (present(fix_idx)) then if (allocated(fix_idx)) then u(fix_idx) = 0.0d0 end if end if if (present(fix_idx)) then if (allocated(fix_idx)) then v(fix_idx) = 0.0d0 end if end if do k = 1, this%itrmax if (present(fix_idx)) then if (allocated(fix_idx)) then dv(fix_idx) = 0.0d0 end if end if dv = this%OmegaSqMatrix%matmul(dv) if (present(fix_idx)) then if (allocated(fix_idx)) then dv(fix_idx) = 0.0d0 end if end if if (norm(LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*dv) & < this%tol) exit u = u + LPF_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*dv v = v + LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*dv end do if (present(RHS)) then if (allocated(this%weights)) then u = u + this%OmegaSqMatrix%tensor_wave_kernel_RHS(RHS=RHS/this%Mmatrix_diag, t=dt, tol=dble(1.0e-25), & itrmax=this%itrmax, cutoff_frequency=cutoff_frequency, weights=this%weights) v = v + LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=RHS/this%Mmatrix_diag, & DampingRatio=this%DampingRatio, itrmax=this%itrmax, tol=this%tol, weights=this%weights) else u = u + this%OmegaSqMatrix%tensor_wave_kernel_RHS(RHS=RHS/this%Mmatrix_diag, t=dt, tol=dble(1.0e-25), & itrmax=this%itrmax, cutoff_frequency=cutoff_frequency) v = v + LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=RHS/this%Mmatrix_diag, & DampingRatio=this%DampingRatio, itrmax=this%itrmax, tol=this%tol) end if end if if (present(fix_idx)) then if (allocated(fix_idx)) then u(fix_idx) = 0.0d0 end if end if ! u = gain_value*u ! v = gain_value*v end subroutine ! ############################################################## ! ############################################################## subroutine getDisplacement_and_Velocity_MPI_WaveKernel(this, u_n, v_n, dt, & fix_idx, cutoff_frequency, debug_mode, u, v, RHS, MPID, FEMDomain) class(WaveKernel_), intent(inout) :: this real(real64), intent(in) :: u_n(:), v_n(:) real(real64), optional, intent(in) :: RHS(:) real(real64), allocatable :: du(:), dv(:) real(real64), allocatable :: u(:), v(:) type(MPI_), intent(inout) :: MPID type(FEMDomain_), intent(inout) :: FEMDomain integer(int32), optional, allocatable, intent(in) :: fix_idx(:) real(real64), intent(in) :: dt logical, optional, intent(in) :: debug_mode real(real64), intent(in) :: cutoff_frequency real(real64) :: error_norm(1), all_error_norm(1) integer(int32) :: j, k ! [CAUTION!!] only undamped is implemented. du = u_n u = u_n v = 0.0d0*v_n do k = 1, this%itrmax !<<< du = this%OmegaSqMatrix%matmul(du) >>> for MPI du = FEMDomain%mpi_matmul(this%OmegaSqMatrix, du, MPID) error_norm = norm(LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*du) call MPID%AllReduce(sendobj=error_norm, recvobj=all_error_norm, count=1, sum=.true.) if (all_error_norm(1) < this%tol) then exit end if u = u + LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*du if (k == 1) then v = -dt*du else v = v - LPF_t_sinc_sqrt_taylor_coefficient(k=k - 1, t=dt, f_c=cutoff_frequency)*du end if end do deallocate (du) dv = v_n ! k=1 u = u + dt*dv v = v + dv do k = 1, this%itrmax !<<< dv = this%OmegaSqMatrix%matmul(dv) >>> for MPI dv = FEMDomain%mpi_matmul(this%OmegaSqMatrix, dv, MPID) error_norm = norm(LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*dv) call MPID%AllReduce(sendobj=error_norm, recvobj=all_error_norm, count=1, sum=.true.) if (all_error_norm(1) < this%tol) then exit end if u = u + LPF_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*dv v = v + LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*dv end do ! RHS (=M^{-1}F) if (present(RHS)) then if (allocated(this%weights)) then u = u + this%OmegaSqMatrix%tensor_wave_kernel_RHS(RHS=RHS/this%Mmatrix_diag, t=dt, tol=dble(1.0e-25), & itrmax=this%itrmax, cutoff_frequency=cutoff_frequency, weights=this%weights) v = v + LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=RHS/this%Mmatrix_diag, & DampingRatio=this%DampingRatio, itrmax=this%itrmax, tol=this%tol, weights=this%weights) else u = u + this%OmegaSqMatrix%tensor_wave_kernel_RHS(RHS=RHS/this%Mmatrix_diag, t=dt, tol=dble(1.0e-25), & itrmax=this%itrmax, cutoff_frequency=cutoff_frequency) v = v + LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=RHS/this%Mmatrix_diag, & DampingRatio=this%DampingRatio, itrmax=this%itrmax, tol=this%tol) end if end if if (present(fix_idx)) then if (allocated(fix_idx)) then u(fix_idx) = 0.0d0 end if end if end subroutine ! ############################################################## ! ############################################################## subroutine getDisplacement_and_Velocity_complex64_WaveKernel(this, u_n, v_n, dt, & fix_idx, cutoff_frequency, debug_mode, u, v, RHS) class(WaveKernel_), intent(inout) :: this complex(real64), intent(in) :: u_n(:), v_n(:) real(real64), optional, intent(in) :: RHS(:) complex(real64), allocatable :: du(:), dv(:) complex(real64), allocatable :: u(:), v(:) integer(int32), optional, allocatable, intent(in) :: fix_idx(:) real(real64), intent(in) :: dt logical, optional, intent(in) :: debug_mode real(real64), optional, intent(in) :: cutoff_frequency integer(int32) :: j, k ! [CAUTION!!] only undamped is implemented. du = u_n u = u_n v = 0.0d0*v_n do k = 1, this%itrmax du = this%OmegaSqMatrix%matmul(du) if (norm(LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*du) & < this%tol) exit u = u + LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*du if (k == 1) then v = -dt*du else v = v - LPF_t_sinc_sqrt_taylor_coefficient(k=k - 1, t=dt, f_c=cutoff_frequency)*du end if end do deallocate (du) dv = v_n ! k=1 u = u + dt*dv v = v + dv do k = 1, this%itrmax dv = this%OmegaSqMatrix%matmul(dv) if (norm(LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*dv) & < this%tol) exit u = u + LPF_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*dv v = v + LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=cutoff_frequency)*dv end do if (present(RHS)) then if (allocated(this%weights)) then u = u + this%OmegaSqMatrix%tensor_wave_kernel_RHS(RHS=RHS/this%Mmatrix_diag, t=dt, tol=dble(1.0e-25), & itrmax=this%itrmax, cutoff_frequency=cutoff_frequency, weights=this%weights) v = v + LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=RHS/this%Mmatrix_diag, & DampingRatio=this%DampingRatio, itrmax=this%itrmax, tol=this%tol, weights=this%weights) else u = u + this%OmegaSqMatrix%tensor_wave_kernel_RHS(RHS=RHS/this%Mmatrix_diag, t=dt, tol=dble(1.0e-25), & itrmax=this%itrmax, cutoff_frequency=cutoff_frequency) v = v + LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=RHS/this%Mmatrix_diag, & DampingRatio=this%DampingRatio, itrmax=this%itrmax, tol=this%tol) end if end if if (present(fix_idx)) then if (allocated(fix_idx)) then u(fix_idx) = 0.0d0 end if end if end subroutine ! ############################################################## ! ############################################################## function getDisplacementWaveKernel(this, u_n, v_n, dt, fix_idx, cutoff_frequency, debug_mode, RHS, gain) result(u) class(WaveKernel_), intent(inout) :: this real(real64), intent(in) :: u_n(:), v_n(:) real(real64), optional, intent(in) :: RHS(:) real(real64), allocatable :: u(:) integer(int32), optional, intent(in) :: fix_idx(:) real(real64), intent(in) :: dt logical, optional, intent(in) :: debug_mode real(real64), optional, intent(in) :: cutoff_frequency, gain real(real64) :: gain_value integer(int32) :: j real(real64) :: h, ddt gain_value = input(default=1.0d0, option=gain) if (present(cutoff_frequency)) then ! if(present(debug_mode) )then ! if(debug_mode)then ! ddt = 1.0d0/(cutoff_frequency*4.0d0) ! ! Hanning Window ! ! u = 0.250d0*exp(-this%DampingRatio/2.0d0*(dt-ddt) )*this%OmegaSqMatrix%tensor_wave_kernel(& ! u0=u_n,v0=v_n,t=dt-ddt,& ! itrmax=this%itrmax,tol=this%tol,fix_idx=fix_idx)& ! + 0.50d0*exp(-this%DampingRatio/2.0d0*(dt) )*this%OmegaSqMatrix%tensor_wave_kernel(& ! u0=u_n,v0=v_n,t=dt,& ! itrmax=this%itrmax,tol=this%tol,fix_idx=fix_idx)& ! + 0.250d0*exp(-this%DampingRatio/2.0d0*(dt-ddt) )*this%OmegaSqMatrix%tensor_wave_kernel(& ! u0=u_n,v0=v_n,t=dt+ddt,& ! itrmax=this%itrmax,tol=this%tol,fix_idx=fix_idx) ! return ! endif ! endif if (allocated(this%weights)) then u = LPF_Damped_cos_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, u_n=u_n, DampingRatio=this%DampingRatio, & fix_idx=fix_idx, itrmax=this%itrmax, tol=this%tol, weights=this%weights) & + LPF_Damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=v_n, DampingRatio=this%DampingRatio, & fix_idx=fix_idx, itrmax=this%itrmax, tol=this%tol, weights=this%weights) else u = LPF_Damped_cos_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, u_n=u_n, DampingRatio=this%DampingRatio, & fix_idx=fix_idx, itrmax=this%itrmax, tol=this%tol) & + LPF_Damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=v_n, DampingRatio=this%DampingRatio, & fix_idx=fix_idx, itrmax=this%itrmax, tol=this%tol) end if ! cutoff = - 10 dB else ! cutoff周波数なし，ただのWKF u = this%OmegaSqMatrix%tensor_wave_kernel( & u0=u_n, v0=v_n, t=dt, h=0.0d0, & itrmax=this%itrmax, tol=this%tol, fix_idx=fix_idx) u = exp(-this%DampingRatio*dt)*u end if if (present(RHS)) then if (allocated(this%weights)) then u = u + this%OmegaSqMatrix%tensor_wave_kernel_RHS(RHS=RHS/this%Mmatrix_diag, t=dt, tol=dble(1.0e-25), & itrmax=this%itrmax, cutoff_frequency=cutoff_frequency, weights=this%weights) else u = u + this%OmegaSqMatrix%tensor_wave_kernel_RHS(RHS=RHS/this%Mmatrix_diag, t=dt, tol=dble(1.0e-25), & itrmax=this%itrmax, cutoff_frequency=cutoff_frequency) end if end if u = gain_value*u end function ! ############################################################## ! ############################################################## function getVelocityWaveKernel(this, u_n, v_n, RHS, dt, fix_idx, cutoff_frequency, debug_mode, gain) result(v) class(WaveKernel_), intent(inout) :: this real(real64), intent(in) :: u_n(:), v_n(:) real(real64), optional, intent(in) :: RHS(:) integer(int32), optional, intent(in) :: fix_idx(:) real(real64), optional, intent(in) :: cutoff_frequency, gain logical, optional, intent(in) :: debug_mode real(real64), allocatable :: v(:) real(real64), intent(in) :: dt integer(int32) :: j real(real64) :: h, ddt, gain_value gain_value = input(default=1.0d0, option=gain) if (present(cutoff_frequency)) then ! if(present(debug_mode) )then ! if( debug_mode )then ! ! cutoff = - 3 dB ! ddt = 1.0d0/(cutoff_frequency*4.0d0) ! ! v = 0.250d0*this%OmegaSqMatrix%tensor_wave_kernel(& ! u0=-h*(dt-ddt)*u_n+v_n, & ! v0=-this%OmegaSqMatrix%matmul(u_n)-h*(dt-ddt)*v_n, & ! h=h,t=(dt-ddt),itrmax=this%itrmax,tol=this%tol,fix_idx=fix_idx)& ! + 0.500d0*this%OmegaSqMatrix%tensor_wave_kernel(& ! u0=-h*dt*u_n+v_n, & ! v0=-this%OmegaSqMatrix%matmul(u_n)-h*dt*v_n, & ! h=h,t=dt,itrmax=this%itrmax,tol=this%tol,fix_idx=fix_idx)& ! + 0.250d0*this%OmegaSqMatrix%tensor_wave_kernel(& ! u0=-h*(dt+ddt)*u_n+v_n, & ! v0=-this%OmegaSqMatrix%matmul(u_n)-h*(dt+ddt)*v_n, & ! h=h,t=(dt+ddt),itrmax=this%itrmax,tol=this%tol,fix_idx=fix_idx) ! ! return ! endif ! endif ! ! if( this%DampingRatio/=0.0d0 )then ! ! cutoff = - 3 dB ! ddt = 1.0d0/(cutoff_frequency*4.0d0) ! ! v = 0.250d0*this%OmegaSqMatrix%tensor_wave_kernel(& ! u0=-h*(dt-ddt)*u_n+v_n, & ! v0=-this%OmegaSqMatrix%matmul(u_n)-h*(dt-ddt)*v_n, & ! h=h,t=(dt-ddt),itrmax=this%itrmax,tol=this%tol,fix_idx=fix_idx)& ! + 0.500d0*this%OmegaSqMatrix%tensor_wave_kernel(& ! u0=-h*dt*u_n+v_n, & ! v0=-this%OmegaSqMatrix%matmul(u_n)-h*dt*v_n, & ! h=h,t=dt,itrmax=this%itrmax,tol=this%tol,fix_idx=fix_idx)& ! + 0.250d0*this%OmegaSqMatrix%tensor_wave_kernel(& ! u0=-h*(dt+ddt)*u_n+v_n, & ! v0=-this%OmegaSqMatrix%matmul(u_n)-h*(dt+ddt)*v_n, & ! h=h,t=(dt+ddt),itrmax=this%itrmax,tol=this%tol,fix_idx=fix_idx) ! return ! endif if (allocated(this%weights)) then v = LPF_damped_cos_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, u_n=v_n, DampingRatio=this%DampingRatio, & fix_idx=fix_idx, itrmax=this%itrmax, tol=this%tol, weights=this%weights) & + LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=-this%OmegaSqMatrix%matmul(u_n), & DampingRatio=this%DampingRatio, fix_idx=fix_idx, itrmax=this%itrmax, tol=this%tol, & weights=this%weights) else v = LPF_damped_cos_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, u_n=v_n, DampingRatio=this%DampingRatio, & fix_idx=fix_idx, itrmax=this%itrmax, tol=this%tol) & + LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=-this%OmegaSqMatrix%matmul(u_n), & DampingRatio=this%DampingRatio, fix_idx=fix_idx, itrmax=this%itrmax, tol=this%tol) end if else ! cutoff周波数なし，ただのWKF v = exp(-this%DampingRatio*dt)*this%OmegaSqMatrix%tensor_wave_kernel( & u0=-this%DampingRatio*dt*u_n + v_n, & v0=-this%OmegaSqMatrix%matmul(u_n) - this%DampingRatio*dt*v_n, & h=0.0d0, t=dt, itrmax=this%itrmax, tol=this%tol, fix_idx=fix_idx) end if if (present(RHS)) then if (allocated(this%weights)) then v = v + LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=RHS/this%Mmatrix_diag, & DampingRatio=this%DampingRatio, itrmax=this%itrmax, tol=this%tol, weights=this%weights) else v = v + LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix=this%OmegaSqMatrix, & dt=dt, f_c=cutoff_frequency, v_n=RHS/this%Mmatrix_diag, & DampingRatio=this%DampingRatio, itrmax=this%itrmax, tol=this%tol) end if end if v = gain_value*v end function ! ############################################################## ! ############################################################## !subroutine hanning_WaveKernel(this,u_n, v_n,cutoff_frequency,dt,fix_idx) ! class(WaveKernel_),intent(inout) :: this ! real(real64),intent(inout) :: u_n(:),v_n(:) ! real(real64),intent(in) :: cutoff_frequency ! integer(int32),optional,intent(in) :: fix_idx(:) ! real(real64),allocatable :: u(:) ! real(real64),intent(in) :: dt ! integer(int32) :: j ! real(real64) :: h ! ! if(.not.allocated(this%hanning_buffer) )then ! this%hanning_buffer = zeros(size(u_n),this%hanning_buffer_size ) ! do j=1,this%hanning_buffer_size ! this%hanning_buffer(:,j) = u_n(:) ! enddo ! endif ! ! ! if(.not.allocated(this%u_in1) )then ! this%u_in1 = u_n ! endif ! ! if(.not.allocated(this%v_in1) )then ! this%v_in1 = v_n ! endif ! ! u_n(fix_idx)=0.0d0 ! !$OMP parallel do ! do j=1,size(u_n) ! u_n(j) = hanning_filter(x_n=u_n(j),& ! fc=cutoff_frequency,& ! sampling_Hz=1.0d0/dt,& ! buf=this%hanning_buffer(:,j) & ! ) ! enddo ! !$OMP end parallel do ! ! v_n(fix_idx)=0.0d0 ! !$OMP parallel do ! do j=1,size(v_n) ! v_n(j) = hanning_filter(x_n=v_n(j),& ! fc=cutoff_frequency,& ! sampling_Hz=1.0d0/dt,& ! buf=this%hanning_buffer(:,j) & ! ) ! enddo ! !$OMP end parallel do ! ! ! !end subroutine !! ############################################################## ! !function hanning_filter_wavekernel(x_n,fc,sampling_Hz,buf) result(ret) ! real(real64),intent(in) :: x_n,fc,sampling_Hz,buf(:) ! real(real64) :: ret ! real(real64) :: t,period_T,dt ! integer(int32) :: i ! type(Math_) :: math ! ! ret = 0.0d0*x_n ! dt = sampling_Hz ! t = - dt*size(buf,1) ! period_T = 1.0d0/fc ! do i=1,size(buf,1) ! t = t + dt ! ret = ret + buf(i)*(0.50d0 - 0.50d0*cos(2.0d0*math%pi*t/period_T) ) ! enddo ! !end function pure function LPF_cos_sqrt_taylor_coefficient(k, t, f_c) result(c_k) integer(int32), intent(in) :: k real(real64), intent(in) :: t real(real64), intent(in) :: f_c real(real64) :: c_k real(real64) :: t_hat, dt type(Math_) :: math if (k == 0) then c_k = 1.0d0 return end if !dt = 1.0d0/(f_c*4.0d0) dt = 1.0d0/f_c/2.0d0/math%pi*acos(sqrt(2.0d0) - 1.0d0) t_hat = 0.250d0*((t - dt)**(2*k)) + 0.50d0*((t)**(2*k)) + 0.250d0*((t + dt)**(2*k)) if (mod(k, 2) == 0) then ! k=even c_k = t_hat/Gamma(2.0d0*k + 1.0d0) else ! k=odd c_k = -1.0d0*t_hat/Gamma(2.0d0*k + 1.0d0) end if end function function LPF_damped_cos_sqrt_taylor_coefficient(k, t, f_c, DampingRatio, weights) result(c_k) integer(int32), intent(in) :: k real(real64), intent(in) :: t real(real64), intent(in) :: f_c real(real64), intent(in) :: DampingRatio(:) real(real64), optional, intent(in) :: weights(:) real(real64) :: dt real(real64), allocatable:: c_k(:) real(real64), allocatable:: t_tilde(:) type(Math_) :: math integer(int32) :: n integer(int32) :: lbd, ubd real(real64), allocatable :: w(:) if (present(weights)) then ubd = size(weights) - (size(weights) + 1)/2 lbd = -(size(weights) - (size(weights) + 1)/2) allocate (w(lbd:ubd)) w(lbd:ubd) = weights(1:size(weights)) end if if (k == 0) then c_k = ones(size(DampingRatio)) return end if t_tilde = zeros(size(DampingRatio)) !dt = 1.0d0/(f_c*4.0d0) dt = 1.0d0/f_c/2.0d0/math%pi*acos(sqrt(2.0d0) - 1.0d0) if (present(weights)) then t_tilde = 0.0d0 do n = lbound(w, 1), ubound(w, 1) t_tilde = t_tilde + w(n)*((t + n*dt)**(2*k))*exp(-DampingRatio*(t + n*dt)) end do else t_tilde = 0.250d0*((t - dt)**(2*k))*exp(-DampingRatio*(t - dt)) & + 0.500d0*((t)**(2*k))*exp(-DampingRatio*(t)) & + 0.250d0*((t + dt)**(2*k))*exp(-DampingRatio*(t + dt)) end if if (mod(k, 2) == 0) then ! k=even c_k = t_tilde/Gamma(2.0d0*k + 1.0d0) else ! k=odd c_k = -1.0d0*t_tilde/Gamma(2.0d0*k + 1.0d0) end if end function pure function LPF_t_sinc_sqrt_taylor_coefficient(k, t, f_c) result(s_k) integer(int32), intent(in) :: k real(real64), intent(in) :: t real(real64), intent(in) :: f_c real(real64) :: s_k real(real64) :: t_hat, dt type(Math_) :: math if (k == 0) then s_k = t return end if !dt = 1.0d0/(f_c*4.0d0) dt = 1.0d0/f_c/2.0d0/math%pi*acos(sqrt(2.0d0) - 1.0d0) t_hat = 0.250d0*((t - dt)**(2*k + 1)) + 0.50d0*((t)**(2*k + 1)) & + 0.250d0*((t + dt)**(2*k + 1)) if (mod(k, 2) == 0) then ! k=even s_k = (2.0d0**(-1 - 2*k))*sqrt(math%pi)/Gamma(1.0d0 + dble(k))/Gamma(1.5d0 + dble(k)) & *t_hat else ! k=odd s_k = -1.0d0* & (2.0d0**(-1 - 2*k))*sqrt(math%pi)/Gamma(1.0d0 + dble(k))/Gamma(1.5d0 + dble(k)) & *t_hat end if end function function LPF_damped_t_sinc_sqrt_taylor_coefficient(k, t, f_c, DampingRatio, weights) result(s_k) integer(int32), intent(in) :: k real(real64), intent(in) :: t real(real64), intent(in) :: f_c real(real64), allocatable :: s_k(:) real(real64), allocatable :: t_tilde(:) real(real64), optional, intent(in) :: weights(:) real(real64) :: dt, buf real(real64), intent(in) :: DampingRatio(:) type(Math_) :: math integer(int32) :: n integer(int32) :: lbd, ubd real(real64), allocatable :: w(:) if (present(weights)) then ubd = size(weights) - (size(weights) + 1)/2 lbd = -(size(weights) - (size(weights) + 1)/2) allocate (w(lbd:ubd)) w(lbd:ubd) = weights(1:) end if if (k == 0) then s_k = t*ones(size(DampingRatio)) return end if t_tilde = zeros(size(DampingRatio)) !dt = 1.0d0/(f_c*4.0d0) dt = 1.0d0/f_c/2.0d0/math%pi*acos(sqrt(2.0d0) - 1.0d0) if (present(weights)) then t_tilde = 0.0d0 do n = lbound(w, 1), ubound(w, 1) t_tilde = t_tilde + w(n)*((t + n*dt)**(2*k + 1))*exp(-DampingRatio*(t + n*dt)) end do ! check solution !if(lbound(weights,1)==-1 .and. ubound(weights,1)==1)then ! if(weights(-1)==0.250d0 .or. weights(0)==0.50d0)then ! if(weights(0)==0.250d0)then ! buf = 0.250d0*((t-dt)**(2*k+1) )*exp(-DampingRatio*(t-dt) ) & ! + 0.500d0*((t )**(2*k+1) )*exp(-DampingRatio*(t ) ) & ! + 0.250d0*((t+dt)**(2*k+1) )*exp(-DampingRatio*(t+dt) ) ! if(buf/=t_tilde)then ! open(10,file="debug.txt") ! write(10,*) buf,t_tilde ! close(10) ! endif ! endif ! endif !endif ! check solution else t_tilde = 0.250d0*((t - dt)**(2*k + 1))*exp(-DampingRatio*(t - dt)) & + 0.500d0*((t)**(2*k + 1))*exp(-DampingRatio*(t)) & + 0.250d0*((t + dt)**(2*k + 1))*exp(-DampingRatio*(t + dt)) end if if (mod(k, 2) == 0) then ! k=even s_k = (2.0d0**(-1 - 2*k))*sqrt(math%pi)/Gamma(1.0d0 + dble(k))/Gamma(1.5d0 + dble(k)) & *t_tilde else ! k=odd s_k = -1.0d0* & (2.0d0**(-1 - 2*k))*sqrt(math%pi)/Gamma(1.0d0 + dble(k))/Gamma(1.5d0 + dble(k)) & *t_tilde end if end function ! ######################################################### function LPF_cos_sqrt_WaveKernelFunction(Omega_sq_matrix, dt, f_c, u_n, itrmax, tol, fix_idx) result(ret) type(CRS_), intent(inout) :: Omega_sq_matrix real(real64), intent(in) :: dt, f_c, u_n(:), tol real(real64), allocatable :: ret(:), du(:) integer(int32), intent(in) :: itrmax, fix_idx(:) integer(int32) :: k k = 0 du = u_n du(fix_idx) = 0.0d0 ret = LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c)*du do k = 1, itrmax du = Omega_sq_matrix%matmul(du) if (norm(LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c)*du) < tol) exit ret = ret + LPF_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c)*du end do end function ! ######################################################### ! ######################################################### function LPF_damped_cos_sqrt_WaveKernelFunction(Omega_sq_matrix,dt,f_c,u_n,DampingRatio, itrmax,tol,fix_idx,weights) result(ret) type(CRS_), intent(inout) :: Omega_sq_matrix real(real64), intent(in) :: dt, f_c, u_n(:), tol, DampingRatio(:) real(real64), allocatable :: ret(:), du(:) integer(int32), intent(in) :: itrmax integer(int32), optional, intent(in) :: fix_idx(:) real(real64), optional, intent(in) :: weights(:) integer(int32) :: k k = 0 du = u_n if (present(fix_idx)) then du(fix_idx) = 0.0d0 end if ret = LPF_damped_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*du do k = 1, itrmax du = Omega_sq_matrix%matmul(du) if (norm(LPF_damped_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*du) < tol) exit ret = ret + LPF_damped_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*du end do end function ! ######################################################### ! ######################################################### function LPF_damped_cos_sqrt_complex64_WaveKernelFunction(Omega_sq_matrix,dt,f_c,u_n,DampingRatio, itrmax,tol,fix_idx) result(ret) type(CRS_), intent(inout) :: Omega_sq_matrix real(real64), intent(in) :: dt, f_c, tol, DampingRatio(:) complex(real64), intent(in) :: u_n(:) complex(real64), allocatable :: ret(:), du(:) integer(int32), intent(in) :: itrmax integer(int32), optional, intent(in) :: fix_idx(:) integer(int32) :: k k = 0 du = u_n if (present(fix_idx)) then du(fix_idx) = 0.0d0 end if ret = LPF_damped_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio)*du do k = 1, itrmax du = Omega_sq_matrix%matmul(du) if (norm(LPF_damped_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio)*du) < tol) exit ret = ret + LPF_damped_cos_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio)*du end do end function ! ######################################################### ! ######################################################### function LPF_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix, dt, f_c, v_n, itrmax, tol, fix_idx) result(ret) type(CRS_), intent(inout) :: Omega_sq_matrix real(real64), intent(in) :: dt, f_c, v_n(:), tol real(real64), allocatable :: ret(:), dv(:) integer(int32), intent(in) :: itrmax, fix_idx(:) integer(int32) :: k k = 0 dv = v_n dv(fix_idx) = 0.0d0 ret = LPF_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c)*dv do k = 1, itrmax dv = Omega_sq_matrix%matmul(dv) if (norm(LPF_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c)*dv) < tol) exit ret = ret + LPF_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c)*dv end do end function ! ######################################################### ! ######################################################### function LPF_damped_t_sinc_sqrt_WaveKernelFunction(Omega_sq_matrix,dt,f_c,v_n,itrmax,tol,fix_idx,DampingRatio,weights) result(ret) type(CRS_), intent(inout) :: Omega_sq_matrix real(real64), intent(in) :: dt, f_c, v_n(:), tol, DampingRatio(:) real(real64), allocatable :: ret(:), dv(:) integer(int32), intent(in) :: itrmax integer(int32), optional, intent(in) :: fix_idx(:) real(real64), optional, intent(in) :: weights(:) integer(int32) :: k k = 0 dv = v_n if (present(fix_idx)) then dv(fix_idx) = 0.0d0 end if ret = LPF_damped_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*dv do k = 1, itrmax dv = Omega_sq_matrix%matmul(dv) if (norm(LPF_damped_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*dv) < tol) exit ret = ret + LPF_damped_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*dv end do end function ! ######################################################### ! ######################################################### function LPF_damped_t_sinc_sqrt_MPI_WaveKernelFunction(Omega_sq_matrix, & dt, f_c, v_n, itrmax, tol, fix_idx, DampingRatio, MPID, FEMDomain, weights) result(ret) type(CRS_), intent(inout) :: Omega_sq_matrix real(real64), intent(in) :: dt, f_c, v_n(:), tol, DampingRatio(:) type(MPI_), intent(inout) :: MPID type(FEMDomain_), intent(inout) :: FEMDomain real(real64), allocatable :: ret(:), dv(:) integer(int32), intent(in) :: itrmax integer(int32), optional, intent(in) :: fix_idx(:) real(real64), optional, intent(in) :: weights(:) integer(int32) :: k k = 0 dv = v_n if (present(fix_idx)) then dv(fix_idx) = 0.0d0 end if ret = LPF_damped_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*dv do k = 1, itrmax call MPID%barrier() !dv = Omega_sq_matrix%matmul(dv) dv = FEMDomain%mpi_matmul(Omega_Sq_Matrix, dv, MPID) if (norm(LPF_damped_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*dv) < tol) then call MPID%barrier() exit end if ret = ret + LPF_damped_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*dv end do call MPID%barrier() end function ! ######################################################### ! ######################################################### function LPF_damped_t_sinc_sqrt_complex64_WaveKernelFunction( & Omega_sq_matrix, dt, f_c, v_n, itrmax, tol, fix_idx, DampingRatio, weights) result(ret) type(CRS_), intent(inout) :: Omega_sq_matrix real(real64), intent(in) :: dt, f_c, tol, DampingRatio(:) complex(real64), intent(in) :: v_n(:) complex(real64), allocatable :: ret(:), dv(:) integer(int32), intent(in) :: itrmax integer(int32), optional, intent(in) :: fix_idx(:) real(real64), optional, intent(in) :: weights(:) integer(int32) :: k k = 0 dv = v_n if (present(fix_idx)) then dv(fix_idx) = 0.0d0 end if ret = LPF_damped_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*dv do k = 1, itrmax dv = Omega_sq_matrix%matmul(dv) if (norm(LPF_damped_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*dv) < tol) exit ret = ret + LPF_damped_t_sinc_sqrt_taylor_coefficient(k=k, t=dt, f_c=f_c, DampingRatio=DampingRatio, weights=weights)*dv end do end function ! ######################################################### subroutine setLPF_WaveKernel(this, weights) class(WaveKernel_), intent(inout) :: this real(real64), intent(in) :: weights(:) if (mod(size(weights), 2) == 0) then print *, "[ERROR] >> setLPF_WaveKernel >> size(weights) should be odd." else this%weights = weights end if end subroutine end module