Adaptive mesh refinement for time dependant problems

NGSolve
,
1

Hi there

I want to implement adaptive mesh refinement for a time dependant problem, say for this example 3.1 Parabolic model problem — NGS-Py 6.2.2404-2-ge6a630133 documentation

I tried to do the mesh.refine() for each time step, it works for the first time step, but for next time steps ndof has already reached the 100000 limit and, as a rest, it doesn’t enter the refinement loop again. Is there any way to reset the mesh to initial mesh after each time step? Or any other suggestions?

Regards
Hamed

2

unrefinement is not supported in NGSolve. It’s on the todo-list.

You can copy the initial mesh before refinement, and start a new refinement loop. You can use the GridFunction on the old mesh as the function of the previous time-step.

Joachim

3

Hi Joachim

Thanks for your reply. Sorry for my dumb question. If gfu is the solution of the previous step on the refined mesh, is this the correct way of using gfu on the old mesh as the new initial condition?

old_mesh=Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.4))
new_mesh = old_mesh.ngmesh.Copy()
new_fes = H1(new_mesh, order=3, autoupdate=True)
new_initial_cond = GridFunction(new_fes, autoupdate=True)
new_initial_cond.Set(gfu)

Regards
Hamed

4

yes, and as you keep refining new_mesh, you project the old gfu into the new space, on every level.

5

Hi Joachim

I am copying my code here. When I run this code, kernel dies. If I remove while fes.ndof < 100000, it works, but obviously doesn’t refine enough. It looks like I am getting out of memory. Any suggestion on what I am doing wrong?

heat3.ipynb (4.5 KB)

from ngsolve import *
from netgen.geom2d import SplineGeometry
from ngsolve.webgui import Draw
from netgen.occ import *
from netgen.webgui import Draw as DrawGeo
shape = Rectangle(2,2).Face().Move((-1,-1,0))
shape.edges.Min(X).name="left"
shape.edges.Max(X).name="right"
shape.edges.Min(Y).name="bottom"
shape.edges.Max(Y).name="top"

time = 0.0
dt = 0.001

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.25))



def Solve(gfu, fes, u, v):
    b = CoefficientFunction((2*y*(1-x*x),-2*x*(1-y*y)))
    a = BilinearForm(fes, symmetric=False)
    a += 0.01*grad(u)*grad(v)*dx + b*grad(u)*v*dx
    a.Assemble()
    m = BilinearForm(fes, symmetric=False)
    m += u*v*dx
    m.Assemble()
    mstar = m.mat.CreateMatrix()
    mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()
    invmstar = mstar.Inverse(freedofs=fes.FreeDofs())
    gaussp = exp(-6*((x+0.5)*(x+0.5)+y*y))-exp(-6*((x-0.5)*(x-0.5)+y*y))
    f = LinearForm(fes)
    f += gaussp*v*dx
    f.Assemble()
    res = dt * f.vec - dt * a.mat * gfu.vec
    gfu.vec.data += invmstar * res
    return gfu

space_flux = HDiv(mesh, order=2, autoupdate=True)
gf_flux = GridFunction(space_flux, "flux", autoupdate=True)

def CalcError(gfu, fes, new_mesh):
    flux = 0.01 * grad(gfu)
    # interpolate finite element flux into H(div) space:
    gf_flux.Set (flux)

    # Gradient-recovery error estimator
    err = 1/0.01*(flux-gf_flux)*(flux-gf_flux)
    elerr = Integrate (err, new_mesh, VOL, element_wise=True)

    maxerr = max(elerr)
    #print ("maxerr = ", maxerr)

    for el in new_mesh.Elements():
        new_mesh.SetRefinementFlag(el, elerr[el.nr] > 0.25*maxerr)
    return new_mesh


def TimeStepping(initial_cond = None, t0 = 0, tend = 2,
                 nsamples = 10):
    if initial_cond:
        fes = H1(Mesh(mesh.ngmesh.Copy()), order=3, dirichlet="bottom|right|left|top", autoupdate=True)
        gfu = GridFunction(fes, autoupdate=True)
        gfu.Set(initial_cond)
    cnt = 0; time = t0
    sample_int = int(floor(tend / dt / nsamples)+1)
    gfut = GridFunction(gfu.space,multidim=0)
    gfut.AddMultiDimComponent(gfu.vec)
    while time < tend - 0.5 * dt:
        new_mesh = Mesh(mesh.ngmesh.Copy())
        fes = H1(new_mesh, order=3, dirichlet="bottom|right|left|top", autoupdate=True)
        u,v = fes.TnT()
        with TaskManager():
            while fes.ndof < 100000:
                Solve(gfu, fes, u, v)
                CalcError(gfu, fes, new_mesh)
                new_mesh.Refine()
        #Solve(gfu, fes, u, v)


        print("\r",time,end="")

        if cnt % sample_int == 0:
            gfut.AddMultiDimComponent(gfu.vec)
        cnt += 1; time = cnt * dt
    return gfut, new_mesh

