Update matrix_free.md

docs/src/matrix_free.md

@@ -338,18 +338,20 @@ This approach enables a simpler and efficient application of the operator in 3D
 using Krylov, LinearAlgebra
 
 # Parameters
-L = 1.0                  # Length of the cubic domain
-Nx, Ny, Nz = 40, 40, 40  # Number of grid points
-Δx = L / (Nx - 1)        # Grid spacing in x
-Δy = L / (Ny - 1)        # Grid spacing in y
-Δz = L / (Nz - 1)        # Grid spacing in z
-wavelength = 0.5         # Wavelength of the wave
-k = 2 * π / wavelength   # Wave number
+L = 1.0                 # Length of the cubic domain
+Nx = 400                # Number of interior grid points in x
+Ny = 400                # Number of interior grid points in y
+Nz = 400                # Number of interior grid points in z
+Δx = L / (Nx + 1)       # Grid spacing in x
+Δy = L / (Ny + 1)       # Grid spacing in y
+Δz = L / (Nz + 1)       # Grid spacing in z
+wavelength = 0.5        # Wavelength of the wave
+k = 2 * π / wavelength  # Wave number
 
 # Create the grid points
-x = 0:Δx:L  # Points in x dimension
-y = 0:Δy:L  # Points in y dimension
-z = 0:Δz:L  # Points in z dimension
+x = 0:Δx:L  # Points in x dimension (Nx + 2)
+y = 0:Δy:L  # Points in y dimension (Ny + 2)
+z = 0:Δz:L  # Points in z dimension (Nz + 2)
 
 # Define a matrix-free Helmholtz operator
 struct HelmholtzOperator
@@ -376,39 +378,30 @@ function LinearAlgebra.mul!(y::Vector, A::HelmholtzOperator, u::Vector)
         for j in 1:A.Ny
             for k in 1:A.Nz
                 if i == 1
-                    # Forward difference on x boundary
-                    dx2 = (U[i+2,j,k] - 2 * U[i+1,j,k] + U[i,j,k]) / (A.Δx)^2
+                    dx2 = (U[i+1,j,k] -2 * U[i,j,k]) / (A.Δx)^2
                 elseif i == A.Nx
-                    # Backward difference on x boundary
-                    dx2 = (U[i,j,k] - 2 * U[i-1,j,k] + U[i-2,j,k]) / (A.Δx)^2
+                    dx2 = (-2 * U[i,j,k] + U[i-1,j,k]) / (A.Δx)^2
                 else
                     # Centered difference for interior points
-                    dx2 = (U[i+1,j,k] - 2 * U[i,j,k] + U[i-1,j,k]) / (A.Δx)^2
+                    dx2 = (U[i+1,j,k] -2 * U[i,j,k] + U[i-1,j,k]) / (A.Δx)^2
                 end
 
                 if j == 1
-                    # Forward difference on y boundary
-                    dy2 = (U[i,j+2,k] - 2 * U[i,j+1,k] + U[i,j,k]) / (A.Δy)^2
+                    dy2 = (U[i,j+1,k] -2 * U[i,j,k]) / (A.Δy)^2
                 elseif j == A.Ny
-                    # Backward difference on y boundary
-                    dy2 = (U[i,j,k] - 2 * U[i,j-1,k] + U[i,j-2,k]) / (A.Δy)^2
+                    dy2 = (-2 * U[i,j,k] + U[i,j-1,k]) / (A.Δy)^2
                 else
-                    # Centered difference for interior points
-                    dy2 = (U[i,j+1,k] - 2 * U[i,j,k] + U[i,j-1,k]) / (A.Δy)^2
+                    dy2 = (U[i,j+1,k] -2 * U[i,j,k] + U[i,j-1,k]) / (A.Δy)^2
                 end
 
                 if k == 1
-                    # Forward difference on z boundary
-                    dz2 = (U[i,j,k+2] - 2 * U[i,j,k+1] + U[i,j,k]) / (A.Δz)^2
+                    dz2 = (U[i,j,k+1] -2 * U[i,j,k]) / (A.Δz)^2
                 elseif k == A.Nz
-                    # Backward difference on z boundary
-                    dz2 = (U[i,j,k] - 2 * U[i,j,k-1] + U[i,j,k-2]) / (A.Δz)^2
+                    dz2 = (-2 * U[i,j,k] + U[i,j,k-1]) / (A.Δz)^2
                 else
-                    # Centered difference for interior points
-                    dz2 = (U[i,j,k+1] - 2 * U[i,j,k] + U[i,j,k-1]) / (A.Δz)^2
+                    dz2 = (U[i,j,k+1] -2 * U[i,j,k] + U[i,j,k-1]) / (A.Δz)^2
                 end
 
-                # Helmholtz operator application
                 Y[i,j,k] = dx2 + dy2 + dz2 + (A.k)^2 * U[i,j,k]
             end
         end
@@ -426,31 +419,22 @@ F = zeros(Nx, Ny, Nz)
 for ii in 1:Nx
     for jj in 1:Ny
         for kk in 1:Nz
-            F[ii,jj,kk] = -2 * k^2 * sin(k * x[ii]) * sin(k * y[jj]) * sin(k * z[kk])
+            F[ii,jj,kk] = -2 * k^2 * sin(k * x[ii+1]) * sin(k * y[jj+1]) * sin(k * z[kk+1])
         end
     end
 end
 
 # Vectorize the right-hand side
 f = vec(F)
 
-# Solve the linear system using GMRES
-u_sol, stats = gmres(A, f, atol=1e-10, rtol=0.0, verbose=1)
+# Solve the linear system using MinAres
+u_sol, stats = minares(A, f, atol=1e-10, rtol=0.0, verbose=1)
 
 # Reshape the solution back to a 3D array for visualization or further processing
 U_sol = reshape(u_sol, Nx, Ny, Nz)
 
-# Compute the exact solution u(x,y,z) = sin(kx) * sin(ky) * sin(kz)
-U_star = zeros(Nx, Ny, Nz)
-for ii in 1:Nx
-    for jj in 1:Ny
-        for kk in 1:Nz
-            U_star[ii,jj,kk] = sin(k * x[ii]) * sin(k * y[jj]) * sin(k * z[kk])
-        end
-    end
-end
-
-norm(U_sol - U_star)
+# Compare U_sol with the exact solution u(x,y,z) = sin(kx) * sin(ky) * sin(kz)
+# U_sol[ii,jj,kk] ≈ sin(k * x[ii+1]) * sin(k * y[jj+1]) * sin(k * z[kk+1])
 ```
 
 Note that preconditioners can be also implemented as abstract operators.