Final version of the Helmholtz equation

docs/src/matrix_free.md

@@ -326,6 +326,7 @@ Combining these, the discretized Helmholtz equation becomes:
 ```
 
 This discretization results in an equation at each grid point, resulting in a large and sparse linear system when assembled across the entire 3D grid.
+To simplify the example, we impose Dirichlet boundary conditions with the solution $u(x, y, z) = 0$ on the boundary of the cubic domain.
 
 Explicitly constructing this large sparse matrix is often impractical and unnecessary.
 Instead, we can define a function that directly applies the Helmholtz operator to the 3D grid, avoiding the need to form the matrix explicitly.
@@ -339,9 +340,9 @@ using Krylov, LinearAlgebra
 
 # Parameters
 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
+Nx = 200                # Number of interior grid points in x
+Ny = 200                # Number of interior grid points in y
+Nz = 200                # 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
@@ -382,7 +383,6 @@ function LinearAlgebra.mul!(y::Vector, A::HelmholtzOperator, u::Vector)
                 elseif i == A.Nx
                     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
                 end
 
@@ -413,28 +413,21 @@ end
 # Create the matrix-free operator for the Helmholtz equation
 A = HelmholtzOperator(Nx, Ny, Nz, Δx, Δy, Δz, k)
 
-# Create the right-hand side with the source term
-# f(x, y, z) = -2k² * sin(kx) * sin(ky) * sin(kz)
-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+1]) * sin(k * y[jj+1]) * sin(k * z[kk+1])
-        end
-    end
-end
-
-# Vectorize the right-hand side
+# Source term f(x, y, z) = -2k² * sin(kx) * sin(ky) * sin(kz)
+F = [-2 * k^2 * sin(k * x[ii+1]) * sin(k * y[jj+1]) * sin(k * z[kk+1]) for ii in 1:Nx, jj in 1:Ny, kk in 1:Nz]
 f = vec(F)
 
 # 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
+# Solution as 3D array
 U_sol = reshape(u_sol, Nx, Ny, Nz)
 
-# 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])
+# Exact solution u(x,y,z) = sin(kx) * sin(ky) * sin(kz)
+U_star = [sin(k * x[ii+1]) * sin(k * y[jj+1]) * sin(k * z[kk+1]) for ii in 1:Nx, jj in 1:Ny, kk in 1:Nz]
+
+# Compute the maximum error between the numerical solution U_sol and the exact solution U_star
+norm(U_sol - U_star, Inf)
 ```
 
 Note that preconditioners can be also implemented as abstract operators.