gfut, new_mesh = TimeStepping((1-y*y)*x)

Draw(gfut, new_mesh, interpolate_multidim=True, animate=True)
6

Hi

I made some changes. It is working smoothly now, but the result is not what I expect. Could someone please check what I am doing wrong?

from ngsolve import *
from netgen.geom2d import SplineGeometry
from ngsolve.webgui import Draw
from netgen.occ import *
from netgen.webgui import Draw as DrawGeo
shape = Rectangle(2,2).Face().Move((-1,-1,0))
shape.edges.Min(X).name="left"
shape.edges.Max(X).name="right"
shape.edges.Min(Y).name="bottom"
shape.edges.Max(Y).name="top"

time = 0.0
dt = 0.1

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.25))
fes_old = H1(mesh, order=3, dirichlet="bottom|right|left|top")
gfu_main = GridFunction(fes_old)
new_mesh = Mesh(mesh.ngmesh.Copy())
fes = H1(new_mesh, order=3, dirichlet="bottom|right|left|top", autoupdate=True)
u,v = fes.TnT()



gfu = GridFunction(fes, autoupdate=True)
gfu.Set((1-y*y)*x)



def Solve(gfu, fes, u, v):
    b = CoefficientFunction((2*y*(1-x*x),-2*x*(1-y*y)))
    gaussp = exp(-6*((x+0.5)*(x+0.5)+y*y))-exp(-6*((x-0.5)*(x-0.5)+y*y))
    a = BilinearForm(fes, symmetric=False)
    a += 0.01*grad(u)*grad(v)*dx + b*grad(u)*v*dx
    a.Assemble()
    m = BilinearForm(fes, symmetric=False)
    m += u*v*dx
    m.Assemble()
    mstar = m.mat.CreateMatrix()
    mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()
    invmstar = mstar.Inverse(freedofs=fes.FreeDofs())
    f = LinearForm(fes)
    f += gaussp*v*dx
    f.Assemble()
    res = dt * f.vec - dt * a.mat * gfu.vec
    gfu.vec.data = invmstar * res
    return gfu

def Solve2(gfu, fes, u, v):
    b = CoefficientFunction((2*y*(1-x*x),-2*x*(1-y*y)))
    gaussp = exp(-6*((x+0.5)*(x+0.5)+y*y))-exp(-6*((x-0.5)*(x-0.5)+y*y))
    a = BilinearForm(fes, symmetric=False)
    a += 0.01*grad(u)*grad(v)*dx + b*grad(u)*v*dx
    a.Assemble()
    m = BilinearForm(fes, symmetric=False)
    m += u*v*dx
    m.Assemble()
    mstar = m.mat.CreateMatrix()
    mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()
    invmstar = mstar.Inverse(freedofs=fes.FreeDofs())
    f = LinearForm(fes)
    f += gaussp*v*dx
    f.Assemble()
    res = dt * f.vec - dt * a.mat * gfu.vec
    gfu.vec.data += invmstar * res
    return gfu


def CalcError(gfu, new_mesh):
    flux = 0.01*grad(gfu)
    # interpolate finite element flux into H(div) space:
    gf_flux.Set (flux)

    # Gradient-recovery error estimator
    err = 1/0.01*(flux-gf_flux)*(flux-gf_flux)
    elerr = Integrate (err, new_mesh, VOL, element_wise=True)

    maxerr = max(elerr)
    #print ("maxerr = ", maxerr)

    for el in new_mesh.Elements():
        new_mesh.SetRefinementFlag(el, elerr[el.nr] > 0.25*maxerr)





time = 0; tend = 2; cnt = 0; nsamples = 10; sample_int = int(floor(tend / dt / nsamples)+1)
gfut = GridFunction(gfu_main.space,multidim=0)
gfut.AddMultiDimComponent(gfu.vec)
while time < tend - 0.5 * dt:
    space_flux = HDiv(new_mesh, order=2, autoupdate=True)
    gf_flux = GridFunction(space_flux, "flux", autoupdate=True)
    with TaskManager():
        while fes.ndof < 100000:
            Solve(gfu, fes, u, v)
            CalcError(gfu, new_mesh)
            new_mesh.Refine()
        gfu_main.Set(Solve2(gfu, fes, u, v))
        print("\r",time,end="")
        new_mesh = Mesh(mesh.ngmesh.Copy())
        fes = H1(new_mesh, order=3, dirichlet="bottom|right|left|top", autoupdate=True)
        gfu = GridFunction(fes, autoupdate=True)
        gfu.Set(gfu_main)
        if cnt % sample_int == 0:
            gfut.AddMultiDimComponent(gfu_main.vec)
        u,v = fes.TnT()
        cnt += 1; time = cnt * dt

