from ngsolve.webgui import Draw

from regpy.util import set_rng_seed
set_rng_seed(42829732344234782348)

Diffusion coefficient problem

We implement the inverse problem to recover the diffusion coefficient given the solution to the PDE. formally that is we consider the problem

\[\nabla a \nabla u = f \text{ in } \Omega\]

\[u = g \text{ on } \partial \Omega.\]

As variational formulation this is given by the bilinear form

\[b(u,v) = \int_\Omega a \nabla(u) \nabla(v) dx\]

and the linear form

\[F(v) = \int_{\partial \Omega} f v ds.\]

As the bilinear form depends linearly on the coefficient \(a\) we may use the provided class SecondOrderEllipticCoefficientPDE.

from regpy.operators.ngsolve import SecondOrderEllipticCoefficientPDE
import ngsolve as ngs

class diffusion(SecondOrderEllipticCoefficientPDE):
    def __init__(self, domain, sol_domain,bdr_val = None,a_bdr_val=None):
        super().__init__(domain, sol_domain, bdr_val=bdr_val,a_bdr_val=a_bdr_val)

    def _bf(self,a,u,v):
        return a*ngs.grad(u)*ngs.grad(v)*ngs.dx

    def _lf(self):
        p = ngs.GridFunction(self.codomain.fes)
        p.Set(-2*ngs.exp(ngs.x+ngs.y))
        lf = ngs.LinearForm(self.codomain.fes)
        lf += p * self.v * ngs.dx
        return lf.Assemble()

For our model problem we use a simple unite_square domain. For the interface to RegPy we then need the NgsVectorSpace from regpy which wraps some functionality around an fes space of NGSolve. Having that in mind we define the mesh and fes spaces for the coefficients and solution spaces and construct the operator by giving it the boundary values. The fes and mesh do not necessarily have to coincide in the setup but we notices that the performance is drastically limited when the meshes are not identical. Note, we use dirichlet boundaries both for the domain and codomain, since we assume these to be known values and thus we can fix them.

from netgen.geom2d import unit_square
from regpy.vecsps import NgsVectorSpace

bdr = "left|top|right|bottom"
mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.1))
fes_domain = ngs.H1(mesh, order=3, dirichlet = bdr)
domain = NgsVectorSpace(fes_domain,bdr = bdr)

bdr = "left|top|right|bottom"
fes_codomain = ngs.H1(mesh, order=3, dirichlet=bdr)
codomain = NgsVectorSpace(fes_codomain, bdr=bdr)

Having the spaces defined and properly wrapped into RegPy vector spaces we can define first the boundary values of the measurements and construct the exact solution and its boundary values. For that we use the convenience method from_ngs of the NgsVectorSpace which converts the CoefficientFunctions to the NgsBaseVector which is a wrapper around the the NGSolve BaseVector and are the vector type used in the definition of operators in RegPy for NGSolve. As you can see you can use definedon to specify a region on which the function should be defined.

bdr_coeff = ngs.sin(ngs.x*4)+2*ngs.y
bdr_val = codomain.from_ngs(bdr_coeff,definedon=codomain.fes.mesh.Boundaries(codomain.bdr))

exact_solution_coeff = 0.5*ngs.exp(-4*(ngs.x-0.5)**2 +4*(ngs.y-0.5)**2)
exact_solution = domain.from_ngs( exact_solution_coeff)

a_bdr_val = domain.from_ngs( exact_solution_coeff, definedon=domain.fes.mesh.Boundaries(domain.bdr))

With the knowledge of the domain, codomain, and boundary values we can now setup our operator.

op = diffusion(
    domain, codomain, bdr_val=bdr_val,a_bdr_val = a_bdr_val
)

To prevent an inverse crime by using the same mesh and fes for the creation of data and reconstruction, we use a refined mesh and new fes to generate data. First we create exact data and add noise which is then interpolated to the spaces for reconstruction.

fine_mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.03))
fine_fes_domain = ngs.H1(fine_mesh, order=4, dirichlet = bdr)
fine_domain = NgsVectorSpace(fine_fes_domain,bdr = bdr)

fine_fes_codomain = ngs.H1(fine_mesh, order=4, dirichlet = bdr)
fine_codomain = NgsVectorSpace(fine_fes_codomain,bdr = bdr)

fine_bdr_val = fine_codomain.from_ngs(bdr_coeff,definedon=fine_codomain.fes.mesh.Boundaries(fine_codomain.bdr))

fine_exact_solution = fine_domain.from_ngs( exact_solution_coeff )
fine_a_bdr_val = fine_domain.from_ngs( exact_solution_coeff, definedon=fine_domain.fes.mesh.Boundaries(fine_domain.bdr) )

