Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add regularization to Laplacian and preconditioners for ConjugateGradientPoissonSolver #3848

Open
wants to merge 10 commits into
base: main
Choose a base branch
from
Open

Add regularization to Laplacian and preconditioners for ConjugateGradientPoissonSolver #3848

Changes from 1 commit
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
4f6ffd7
Add regularization to Poisson operator for CG solve
glwagner Oct 16, 2024
19f0d68
Add custom gauge condition to FFT solver
glwagner Oct 16, 2024
b62951c
Better implementation
glwagner Oct 16, 2024
d382923
Fix sign error
glwagner Oct 16, 2024
f132a5f
Cosmetic fix
glwagner Oct 16, 2024
42c3993
Add regularization consistently in preconditioner
glwagner Oct 16, 2024
b667af0
Add temporary interface for building fft solver
glwagner Oct 17, 2024
a12c2f0
Merge branch 'main' into glw/regularized-laplacian
glwagner Oct 17, 2024
c653542
Change name of preconditioner
glwagner Oct 17, 2024
3f62872
Change sign of regularization, plus cleaner code design
glwagner Oct 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 13 additions & 20 deletions src/Solvers/conjugate_gradient_poisson_solver.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,8 @@ end

@kernel function laplacian!(∇²ϕ, grid, ϕ)
i, j, k = @index(Global, NTuple)
@inbounds ∇²ϕ[i, j, k] = ∇²ᶜᶜᶜ(i, j, k, grid, ϕ)
active = !inactive_cell(i, j, k, grid)
@inbounds ∇²ϕ[i, j, k] = ∇²ᶜᶜᶜ(i, j, k, grid, ϕ) * active
end

struct RegularizedLaplacian{D}
Expand All @@ -47,7 +48,10 @@ function (L::RegularizedLaplacian)(Lϕ, ϕ)
if !isnothing(L.δ)
# Add regularization
ϕ̄ = mean(ϕ)
parent(Lϕ) .+= L.δ * ϕ̄
ΔLϕ = L.δ * ϕ̄
grid = ϕ.grid
arch = architecture(grid)
launch!(arch, grid, :xyz, subtract_and_mask!, Lϕ, grid, ΔLϕ)
end

return nothing
Expand Down Expand Up @@ -136,15 +140,19 @@ const FFTBasedPreconditioner = RegularizedPreconditioner{<:SolverWithFFT}
function precondition!(p, regularized::FFTBasedPreconditioner, r, args...)
solver = regularized.preconditioner
compute_preconditioner_rhs!(solver, r)
p = solve!(p, solver)
solve!(p, solver)
regularize_poisson_solution!(p, regularized)
return p
end

function regularize_poisson_solution!(p, regularized)
δ = regularized.regularization
rhs = regularized.rhs
mean_p = mean(p)
Copy link
Contributor

@Yixiao-Zhang Yixiao-Zhang Oct 18, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I may have found why the solver blows up for positive regularization.

Here, p is the solution in the whole grid by the FFT Poisson solver, but mean_p calculates the mean of p in the immersed boundary grid. In my test, mean_p is usually order 1e-5, which is much larger than round-off errors. I also calculated the eigenvalues of the preconditioner with regularization=0, and I found both negative and positive eigenvalues.

Removing the subtraction of mean_p enables convergence for positive regularization according to my test.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is interesting! I'll test it. I don't completely understand why this is what's needed, but let's think about it.

Copy link
Member Author

@glwagner glwagner Oct 18, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm still getting NaNs. What regularization did you use? Did you make any other changes? EDIT: the regularization does seem to matter, now I get non-blow-up with regularization ~ 1/L^2.

Copy link
Contributor

@Yixiao-Zhang Yixiao-Zhang Oct 18, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tested the code in a simpler case without FlatExtrapolationOpenBoundaryCondition. That may explain why the result is different. Besides, I usually check the distribution of eigenvalues to judge whether the solver is numerically stable than running the simulation. I am looking into the case with FlatExtrapolationOpenBoundaryCondition and aim to understand how it changes the distribution of eigenvalues.