Draw(gfut, mesh, interpolate_multidim=True, animate=True)

heat3.ipynb (6.1 KB)

7

Hi Joachim

Thanks for your help! I was able to fix it eventually. I am putting the code here so other people can also use it.

from ngsolve import *
from netgen.geom2d import SplineGeometry
from ngsolve.webgui import Draw
from netgen.occ import *
from netgen.webgui import Draw as DrawGeo
shape = Rectangle(2,2).Face().Move((-1,-1,0))
shape.edges.Min(X).name="left"
shape.edges.Max(X).name="right"
shape.edges.Min(Y).name="bottom"
shape.edges.Max(Y).name="top"

time = 0.0
dt = 0.01

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.25))
fes_main = H1(mesh, order=3, dirichlet="bottom|right|left|top")
gfu_main = GridFunction(fes_main)
new_mesh = Mesh(mesh.ngmesh.Copy())
fes = H1(new_mesh, order=3, dirichlet="bottom|right|left|top", autoupdate=True)
u,v = fes.TnT()

gfu_new = GridFunction(fes, autoupdate=True)
gfu_old = GridFunction(fes, autoupdate=True)
gfu_old.Set((1-y*y)*x)



def Solve(gfu_old, fes, u, v):
    b = CoefficientFunction((2*y*(1-x*x),-2*x*(1-y*y)))
    gaussp = exp(-6*((x+0.5)*(x+0.5)+y*y))-exp(-6*((x-0.5)*(x-0.5)+y*y))
    a = BilinearForm(fes, symmetric=False)
    a += 0.01*grad(u)*grad(v)*dx + b*grad(u)*v*dx
    a.Assemble()
    m = BilinearForm(fes, symmetric=False)
    m += u*v*dx
    m.Assemble()
    mstar = m.mat.CreateMatrix()
    mstar.AsVector().data = m.mat.AsVector() + dt * a.mat.AsVector()
    invmstar = mstar.Inverse(freedofs=fes.FreeDofs())
    f = LinearForm(fes)
    f += gaussp*v*dx
    f.Assemble()
    res = dt * f.vec - dt * a.mat * gfu_old.vec
    gfu_new.vec.data = gfu_old.vec + invmstar * res
    return gfu_new


def CalcError():
    flux = 0.01*grad(gfu_new)
    # interpolate finite element flux into H(div) space:
    gf_flux.Set(flux)

    # Gradient-recovery error estimator
    err = 1/0.01*(flux-gf_flux)*(flux-gf_flux)
    elerr = Integrate (err, new_mesh, VOL, element_wise=True)

    maxerr = max(elerr)
    #print ("maxerr = ", maxerr)

    for el in new_mesh.Elements():
        new_mesh.SetRefinementFlag(el, elerr[el.nr] > 0.25*maxerr)





cnt = 1;time = cnt*dt; tend = 2;  nsamples = 10; sample_int = int(floor(tend / dt / nsamples)+1)
gfut = GridFunction(gfu_main.space,multidim=0)
gfu_main.Set(gfu_old)
gfut.AddMultiDimComponent(gfu_main.vec)
while time < tend - 0.5 * dt:
    space_flux = HDiv(new_mesh, order=2, autoupdate=True)
    gf_flux = GridFunction(space_flux, "flux", autoupdate=True)
    with TaskManager():
        while fes.ndof < 100000:
            Solve(gfu_old, fes, u, v)
            CalcError()
            new_mesh.Refine()
    gfu_main.Set(Solve(gfu_old, fes, u, v))
    print("\r",time,end="")
    new_mesh = Mesh(mesh.ngmesh.Copy())
    fes = H1(new_mesh, order=3, dirichlet="bottom|right|left|top", autoupdate=True)
    gfu_old = GridFunction(fes, autoupdate=True)
    gfu_new = GridFunction(fes, autoupdate=True)
    gfu_old.Set(gfu_main)
    if cnt % sample_int == 0:
        gfut.AddMultiDimComponent(gfu_old.vec)
    u,v = fes.TnT()
    cnt += 1;    time = cnt * dt

Draw(gfut, mesh, interpolate_multidim=True, animate=True)