fine_op = diffusion(
    fine_domain, fine_codomain, bdr_val=fine_bdr_val,a_bdr_val = fine_a_bdr_val
)

fine_exact_data = fine_op(fine_exact_solution)
fine_noise = 0.01 * fine_codomain.randn()
fine_data = fine_exact_data+fine_noise

# Interpolate into the original space
data = codomain.from_ngs(fine_codomain.to_gf(fine_data))
noise = codomain.from_ngs(fine_codomain.to_gf(fine_noise))

Setting up the inversion

Define a an initial guess by choosing the constant function on the domain and satisfying the apriori known boundary conditions. Then define a regularization setting by choosing appropriate norms on both the domain and codomain. We choose iterative regularized Gauss-Newton method as our regularization scheme. Combined with a stopping rule that is composed of a max iteration and discrepancy principle.

from regpy.solvers import Setting
from regpy.solvers.nonlinear import IrgnmCG
import regpy.stoprules as rules
from regpy.hilbert import Hm0

gfu_constant = ngs.GridFunction(domain.fes)
gfu_constant.Set(1)
gfu_init = ngs.GridFunction(domain.fes)
gfu_init.vec.data = ngs.Projector(domain.fes.FreeDofs(), range=True).Project(gfu_constant.vec) + domain.to_gf(a_bdr_val).vec
init = domain.from_ngs(gfu_init)

setting = Setting(op=op, penalty=Hm0, data_fid=Hm0, data=data)


irgnm = IrgnmCG(setting, regpar = 0.5,init = init)

stoprule = (
        rules.CountIterations(30) +
        rules.Discrepancy(setting.h_codomain.norm,data,setting.h_codomain.norm(noise),tau = 1.3)
)

The inversion

Do the inversion by calling irgnm.run(stoprule).

results = ngs.GridFunction(domain.fes, multidim = 0)
for reco, reco_data in irgnm.until(stoprule):
    results.AddMultiDimComponent(reco.vec.Copy())

