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)
)

2026-01-29 15:10:02,925 WARNING  Discrepancy          :: Initialization of Discrepancy with three positional arguments (norm, data, noise_level) deprecated. Use one positional argument (noiselevel) and specifiy norm and data via a keyword argument setting!

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())

2026-01-29 15:10:02,947 INFO     CombineRules         :: it. 0>=1000 | ((rel X:--(x=0)! | rel Y:--(y=0) ) & kappa:1.0e+00) | it. 0>=1000
2026-01-29 15:10:02,950 INFO     CombineRules         :: it. 1>=1000 | ((rel X:8.8e-02>=2.5e-01  | rel Y:1.5e-01>=2.5e-01 [Rule RelTolXStop(0.3333333333333333) triggered.Rule RelTolYStop(0.3333333333333333) triggered.]) & kappa:1.0e+00[All rules triggered.]) | it. 1>=1000
2026-01-29 15:10:02,956 INFO     CombineRules         :: it. 0>=30 | discr./noiselevel = 3.37< 1.30
2026-01-29 15:10:02,968 INFO     CombineRules         :: it. 0>=1000 | ((rel X:--(x=0)! | rel Y:4.9e+14>=2.5e-01 ) & kappa:1.0e+00) | it. 0>=1000
2026-01-29 15:10:02,971 INFO     CombineRules         :: it. 1>=1000 | ((rel X:2.5e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.0e+00[All rules triggered.]) | it. 1>=1000
2026-01-29 15:10:02,975 INFO     CombineRules         :: it. 1>=30 | discr./noiselevel = 3.01< 1.30
2026-01-29 15:10:02,986 INFO     CombineRules         :: it. 0>=1000 | ((rel X:--(x=0)! | rel Y:1.0e+15>=2.5e-01 ) & kappa:1.0e+00) | it. 0>=1000
2026-01-29 15:10:02,989 INFO     CombineRules         :: it. 1>=1000 | ((rel X:5.0e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.1e+00) | it. 1>=1000
2026-01-29 15:10:02,992 INFO     CombineRules         :: it. 2>=1000 | ((rel X:9.6e-02>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.0e+00[All rules triggered.]) | it. 2>=1000
2026-01-29 15:10:02,998 INFO     CombineRules         :: it. 2>=30 | discr./noiselevel = 2.67< 1.30
2026-01-29 15:10:03,010 INFO     CombineRules         :: it. 0>=1000 | ((rel X:--(x=0)! | rel Y:5.6e+15>=2.5e-01 ) & kappa:1.0e+00) | it. 0>=1000
2026-01-29 15:10:03,015 INFO     CombineRules         :: it. 1>=1000 | ((rel X:8.8e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.1e+00) | it. 1>=1000
2026-01-29 15:10:03,018 INFO     CombineRules         :: it. 2>=1000 | ((rel X:1.9e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.1e+00[All rules triggered.]) | it. 2>=1000
2026-01-29 15:10:03,028 INFO     CombineRules         :: it. 3>=30 | discr./noiselevel = 2.35< 1.30
2026-01-29 15:10:03,035 INFO     CombineRules         :: it. 0>=1000 | ((rel X:--(x=0)! | rel Y:1.1e+16>=2.5e-01 ) & kappa:1.0e+00) | it. 0>=1000
2026-01-29 15:10:03,039 INFO     CombineRules         :: it. 1>=1000 | ((rel X:1.8e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.2e+00) | it. 1>=1000
2026-01-29 15:10:03,042 INFO     CombineRules         :: it. 2>=1000 | ((rel X:3.8e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.1e+00) | it. 2>=1000
2026-01-29 15:10:03,044 INFO     CombineRules         :: it. 3>=1000 | ((rel X:1.5e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.2e+00[All rules triggered.]) | it. 3>=1000
2026-01-29 15:10:03,056 INFO     CombineRules         :: it. 4>=30 | discr./noiselevel = 2.07< 1.30
2026-01-29 15:10:03,064 INFO     CombineRules         :: it. 0>=1000 | ((rel X:--(x=0)! | rel Y:1.7e+16>=2.5e-01 ) & kappa:1.0e+00) | it. 0>=1000
2026-01-29 15:10:03,067 INFO     CombineRules         :: it. 1>=1000 | ((rel X:3.1e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.2e+00) | it. 1>=1000
2026-01-29 15:10:03,069 INFO     CombineRules         :: it. 2>=1000 | ((rel X:8.3e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.3e+00) | it. 2>=1000
2026-01-29 15:10:03,071 INFO     CombineRules         :: it. 3>=1000 | ((rel X:3.3e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.3e+00) | it. 3>=1000
2026-01-29 15:10:03,073 INFO     CombineRules         :: it. 4>=1000 | ((rel X:1.2e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.2e+00[All rules triggered.]) | it. 4>=1000
2026-01-29 15:10:03,084 INFO     CombineRules         :: it. 5>=30 | discr./noiselevel = 1.84< 1.30
2026-01-29 15:10:03,093 INFO     CombineRules         :: it. 0>=1000 | ((rel X:--(x=0)! | rel Y:2.8e+16>=2.5e-01 ) & kappa:1.0e+00) | it. 0>=1000
2026-01-29 15:10:03,095 INFO     CombineRules         :: it. 1>=1000 | ((rel X:5.3e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.4e+00) | it. 1>=1000
2026-01-29 15:10:03,097 INFO     CombineRules         :: it. 2>=1000 | ((rel X:1.3e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.4e+00) | it. 2>=1000
2026-01-29 15:10:03,100 INFO     CombineRules         :: it. 3>=1000 | ((rel X:6.2e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.5e+00) | it. 3>=1000
2026-01-29 15:10:03,101 INFO     CombineRules         :: it. 4>=1000 | ((rel X:2.7e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.3e+00) | it. 4>=1000
2026-01-29 15:10:03,103 INFO     CombineRules         :: it. 5>=1000 | ((rel X:1.0e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.2e+00[All rules triggered.]) | it. 5>=1000
2026-01-29 15:10:03,111 INFO     CombineRules         :: it. 6>=30 | discr./noiselevel = 1.64< 1.30
2026-01-29 15:10:03,118 INFO     CombineRules         :: it. 0>=1000 | ((rel X:--(x=0)! | rel Y:8.9e+16>=2.5e-01 ) & kappa:1.0e+00) | it. 0>=1000
2026-01-29 15:10:03,120 INFO     CombineRules         :: it. 1>=1000 | ((rel X:8.5e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.5e+00) | it. 1>=1000
2026-01-29 15:10:03,122 INFO     CombineRules         :: it. 2>=1000 | ((rel X:2.0e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.6e+00) | it. 2>=1000
2026-01-29 15:10:03,124 INFO     CombineRules         :: it. 3>=1000 | ((rel X:1.0e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.7e+00) | it. 3>=1000
2026-01-29 15:10:03,126 INFO     CombineRules         :: it. 4>=1000 | ((rel X:5.0e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.5e+00) | it. 4>=1000
2026-01-29 15:10:03,129 INFO     CombineRules         :: it. 5>=1000 | ((rel X:2.3e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.3e+00[All rules triggered.]) | it. 5>=1000
2026-01-29 15:10:03,135 INFO     CombineRules         :: it. 7>=30 | discr./noiselevel = 1.46< 1.30
2026-01-29 15:10:03,143 INFO     CombineRules         :: it. 0>=1000 | ((rel X:--(x=0)! | rel Y:5.1e+16>=2.5e-01 ) & kappa:1.0e+00) | it. 0>=1000
2026-01-29 15:10:03,145 INFO     CombineRules         :: it. 1>=1000 | ((rel X:1.3e+01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.7e+00) | it. 1>=1000
2026-01-29 15:10:03,147 INFO     CombineRules         :: it. 2>=1000 | ((rel X:3.0e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.8e+00) | it. 2>=1000
2026-01-29 15:10:03,149 INFO     CombineRules         :: it. 3>=1000 | ((rel X:1.6e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:2.1e+00) | it. 3>=1000
2026-01-29 15:10:03,151 INFO     CombineRules         :: it. 4>=1000 | ((rel X:8.2e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.7e+00) | it. 4>=1000
2026-01-29 15:10:03,153 INFO     CombineRules         :: it. 5>=1000 | ((rel X:4.2e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.5e+00) | it. 5>=1000
2026-01-29 15:10:03,154 INFO     CombineRules         :: it. 6>=1000 | ((rel X:2.2e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.4e+00[All rules triggered.]) | it. 6>=1000
2026-01-29 15:10:03,158 INFO     CombineRules         :: it. 8>=30 | discr./noiselevel = 1.31< 1.30
2026-01-29 15:10:03,167 INFO     CombineRules         :: it. 0>=1000 | ((rel X:--(x=0)! | rel Y:7.1e+16>=2.5e-01 ) & kappa:1.0e+00) | it. 0>=1000
2026-01-29 15:10:03,172 INFO     CombineRules         :: it. 1>=1000 | ((rel X:2.0e+01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.9e+00) | it. 1>=1000
2026-01-29 15:10:03,175 INFO     CombineRules         :: it. 2>=1000 | ((rel X:4.5e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:2.0e+00) | it. 2>=1000
2026-01-29 15:10:03,177 INFO     CombineRules         :: it. 3>=1000 | ((rel X:2.3e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:2.2e+00) | it. 3>=1000
2026-01-29 15:10:03,178 INFO     CombineRules         :: it. 4>=1000 | ((rel X:1.3e+00>=2.5e-01  | rel Y:--(y=0) ) & kappa:2.2e+00) | it. 4>=1000
2026-01-29 15:10:03,180 INFO     CombineRules         :: it. 5>=1000 | ((rel X:7.1e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.7e+00) | it. 5>=1000
2026-01-29 15:10:03,182 INFO     CombineRules         :: it. 6>=1000 | ((rel X:4.3e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.8e+00) | it. 6>=1000
2026-01-29 15:10:03,184 INFO     CombineRules         :: it. 7>=1000 | ((rel X:2.9e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.9e+00) | it. 7>=1000
2026-01-29 15:10:03,185 INFO     CombineRules         :: it. 8>=1000 | ((rel X:1.7e-01>=2.5e-01  | rel Y:--(y=0) ) & kappa:1.5e+00[All rules triggered.]) | it. 8>=1000
2026-01-29 15:10:03,189 INFO     CombineRules         :: it. 9>=30 | discr./noiselevel = 1.16< 1.30
Rule Discrepancy(noiselevel=0.23242976722781797, tau=1.3) triggered.

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)
)

2026-01-29 15:10:03,291 WARNING  Setting              :: Overwriting existing initial guess in setting!
2026-01-29 15:10:03,299 WARNING  Discrepancy          :: Initialization of Discrepancy with three positional arguments (norm, data, noise_level) deprecated. Use one positional argument (noiselevel) and specifiy norm and data via a keyword argument setting!

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

2026-01-29 15:10:03,595 INFO     CombineRules         :: it. 0>=500 | discr./noiselevel = 2.52< 1.30
2026-01-29 15:10:03,618 INFO     CombineRules         :: it. 1>=500 | discr./noiselevel = 2.31< 1.30
2026-01-29 15:10:03,631 INFO     CombineRules         :: it. 2>=500 | discr./noiselevel = 1.97< 1.30
2026-01-29 15:10:03,644 INFO     CombineRules         :: it. 3>=500 | discr./noiselevel = 1.86< 1.30
2026-01-29 15:10:03,657 INFO     CombineRules         :: it. 4>=500 | discr./noiselevel = 1.76< 1.30
2026-01-29 15:10:03,674 INFO     CombineRules         :: it. 5>=500 | discr./noiselevel = 1.69< 1.30
2026-01-29 15:10:03,695 INFO     CombineRules         :: it. 6>=500 | discr./noiselevel = 1.63< 1.30
2026-01-29 15:10:03,719 INFO     CombineRules         :: it. 7>=500 | discr./noiselevel = 1.58< 1.30
2026-01-29 15:10:03,741 INFO     CombineRules         :: it. 8>=500 | discr./noiselevel = 1.53< 1.30
2026-01-29 15:10:03,755 INFO     CombineRules         :: it. 9>=500 | discr./noiselevel = 1.49< 1.30
2026-01-29 15:10:03,771 INFO     CombineRules         :: it. 10>=500 | discr./noiselevel = 1.46< 1.30
2026-01-29 15:10:03,787 INFO     CombineRules         :: it. 11>=500 | discr./noiselevel = 1.42< 1.30
2026-01-29 15:10:03,801 INFO     CombineRules         :: it. 12>=500 | discr./noiselevel = 1.40< 1.30
2026-01-29 15:10:03,815 INFO     CombineRules         :: it. 13>=500 | discr./noiselevel = 1.37< 1.30
2026-01-29 15:10:03,834 INFO     CombineRules         :: it. 14>=500 | discr./noiselevel = 1.34< 1.30
2026-01-29 15:10:03,847 INFO     CombineRules         :: it. 15>=500 | discr./noiselevel = 1.32< 1.30
2026-01-29 15:10:03,863 INFO     CombineRules         :: it. 16>=500 | discr./noiselevel = 1.30< 1.30
Rule Discrepancy(noiselevel=0.23242976722781797, tau=1.3) triggered.

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

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

BaseWebGuiScene