multi-grid with vector valued space

More
2 years 7 months ago #4312 by walker_tr9
This is related to a recent post of mine.

If I define a FE space like this:  H1(mesh, order=1, dim=2)
then when if I solve my elliptic system with geometric multi-grid, it is fine.  However, if I use "h1amg", then it "silently" crashes; it just exits and says nothing.

However, if I define the space like this:
H1(mesh, order=1) ** 2
then the reverse happens.  I can use "h1amg" no problem, but if I use "multigrid", then it "silently" crashes again.

What is the deal here?
More
2 years 7 months ago #4326 by joachim
The "h1amg" is now (will be online tomorrow)  throwing an exception when applied to H1(dim=2)

The geometric multigrid is working for both versions, H1(dim=2) as well as H1()**2 in my tests. Can you send your example ? 
 
More
2 years 7 months ago #4329 by walker_tr9
Do I need to update my NGsolve?

It will take some time to produce a "standalone" example.
More
2 years 7 months ago #4333 by walker_tr9
Here is a standalone example. It seems to crash when doing the assembly...

import numpy as np

import netgen.gui
from ngsolve import *

from netgen.geom2d import SplineGeometry

# make the mesh
geo = SplineGeometry()
geo.AddRectangle( (0,0), (1,1), bc = "Gamma")
mesh = Mesh(geo.GenerateMesh(maxh=1.0))
for r in range(4):
mesh.Refine()

Draw(mesh)
print(mesh.GetBoundaries())

zero_vec = CoefficientFunction((0.0, 0.0))

fes_V = H1(mesh, order=1) ** 2

a_form = BilinearForm (fes_V, symmetric=False)

# setup weak formulation
qq = fes_V.TrialFunction()
pp = fes_V.TestFunction()

a_form += InnerProduct(grad(qq),grad(pp)) * dx
a_form += InnerProduct(qq,pp) * dx

f_form = LinearForm(fes_V)

f_vec = CoefficientFunction((x, y))
f_form += InnerProduct(f_vec,pp) * dx

# this does NOT work (silently crashes; it doesn't make it to the "Attempt to solve..." statement below)
c = Preconditioner(a_form, "multigrid") # Register 'c' to 'a' BEFORE assembly

# # this does work
# c = Preconditioner(a_form, "h1amg") # Register 'c' to 'a' BEFORE assembly

print("Attempt to assemble...")
a_form.Assemble()
f_form.Assemble()

print("Attempt to solve...")
with TaskManager():
# Conjugate gradient solver
inv = CGSolver(a_form.mat, c.mat, maxsteps=5000) #, precision=1e-12)
# do the solve
vec_data = inv * f_form.vec

print("We solved it!")

input("Press Enter to finish...")
More
2 years 5 months ago #4402 by walker_tr9
Did this ever get fixed? I still get the error with that sample code I gave.
Time to create page: 0.121 seconds