Composing a solution from subregions

General
1

Hello everyone; I would like to be able to split my simulation region into parts; solve those partially, combine those to get the solution.
The reason is that the problem is composed of repeating structures and I hope to increase the calculation speed.
Can anyone point to an example; I could not find any.
Thanks in advance

2

Have a look into this:

https://jschoeberl.github.io/iFEM/subspacecorrection/minimaldd.html

and also section 72 ff

Hope this gives some ideas.

Joachim

3

Thank you very much; I was not aware of that resource. That should at least keep me busy for some time.

4

I had a look at the given example. Since I am using a more complex geometry I tried to extend the example. As a first step I wanted to combine sets of 3 rectangles to form 3 regions.

    regions = []
    key = lambda x: x[-1]
    for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key):
        reg = reduce(operator.add, (mesh.Materials(m) for m in g))
        regions.append(reg)

    fes = FESpace([Compress(H1(mesh, definedon=dom, dirichlet="bot")) for dom in regions])

    a = BilinearForm(fes)
    f = LinearForm(fes)

    for (ui,vi) in zip(u,v):
        a += grad(ui)*grad(vi)*dx + 1*ui*vi*dx
        f += y*x*vi*dx

    gfu = GridFunction(fes)
    gfu.vec.data = a.mat.Inverse(inverse="sparsecholesky", freedofs=fes.FreeDofs())*f.vec
    gftot = CF(list(gfu.components))

This clearly does not work; (simply check by actually composing the the geometry to 3 larger rectangles).
What am I missing here.
Best to everyone

5

if you post a complete example, your chance for help increases
Joachim

6

My apologies; posting a minimal working example should have been a no-brainer.
Here it is and thanks for your quick reply


from ngsolve import *
from netgen.geom2d import CSG2d, Rectangle
from itertools import groupby
from functools import reduce
import operator


geo = CSG2d()

mx, my = 3,3
for i in range(mx): 
    for j in range(my):
        rect = Rectangle(pmin=(i/mx,j/my), \
                        pmax=((i+1)/mx,(j+1)/my), \
                        mat='mat'+str(i)+str(j), \
                        bc = 'default', \
                        bottom = 'bot' if j == 0 else None)
        geo.Add(rect)
        
mesh = Mesh(geo.GenerateMesh(maxh=0.04))

key = lambda x: x[-1]

mats = ['|'.join(m for m in g) for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key)]
regions = [reg:=mesh.Materials(mat) for mat in mats]

# equivalent result using this construction of the regions
# mats = [list(g) for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key)]
# regions = []
# for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key):
#     reg = reduce(operator.add, (mesh.Materials(m) for m in g))
#     regions.append(reg)

fes = FESpace([Compress(H1(mesh, definedon=dom, dirichlet="bot")) for dom in regions])

print ("ndof =", fes.ndof)

u, v = fes.TnT()

a = BilinearForm(fes)
f = LinearForm(fes)

for (ui,vi) in zip(u,v):
    a += grad(ui)*grad(vi)*dx + 1*ui*vi*dx
    f += y*x*vi*dx
        
a.Assemble()
f.Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(inverse="sparsecholesky", freedofs=fes.FreeDofs())*f.vec
gftot = CF(list(gfu.components))
Draw(gftot)
7

Using the method mentioned in your initial response, I can get it to work.
Still, I am confused what went wrong with the other approach


from ngsolve import *
from netgen.geom2d import CSG2d, Rectangle
from itertools import groupby
from functools import reduce
import operator
from ngsolve.krylovspace import CGSolver


geo = CSG2d()

mx, my = 3,3
for i in range(mx): 
    for j in range(my):
        rect = Rectangle(pmin=(i/mx,j/my), \
                        pmax=((i+1)/mx,(j+1)/my), \
                        mat='mat'+str(i)+str(j), \
                        bc = 'default', \
                        bottom = 'bot' if j == 0 else None)
        geo.Add(rect)
        
mesh = Mesh(geo.GenerateMesh(maxh=0.04))

key = lambda x: x[-1]
mats = [list(g) for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key)]

regions = []
for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key):
    reg = reduce(operator.add, (mesh.Materials(m) for m in g))
    regions.append(reg)


fes = H1(mesh, order=1, dirichlet='bot')

print ("ndof =", fes.ndof)

u, v = fes.TnT()

a = BilinearForm(fes)
f = LinearForm(fes)
f += x*v*dx
a += grad(u)*grad(v)*dx
f.Assemble()
a.Assemble()

invs = []
for dom in regions:
    domaindofs = fes.GetDofs(dom) & fes.FreeDofs()
    print ("num dofs:", domaindofs.NumSet())
    invs.append(a.mat.Inverse(domaindofs))

invs = reduce(operator.add, invs)

inv = CGSolver(mat=a.mat, pre=invs)
print ("Iterations: ", len(inv.errors))
gfu = GridFunction(fes)
gfu.vec.data = inv * f.vec
gftot = CF(gfu)
Draw(gftot)
8

Some apologies at this point: I got it to work; since I use custom drawing functions, drawing the result was wrong. it is still a bit surprising but with the following replacement I got it to work:

gftot = reduce(operator.add, (CF(comp) for comp in gfu.components))
# instead of 
gftot = CF(list(gfu.components))
9

CF (list) is domain wise, CF(tuple) vectorial. it has some legacy reasons and is a quite common pitfall (as seen in documentation)