from ngsolve.webgui import Draw

from regpy.util import set_rng_seed
set_rng_seed(42829732344234782348)

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[1], line 1
----> 1 from ngsolve.webgui import Draw
      2
      3 from regpy.util import set_rng_seed
      4 set_rng_seed(42829732344234782348)

File /usr/local/lib/python3.12/site-packages/ngsolve/__init__.py:38
     35     except importlib.metadata.PackageNotFoundError:
     36         pass
---> 38 from .ngslib import __version__, ngstd, bla, la, fem, comp, solve
     40 from netgen import Redraw, TimeFunction
     42 from pyngcore import BitArray, TaskManager, SetNumThreads, PajeTrace, Timers, Timer

ImportError: libmkl_rt.so.2: cannot open shared object file: No such file or directory

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

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[2], line 1
----> 1 from regpy.operators.ngsolve import SecondOrderEllipticCoefficientPDE
      2 import ngsolve as ngs

File /usr/local/lib/python3.12/site-packages/regpy/operators/ngsolve.py:5
      1 r"""PDE forward operators using NGSolve
      2 """
      3 import types
----> 5 import ngsolve as ngs
      6 from pyngcore.pyngcore import BitArray
      7 import numpy as np

File /usr/local/lib/python3.12/site-packages/ngsolve/__init__.py:38
     35     except importlib.metadata.PackageNotFoundError:
     36         pass
---> 38 from .ngslib import __version__, ngstd, bla, la, fem, comp, solve
     40 from netgen import Redraw, TimeFunction
     42 from pyngcore import BitArray, TaskManager, SetNumThreads, PajeTrace, Timers, Timer

ImportError: libmkl_rt.so.2: cannot open shared object file: No such file or directory

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 1
----> 1 class diffusion(SecondOrderEllipticCoefficientPDE):
      2     def __init__(self, domain, sol_domain,bdr_val = None,a_bdr_val=None):
      3         super().__init__(domain, sol_domain, bdr_val=bdr_val,a_bdr_val=a_bdr_val)
      4

NameError: name 'SecondOrderEllipticCoefficientPDE' is not defined

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

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[4], line 2
      1 from netgen.geom2d import unit_square
----> 2 from regpy.vecsps import NgsVectorSpace

ImportError: cannot import name 'NgsVectorSpace' from 'regpy.vecsps' (/usr/local/lib/python3.12/site-packages/regpy/vecsps/__init__.py)

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)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 2
      1 bdr = "left|top|right|bottom"
----> 2 mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.1))
      3 fes_domain = ngs.H1(mesh, order=3, dirichlet = bdr)
      4 domain = NgsVectorSpace(fes_domain,bdr = bdr)
      5

NameError: name 'ngs' is not defined

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 bdr_coeff = ngs.sin(ngs.x*4)+2*ngs.y
      2 bdr_val = codomain.from_ngs(bdr_coeff,definedon=codomain.fes.mesh.Boundaries(codomain.bdr))
      3
      4 exact_solution_coeff = 0.5*ngs.exp(-4*(ngs.x-0.5)**2 +4*(ngs.y-0.5)**2)

NameError: name 'ngs' is not defined

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
)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 1
----> 1 op = diffusion(
      2     domain, codomain, bdr_val=bdr_val,a_bdr_val = a_bdr_val
      3 )

NameError: name 'diffusion' is not defined

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 1
----> 1 fine_mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.03))
      2 fine_fes_domain = ngs.H1(fine_mesh, order=4, dirichlet = bdr)
      3 fine_domain = NgsVectorSpace(fine_fes_domain,bdr = bdr)
      4

NameError: name 'ngs' is not defined

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[9], line 6
      2 from regpy.solvers.nonlinear import IrgnmCG
      3 import regpy.stoprules as rules
      4 from regpy.hilbert import Hm0
      5
----> 6 gfu_constant = ngs.GridFunction(domain.fes)
      7 gfu_constant.Set(1)
      8 gfu_init = ngs.GridFunction(domain.fes)
      9 gfu_init.vec.data = ngs.Projector(domain.fes.FreeDofs(), range=True).Project(gfu_constant.vec) + domain.to_gf(a_bdr_val).vec

NameError: name 'ngs' is not defined

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[10], line 1
----> 1 results = ngs.GridFunction(domain.fes, multidim = 0)
      2 for reco, reco_data in irgnm.until(stoprule):
      3     results.AddMultiDimComponent(reco.vec.Copy())

NameError: name 'ngs' is not defined

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

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[11], line 1
----> 1 Draw(exact_solution_coeff, op.domain.fes.mesh)
      2
      3 # Draw reconstructed solution
      4 Draw(results, settings={"Multidim": {"animate": True}})

NameError: name 'Draw' is not defined

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[13], line 1
----> 1 landweber = Landweber(setting,init = init, data=data)
      2
      3 stoprule = (
      4         rules.CountIterations(500) +

NameError: name 'setting' is not defined

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[14], line 1
----> 1 results = ngs.GridFunction(domain.fes, multidim = 0)
      2 for reco, reco_data in landweber.until(stoprule):
      3     results.AddMultiDimComponent(reco.vec.Evaluate().Copy())

NameError: name 'ngs' is not defined

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

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[15], line 1
----> 1 Draw(exact_solution_coeff, op.domain.fes.mesh, "exact")
      2
      3 # Draw reconstructed solution
      4 Draw(results, settings={"Multidim": {"animate": True, "speed" : 4}})

NameError: name 'Draw' is not defined