Merge pull request #556 from kazulagi/work

bugfix >> SparseClass >> tensor_wave_kernel_RHS_LPF_real_64_crs >> t_…
kazulagi · Jan 27, 2025 · e558ffe · e558ffe
2 parents 42fdcc4 + 0761c3d
commit e558ffe
Show file tree

Hide file tree

Showing 3 changed files with 14 additions and 2 deletions.
diff --git 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
@@ -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
@@ -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(:)