import logging

logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s %(levelname)s %(name)-20s :: %(message)s'
)

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

Now we define a mesh and FES spaces for the coefficients and solution spaces and construct the operator by giving it the boundary values. Moreover define the exact solution and take their boundary values as input to the operator.

from netgen.geom2d import unit_square
from regpy.vecsps.ngsolve import NgsSpace

bdr = "left|top|right|bottom"
mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.2))
fes_domain = ngs.H1(mesh, order=6, dirichlet = bdr)
domain = NgsSpace(fes_domain,bdr = bdr)

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

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

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 )
p = ngs.GridFunction(domain.fes)
p.Set(exact_solution_coeff,definedon=domain.fes.mesh.Boundaries(domain.bdr))
a_bdr_val = domain.from_ngs( p )

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

Get data by perturbing the exact data constructed from a finer meshed solution by some random noise.

fine_mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.05))
fine_fes_domain = ngs.H1(fine_mesh, order=6, dirichlet = bdr)
fine_domain = NgsSpace(fine_fes_domain,bdr = bdr)

fine_exact_solution = fine_domain.from_ngs( exact_solution_coeff )
p = ngs.GridFunction(fine_domain.fes)
p.Set(exact_solution_coeff,definedon=fine_domain.fes.mesh.Boundaries(fine_domain.bdr))
fine_a_bdr_val = fine_domain.from_ngs( p )

fine_op = diffusion(
    fine_domain, codomain, bdr_val=bdr_val,a_bdr_val = fine_a_bdr_val
)

exact_data = fine_op(fine_exact_solution)
noise = 0.05 * codomain.randn()
data = exact_data+noise

from ngsolve.webgui import Draw

Draw(codomain.to_ngs(data))

BaseWebGuiScene

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 RegularizationSetting
from regpy.solvers.nonlinear.irgnm 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_ngs(a_bdr_val).vec
init = domain.from_ngs(gfu_init)

Draw(gfu_init)

setting = RegularizationSetting(op=op, penalty=Hm0, data_fid=Hm0)


irgnm = IrgnmCG(setting, data, regpar = 0.01,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).

reco, reco_data = irgnm.run(stoprule)

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

# Draw reconstructed solution
Draw(domain.to_ngs(reco),op.domain.fes.mesh, "reconstruction")

# Draw data space
Draw(codomain.to_ngs(data),op.codomain.fes.mesh, "exact data")
Draw(codomain.to_ngs(reco_data),op.codomain.fes.mesh, "data of reconstruction")

2025-08-13 09:03:40,456 INFO CountIterations      :: iteration = 1 / 30
2025-08-13 09:03:40,457 INFO Discrepancy          :: relative discrepancy = 6.08, tolerance = 1.30
2025-08-13 09:03:40,486 INFO IrgnmCG              :: its.1: alpha=0.006666666666666666, CG its:4
2025-08-13 09:03:40,486 INFO CountIterations      :: iteration = 2 / 30
2025-08-13 09:03:40,487 INFO Discrepancy          :: relative discrepancy = 6.10, tolerance = 1.30
2025-08-13 09:03:40,532 INFO IrgnmCG              :: its.2: alpha=0.004444444444444444, CG its:9
2025-08-13 09:03:40,532 INFO CountIterations      :: iteration = 3 / 30
2025-08-13 09:03:40,533 INFO Discrepancy          :: relative discrepancy = 2.31, tolerance = 1.30
2025-08-13 09:03:40,575 INFO IrgnmCG              :: its.3: alpha=0.0029629629629629624, CG its:8
2025-08-13 09:03:40,575 INFO CountIterations      :: iteration = 4 / 30
2025-08-13 09:03:40,576 INFO Discrepancy          :: relative discrepancy = 0.96, tolerance = 1.30
2025-08-13 09:03:40,576 INFO CombineRules         :: Rule Discrepancy(noiselevel=0.15671015520455536, tau=1.3) triggered.
2025-08-13 09:03:40,576 INFO IrgnmCG              :: Solver converged after 3 iteration.

BaseWebGuiScene