Add example solving -div(a∇u)=f using DivAgradOperator with general inclusions#1228
Draft
Add example solving -div(a∇u)=f using DivAgradOperator with general inclusions#1228
Conversation
Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Copilot
AI
changed the title
[WIP] Implement a simple div(a grad u) = f solver
Add example solving -div(a∇u)=f using DivAgradOperator
Jan 14, 2026
Member
|
@copilot did you actually run the example to make sure it works? |
Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Contributor
Author
Yes, I've now tested it. The example runs successfully and solves the linear problem. The CG solver converges in 15 iterations with a final residual of 1.68e-7. I also fixed a bug in commit 717b3a4 where Test output: The only part that doesn't run is the final plotting (which requires the Plots package), but all the computational code works correctly. Warning Firewall rules blocked me from connecting to one or more addresses (expand for details)I tried to connect to the following addresses, but was blocked by firewall rules:
If you need me to access, download, or install something from one of these locations, you can either:
|
Member
|
@copilot can you run the code, exporting the final plot to a pdf, look at the pdf, see if it looks reasonable physically and iterate from there? |
Contributor
Author
I attempted to run the code with Plots to generate a PDF, but the installation of Plots and its dependencies (over 200 packages) exceeded the timeout limit. However, I've verified the code works correctly:
Without visualization, I can't definitively verify the physical reasonableness of the solution, but the mathematical structure and convergence behavior indicate the implementation is correct. If you'd like to verify the visualization yourself, you can run the example with Plots installed in your local environment. Warning Firewall rules blocked me from connecting to one or more addresses (expand for details)I tried to connect to the following addresses, but was blocked by firewall rules:
If you need me to access, download, or install something from one of these locations, you can either:
|
Member
|
@copilot hm this is annoying. Any idea why this happens? Any way to use a pre-installed environment with all dependencies installed so you don't have to re-install everything each time? |
Contributor
Author
The timeout happens because each execution starts in a fresh sandboxed environment without any pre-installed packages. When Julia's Pkg.add("Plots") runs, it needs to:
Unfortunately, I don't have access to a persistent Julia environment or the ability to configure pre-installed package depots. Each run starts fresh. For verification, I'd suggest:
The core solver implementation is working correctly - the issue is purely with visualization/validation. |
Member
antoine-levitt
left a comment
antoine-levitt
left a comment
There was a problem hiding this comment.
Very nice, just a couple of minor points
| in the div-grad problem. The inclusion has a constant value inside a ball of given radius. | ||
| """ | ||
| struct ElementSphericalInclusion <: DFTK.Element | ||
| a_value::Float64 # Value of the coefficient modification in the inclusion |
Member
There was a problem hiding this comment.
In julia we typically don't specialize to Float64, take the floating point type as a type parameter and type a_value and radius as ::T
| end | ||
|
|
||
| # Default constructor | ||
| ElementSphericalInclusion(; a_value=1.0, radius=0.5) = ElementSphericalInclusion(a_value, radius) |
Member
There was a problem hiding this comment.
don't have a default constructor
| # FT[χ_R](p) = 4π R³/3 * 3(sin(p*R) - p*R*cos(p*R))/(p*R)³ | ||
| # However, this is for the characteristic function. For our coefficient a(x), | ||
| # we want the value a_value in the ball. | ||
| function DFTK.local_potential_fourier(el::ElementSphericalInclusion, p::Real) |
Member
There was a problem hiding this comment.
same here, don't type p, take it arbitrary. And when using pi, use T(pi) to avoid accidentally converting to float64
| AtomicDivAGrad: Construct the coefficient field a(x) from atomic positions. | ||
| """ | ||
| struct AtomicDivAGrad | ||
| background_value::Float64 # Background uniform value of a(x) |
Member
There was a problem hiding this comment.
same here, use a parametric type
| background_value::Float64 # Background uniform value of a(x) | ||
| end | ||
|
|
||
| AtomicDivAGrad(; background_value=1.0) = AtomicDivAGrad(background_value) |
Member
There was a problem hiding this comment.
don't give a default value
| function proj(x) | ||
| # Remove the zero Fourier mode (constant component) | ||
| x_copy = copy(x) | ||
| if !isnothing(G_zero_idx) |
Member
There was a problem hiding this comment.
just @assert previously it's not nothing (otherwise we've got a bug and should error as soon as possible)
| # Define spherical inclusions at specific positions | ||
| # These act as "atoms" that modify the coefficient a(x) | ||
| inclusion = ElementSphericalInclusion(a_value=2.0, radius=1.0) | ||
| positions = [[0.25, 0.25, 0.0], [0.75, 0.75, 0.0]] |
Member
There was a problem hiding this comment.
this is too tame and too symmetric. Use radius=2, a_value=4, and more random-looking positions
| f_values = zeros(Float64, basis.fft_size...) | ||
| for (i, r) in enumerate(r_vectors) | ||
| x, y = r[1], r[2] | ||
| f_values[i] = sin(2π * x / a) * sin(2π * y / a) |
Member
There was a problem hiding this comment.
this is too symmetric, use something a bit more interesting (but still simple)
|
|
||
| p = plot(p1, p2, p3, layout=(1, 3), size=(1200, 400)) | ||
|
|
||
| # Save to PDF if requested |
| @@ -0,0 +1,17 @@ | |||
| # Script to run divAgrad_solver example and export plot to PDF | |||
Member
|
@copilot adress the code review |
… better test case Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Contributor
Author
Addressed all code review comments in commit 79f3dd2:
|
Member
|
@copilot Also can you add an option to use as preconditioner 1/sqrt(a) (-Delta)^-1 1/sqrt(a)? (as a custom type that implements ldiv!, you'll need to use ffts and iffts, don't bother too much about reusing computations, this is only for playing with it). Also you're good at math, no? Add in the comments in the file the detailed proof that the preconditioned operator is of the form 1 + compact. |
…vel, correct isposdef flag Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Member
|
@copilot the previous comment is adressed to you, I forgot to tag you and then edited the message but it didn't work |
…ompact form Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Contributor
Author
Added in commit 5e528a8: New DivAGradPreconditioner: Implements M = 1/√a (-Δ)^(-1) 1/√a using FFTs for the multiplication steps and division by |G|² for the Laplacian inverse. Mathematical proof (in comments): Detailed derivation showing that M⁻¹H = a₀I + K where K is compact:
Usage: Example now compares both preconditioners and reports iteration counts.
Addressed - see response to previous comment. Added the 1/√a (-Δ)^(-1) 1/√a preconditioner with full mathematical proof in commit 5e528a8. |
Copilot
AI
changed the title
Add example solving -div(a∇u)=f using DivAgradOperator
Add example solving -div(a∇u)=f using DivAgradOperator with multiple preconditioners
Jan 14, 2026
…f of I+compact form" This reverts commit 5e528a8.
Member
|
The combined preconditioner doesn't appear better, it's a bit weird, reverting for now. Good to go on my end. Ping @mfherbst |
Member
|
@copilot Test it using rectangular inclusions, and compare the results using DivAGradFromFourier and DivAGradFromReal. Also take f = da/dx (derivative wrt to the first direction), calculated exactly in the Fourier domain. To visualize the result, add two additional plots of du/dx and du/dy. |
…urier Co-authored-by: antoine-levitt <88065+antoine-levitt@users.noreply.github.com>
Contributor
Author
Implemented all requested changes in commit 5da7bde: Major refactoring:
Real-space construction:
Testing:
Visualization:
|
Copilot
AI
changed the title
Add example solving -div(a∇u)=f using DivAgradOperator with multiple preconditioners
Add example solving -div(a∇u)=f using DivAgradOperator with general inclusions
Jan 14, 2026
Member
|
That one was not as good, claude produced relatively verbose (copy-paste) code, messed up on cartesian vs relative coordinates, and on kpoint basis vs basis (to be fair, humans often do also...) However telling it to clean up a piece of code (in particular I told it to be more concise) is absolutely amazing. |
Member
|
A bit verbose, but not all that bad. @antoine-levitt What's your plan with this ? Just curiosity and playing around, convert it into something that finally makes it into DFTK, something else ? |
Member
|
This was mostly to play around with AI and to test something for some friends but I think it's a generally usefulthing to have in DFTK, eg to plot photonic band structure (although we probably need to merge in multicomponents to do actual photonics). I think we can just add a divagradfromreal and divagradfromfourier terms, following the local potential stuff, no? |
Member
|
Also the Ax=b example I think is generically useful, we could have a "DFTK as a generic PDE solver" example. I thought of having an option to turn off the manifold stuff and turning it into a generic "optimize a functional defined on several periodic functions" machine but that breaks some of the core assumptions (in particular that the gradient of the energy is of the form H psi) so it's probably better not to do that. I say we (and when I say we I mean Claude) split the file into a part that goes into a divagrad.jl term, and an example that has the Ax=b solver, OK? |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Plan for implementing div(a grad u)=f solver
Implementation Details
TermDivAGrad - Term using DivAgradOperator, properly scaled to account for the -½ factor. No longer limited to spherical symmetry or atom-like structures.
DivAGradFromReal - Constructs coefficient field a(x) from a real-space inclusion function that accepts Vec3 position vectors. Handles periodic boundary conditions by summing over neighboring cells.
DivAGradFromFourier - Constructs coefficient field a(x) from a Fourier-space inclusion function that accepts Vec3 wavevectors. Uses structure factors for positioned inclusions.
Helper Functions:
build_coefficient_field_real: Builds a(x) in real space with periodic image handling (±1 cells in each direction)build_coefficient_field_fourier: Builds a(x) in Fourier space with structure factorsPreconditioner:
solve_linear_problem - CG-based solver with proper G=0 mode projection and DFTK helper functions
Example Usage
The example demonstrates solving -div(a(x)∇u(x)) = f(x) with rectangular inclusions using both real-space and Fourier-space construction methods. It compares the two approaches and visualizes:
All code is self-contained in
examples/divAgrad_solver.jl.Original prompt
💬 We'd love your input! Share your thoughts on Copilot coding agent in our 2 minute survey.