I just ran your demo and got NaNs as well.

Copy link
Member Author

@glwagner glwagner Oct 18, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just to write out the idea, we are trying to solve a "regularized" Poisson equation:

$$ \tilde L \phi = R $$

where $\tilde L = \nabla^2 - \delta V^{-1} \int \mathrm{d V}$ is the regularized Laplace operator (on an immersed boundary grid) with the singularity removed.

We're proposing to use a preconditioner that directly inverts (with an FFT --- fast) the "non-immersed" (but regularized) Laplace operator,

$$ L_u = \nabla_u^2 - \delta V_\star^{-1} \int_{\Omega_\star} \mathrm{d V} $$

where $\nabla_u^2$ applies on the "underlying grid" (ie ignoring the immersed boundary). The present discussion is how to define the mean, which is taken over the domain $\Omega_\star$ with volume $V_\star$. The current code takes $V_\star = V$, but we can also consider taking $V_\star = V_u$ and $\Omega_\star = \Omega_u$.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

By the way @Yixiao-Zhang I think this example really meaningfully exposes the effect of regularization, because it's only with this kind of open boundary condition that $\bar R \ne 0$. I believe that with impenetrable (and periodic of course) boundary conditions $\bar R$ is very small.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did not realized $\overline{R} \neq 0$ for this kind of open boundary condition. However, I think the current pressure solver cannot handle FlatExtrapolationOpenBoundaryCondition properly. As far as I understand, the pressure solver is unaware of the open boundary condition, but the open boundary condition should affect how fill_halo_region! works for rhs. The original boundary condition for rhs is $\hat{n}\cdot \nabla r = 0$, but it should be $(\hat{n}\cdot \nabla)^2 r = 0$ for open boundary conditions? (I am not sure.)

Right now, rhs is defined in this way:

rhs = CenterField(grid)

The good news is that the Laplace operator becomes negative definite if we take the open boundary into account. So, there is no need to regularize the Poisson solver. (The preconditioner is still semi-definite, though.)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the pressure solver is unaware of the open boundary condition

This isn't quite true. The inhomogeneous part of any boundary condition is incorporated into the Poisson equation RHS (this happens implicitly by 1) filling halo regions on the predictor velocity prior to the pressure solve and then 2) taking the divergence of the predictor velocity.) This is why the operator / FFT-solve must use homogeneous Neumann conditions for the algorithm to work.


if !isnothing(δ)
mean_rhs = mean(rhs)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In principle, the preconditioner solves the Poisson equation in the underlying grid. Following this logic, it is the mean of rhs calculated in the underlying grid that should be subtracted. I have not tested which is better, but I feel the mean in the underlying grid is more consistent mathematically. (I assume mean(rhs) here is the mean of rhs calculated in the immersed boundary grid.)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Correct, mean(rhs) accounts for the immersed boundary. My objective with this was to define a preconditioner that (approximately) solves the same problem as the CG iteration. This means using mean over the whole grid.

Put another way, if we use mean over the whole grid, then the mean pressure produced by the preconditioner would be different than the mean pressure that would minimize the conjugate gradient residual?

It's worth testing though. To implement this properly, we need to know the volume of the immersed region which will require an additional computation.

Δp = mean_p - mean_rhs / δ
Δp = mean_p + mean_rhs / δ
else
Δp = mean_p
end
Expand Down Expand Up @@ -175,22 +183,7 @@ Base.summary(::AsymptoticPoissonPreconditioner) = "AsymptoticPoissonPrecondition
arch = architecture(p)
fill_halo_regions!(r)
launch!(arch, grid, :xyz, _asymptotic_poisson_precondition!, p, grid, r)

δ = preconditioner.regularization
rhs = preconditioner.rhs
mean_p = mean(p)

if !isnothing(δ)
mean_rhs = mean(rhs)
Δp = mean_p - mean_rhs / δ
else
Δp = mean_p
end

grid = p.grid
arch = architecture(grid)
launch!(arch, grid, :xyz, subtract_and_mask!, p, grid, Δp)

regularize_poisson_solution!(p, preconditioner)
return p
end

Expand Down