From 0761c3de505c7a7d03c8149790ad5655d50bdf9f Mon Sep 17 00:00:00 2001
From: kazulagi <you@example.com>
Date: Mon, 27 Jan 2025 16:04:42 +0900
Subject: [PATCH] bugfix >> SparseClass >>
 tensor_wave_kernel_RHS_LPF_real_64_crs >> t_hat should be initialized before
 assembled.

---
 src/SparseClass/SparseClass.f90                   | 12 ++++++++++++
 src/SpectreAnalysisClass/SpectreAnalysisClass.f90 |  2 +-
 src/WaveKernelClass/WaveKernelClass.f90           |  2 +-
 3 files changed, 14 insertions(+), 2 deletions(-)

diff --git a/src/SparseClass/SparseClass.f90 b/src/SparseClass/SparseClass.f90
index c21700ba..2dfb4a23 100755
--- a/src/SparseClass/SparseClass.f90
+++ b/src/SparseClass/SparseClass.f90
@@ -3693,6 +3693,7 @@ function tensor_wave_kernel_RHS_LPF_real_64_crs(this, RHS, t, tol, itrmax, fix_i
       real(real64), allocatable :: w(:)
 
       integer(int32) :: lbd, ubd
+      
 
       if (present(weights)) then
          ubd = size(weights) - (size(weights) + 1)/2
@@ -3720,11 +3721,13 @@ function tensor_wave_kernel_RHS_LPF_real_64_crs(this, RHS, t, tol, itrmax, fix_i
       !ddt = 1.0d0/4.0d0/cutoff_frequency
       ddt = 1.0d0/cutoff_frequency/2.0d0/math%pi*acos(sqrt(2.0d0) - 1.0d0)
 
+      
       !u = -1.0d0/2.0d0*t*t*RHS
 
       u = 0.0d0*RHS
 
       hat_t = 0.0d0
+      !print *,"dbg",lbound(w, 1), ubound(w, 1) 
       do m = lbound(w, 1), ubound(w, 1)
          hat_t = hat_t + w(m)*(t + m*ddt)**2
       end do
@@ -3735,21 +3738,30 @@ function tensor_wave_kernel_RHS_LPF_real_64_crs(this, RHS, t, tol, itrmax, fix_i
       du = -1.0d0/2.0d0*RHS
       u = u + hat_t*du
 
+
+      
+
       do n = 1, itr_max
+         
+         print *, "debug 2025_01_27",n,maxval(u),minval(u)
+         hat_t = 0.0d0
          do m = lbound(w, 1), ubound(w, 1)
             hat_t = hat_t + w(m)*(t + m*ddt)**(2*n + 2)
          end do
+
          !hat_t = 0.250d0*(t-ddt)**(2*n+2) &
          !     + 0.50d0*(t)**(2*n+2)&
          !     + 0.250d0*(t+ddt)**(2*n+2)
 
          du = -1.0d0/dble(2*n + 1)/dble(2*n + 2)*this%matmul(du)
          u = u + hat_t*du
+         
          if (present(debug)) then
             if (debug) then
                print *, n, norm(hat_t*du)
             end if
          end if
+
          if (norm(hat_t*du) < itr_tol) exit
       end do
 
diff --git a/src/SpectreAnalysisClass/SpectreAnalysisClass.f90 b/src/SpectreAnalysisClass/SpectreAnalysisClass.f90
index 2b3a20e0..ce12d450 100755
--- a/src/SpectreAnalysisClass/SpectreAnalysisClass.f90
+++ b/src/SpectreAnalysisClass/SpectreAnalysisClass.f90
@@ -1336,7 +1336,7 @@ function fk_SpectreAnalysis(this,t,x,position,k_axis) result(z)
 
                f = ((1.0d0/dt))/(num_fft)*(k-1)
                w = 2.0d0*math%PI*f
-               
+               if (abs(U(freq,m)) ==0.0d0) cycle
                z(freq,k) = z(freq,k) &
                   + U(freq,m)&
                      *exp(math%i    &
diff --git a/src/WaveKernelClass/WaveKernelClass.f90 b/src/WaveKernelClass/WaveKernelClass.f90
index 79ec9edf..39129ab2 100644
--- a/src/WaveKernelClass/WaveKernelClass.f90
+++ b/src/WaveKernelClass/WaveKernelClass.f90
@@ -365,7 +365,7 @@ subroutine movingAverage_WaveKernel(this, u_n, v_n, fix_idx)
 
 ! ##############################################################
    subroutine getDisplacement_and_Velocity_WaveKernel(this, u_n, v_n, dt, &
-                                                      fix_idx, cutoff_frequency, debug_mode, u, v, RHS)
+                              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(:)