2025-12-19 19:09:39,518 WARNING  Setting              :: Setting does not contain any explicit data.
2025-12-19 19:09:39,525 INFO     IrgnmCG.CountIterations :: it. 0>=1000
2025-12-19 19:09:39,528 INFO     IrgnmCG.CountIterations :: it. 1>=1000
2025-12-19 19:09:39,536 WARNING  Setting              :: Setting does not contain any explicit data.
2025-12-19 19:09:39,542 INFO     IrgnmCG.CountIterations :: it. 0>=1000
2025-12-19 19:09:39,545 INFO     IrgnmCG.CountIterations :: it. 1>=1000
2025-12-19 19:09:39,551 WARNING  Setting              :: Setting does not contain any explicit data.
2025-12-19 19:09:39,557 INFO     IrgnmCG.CountIterations :: it. 0>=1000
2025-12-19 19:09:39,560 INFO     IrgnmCG.CountIterations :: it. 1>=1000
2025-12-19 19:09:39,563 INFO     IrgnmCG.CountIterations :: it. 2>=1000
2025-12-19 19:09:39,569 WARNING  Setting              :: Setting does not contain any explicit data.
2025-12-19 19:09:39,575 INFO     IrgnmCG.CountIterations :: it. 0>=1000
2025-12-19 19:09:39,577 INFO     IrgnmCG.CountIterations :: it. 1>=1000
2025-12-19 19:09:39,580 INFO     IrgnmCG.CountIterations :: it. 2>=1000
2025-12-19 19:09:39,585 WARNING  Setting              :: Setting does not contain any explicit data.
2025-12-19 19:09:39,590 INFO     IrgnmCG.CountIterations :: it. 0>=1000
2025-12-19 19:09:39,592 INFO     IrgnmCG.CountIterations :: it. 1>=1000
2025-12-19 19:09:39,594 INFO     IrgnmCG.CountIterations :: it. 2>=1000
2025-12-19 19:09:39,596 INFO     IrgnmCG.CountIterations :: it. 3>=1000
2025-12-19 19:09:39,601 WARNING  Setting              :: Setting does not contain any explicit data.
2025-12-19 19:09:39,606 INFO     IrgnmCG.CountIterations :: it. 0>=1000
2025-12-19 19:09:39,608 INFO     IrgnmCG.CountIterations :: it. 1>=1000
2025-12-19 19:09:39,610 INFO     IrgnmCG.CountIterations :: it. 2>=1000
2025-12-19 19:09:39,612 INFO     IrgnmCG.CountIterations :: it. 3>=1000
2025-12-19 19:09:39,615 INFO     IrgnmCG.CountIterations :: it. 4>=1000
2025-12-19 19:09:39,619 WARNING  Setting              :: Setting does not contain any explicit data.
2025-12-19 19:09:39,624 INFO     IrgnmCG.CountIterations :: it. 0>=1000
2025-12-19 19:09:39,627 INFO     IrgnmCG.CountIterations :: it. 1>=1000
2025-12-19 19:09:39,629 INFO     IrgnmCG.CountIterations :: it. 2>=1000
2025-12-19 19:09:39,631 INFO     IrgnmCG.CountIterations :: it. 3>=1000
2025-12-19 19:09:39,633 INFO     IrgnmCG.CountIterations :: it. 4>=1000
2025-12-19 19:09:39,638 WARNING  Setting              :: Setting does not contain any explicit data.
2025-12-19 19:09:39,643 INFO     IrgnmCG.CountIterations :: it. 0>=1000
2025-12-19 19:09:39,645 INFO     IrgnmCG.CountIterations :: it. 1>=1000
2025-12-19 19:09:39,648 INFO     IrgnmCG.CountIterations :: it. 2>=1000
2025-12-19 19:09:39,650 INFO     IrgnmCG.CountIterations :: it. 3>=1000
2025-12-19 19:09:39,652 INFO     IrgnmCG.CountIterations :: it. 4>=1000
2025-12-19 19:09:39,655 INFO     IrgnmCG.CountIterations :: it. 5>=1000
2025-12-19 19:09:39,659 WARNING  Setting              :: Setting does not contain any explicit data.
2025-12-19 19:09:39,665 INFO     IrgnmCG.CountIterations :: it. 0>=1000
2025-12-19 19:09:39,667 INFO     IrgnmCG.CountIterations :: it. 1>=1000
2025-12-19 19:09:39,669 INFO     IrgnmCG.CountIterations :: it. 2>=1000
2025-12-19 19:09:39,672 INFO     IrgnmCG.CountIterations :: it. 3>=1000
2025-12-19 19:09:39,674 INFO     IrgnmCG.CountIterations :: it. 4>=1000
2025-12-19 19:09:39,676 INFO     IrgnmCG.CountIterations :: it. 5>=1000
2025-12-19 19:09:39,679 INFO     IrgnmCG.CountIterations :: it. 6>=1000
2025-12-19 19:09:39,684 WARNING  Setting              :: Setting does not contain any explicit data.
2025-12-19 19:09:39,691 INFO     IrgnmCG.CountIterations :: it. 0>=1000
2025-12-19 19:09:39,694 INFO     IrgnmCG.CountIterations :: it. 1>=1000
2025-12-19 19:09:39,696 INFO     IrgnmCG.CountIterations :: it. 2>=1000
2025-12-19 19:09:39,699 INFO     IrgnmCG.CountIterations :: it. 3>=1000
2025-12-19 19:09:39,701 INFO     IrgnmCG.CountIterations :: it. 4>=1000
2025-12-19 19:09:39,704 INFO     IrgnmCG.CountIterations :: it. 5>=1000
2025-12-19 19:09:39,706 INFO     IrgnmCG.CountIterations :: it. 6>=1000
2025-12-19 19:09:39,708 INFO     IrgnmCG.CountIterations :: it. 7>=1000
2025-12-19 19:09:39,711 INFO     IrgnmCG.CountIterations :: it. 8>=1000

Draw(exact_solution_coeff, op.domain.fes.mesh)

# Draw reconstructed solution
Draw(results, settings={"Multidim": {"animate": True}})

BaseWebGuiScene

Alternative reconstruction with Landweber

Now if we want to use a different solver we can simply setup another solver with the same setting and a new stop rule. Note that Landweber requires more iterations hence we use a higher maximal iteration count but keep the discrepancy principle to possibly stop early.

from regpy.solvers.nonlinear import Landweber

landweber = Landweber(setting,init = init, data=data)

stoprule = (
        rules.CountIterations(500) +
        rules.Discrepancy(setting.h_codomain.norm,data,setting.h_codomain.norm(noise),tau = 1.3)
)

results = ngs.GridFunction(domain.fes, multidim = 0)
for reco, reco_data in landweber.until(stoprule):
    results.AddMultiDimComponent(reco.vec.Evaluate().Copy())

Draw(exact_solution_coeff, op.domain.fes.mesh, "exact")

# Draw reconstructed solution
Draw(results, settings={"Multidim": {"animate": True, "speed" : 4}})

BaseWebGuiScene