Ctrl+K

Solver Layers

19.1. Solver Layers#

Current situation:

  • We implement a time-stepper

  • inside the time-stepper we run some Newton-like method to solve a non-linear equation.

  • Newton’s method computes residual and tangent, and uses some direct solver to invert.

Layer approach:

  • We create a Preconditoner object for the BilinearForm object

  • We create a LinearSolver object from the BilinearForm.mat and the Preconditioner object

  • We create a NonlinearSolver from the BilinearForm and the LinearSolver object (or LinearSolver-class+args)

  • We create a TimeStepper from a BilinearForm object and NonLinearSolver class+args

19.1.1. Example:#

setting up the model:

from ngsolve import *
from ngsolve.webgui import Draw

from netgen.occ import *
shape = Rectangle(1,0.1).Face()
shape.edges.Max(X).name="right"
shape.edges.Min(X).name="left"
shape.edges.Max(Y).name="top"
shape.edges.Min(Y).name="bot"
# shape.vertices.Min(X+Y).maxh=0.01
# shape.vertices.Min(X-Y).maxh=0.01
shape.vertices.Min(X+Y).hpref=1
shape.vertices.Min(X-Y).hpref=1

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

E, nu = 210, 0.2
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))

def C(u):
    F = Id(2) + Grad(u)
    return F.trans * F

def NeoHooke (C):
    return 0.5*mu*(Trace(C-Id(2)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)

factor = Parameter(0.5)
force = CoefficientFunction( (0,factor) )

fes = H1(mesh, order=4, dirichlet="left", dim=mesh.dim)
u,v  = fes.TnT()

a = BilinearForm(fes, symmetric=True)
a += Variation(NeoHooke(C(u)).Compile()*dx)
a += Variation((-InnerProduct(force,u)).Compile()*dx)
gfu = GridFunction(fes)
gfu.vec[:] = 0

19.1.1.1. Usual approach#

Call Newton’s method, with a couple of arguments:

solvers.Newton(a=a, u=gfu, inverse="sparsecholesky", printing=False)
scene = Draw (C(gfu)[0,0]-1, mesh, deformation=gfu, min=-0.1, max=0.1)

19.1.1.2. Solver layers approach:#

pre = preconditioners.MultiGrid(a)
linsolve = solvers.CGSolver(a.mat, pre, maxiter=20)
nlsolve = nonlinearsolvers.NewtonSolver(a=a, u=gfu, solver=linsolve)

gfu.vec[:] = 0
nlsolve.Solve(printing=False);

scene = Draw (C(gfu)[0,0]-1, mesh, deformation=gfu, min=-0.1, max=0.1);

Variation keeps the expression for the energy. This is not yet supported for timestepping. We can use symbolic differentiation to convert it to variational forms:

a = BilinearForm(NeoHooke(C(u)).Diff(u, v)*dx - force*v*dx)
gfu = GridFunction(fes)
a.Assemble()  # we need that to construct a CGSolver object outside (TODO: solver_cls, solver_args)

pre = preconditioners.MultiGrid(a)
linsolve = solvers.CGSolver(a.mat, pre, maxiter=20)
nlsolve = nonlinearsolvers.NewtonSolver(a=a, u=gfu, solver=linsolve)

# nlsolve = nonlinearsolvers.NewtonSolver(a=a, u=gfu, lin_solver_cls=solvers.CGSolver, lin_solver_args={ "pre":pre })

gfu.vec[:] = 0
nlsolve.Solve(printing=False);

Draw (C(gfu)[0,0]-1, mesh, deformation=gfu, min=-0.1, max=0.1);

19.1.2. Timesteppers#

  • User provides the time-dependent equation as a variational form

  • Brand new: time-derivative u.dt

force = CF( (0,-1))
eq = NeoHooke(C(u)).Diff(u,v) * dx - force * v * dx + u.dt.dt * v * dx

ts = timestepping.Newmark(equation=eq, dt=1e-1)

gfu = GridFunction(fes)
scene = Draw(gfu, deformation=True, center=(0,-0.3), radius=0.6)

def callback(t, gfu):
    # print("t = ", t, "displacement = ", gfu(mesh(1,0.05)))
    scene.Redraw()
ts.Integrate(gfu, start_time=0, end_time=10, callback=callback);

19.1.2.1. New features to work on exptression trees#

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
fes = H1(mesh, order=3, dirichlet=".*")
u,v = fes.TnT()

equ = u.dt*v*dx + 1e-3*u*v*dx + grad(u)*grad(v)*dx - v*dx
print (equ)

coef binary operation '*', real
  coef trial-function diffop = Id order-dt=1, real
  coef test-function diffop = Id, real
 VOL
coef binary operation '*', real
  coef scale 0.001, real
    coef trial-function diffop = Id, real
  coef test-function diffop = Id, real
 VOL
coef innerproduct, fix size = 2, real
  coef trial-function diffop = grad, real, dim=2
  coef test-function diffop = grad, real, dim=2
 VOL
-1*(coef test-function diffop = Id, real
) VOL

we can extract all proxies (trial- or test-functions) from the variational form:

allproxies = equ.GetProxies(trial=True)
print ("proxies:")
for proxy in allproxies: print(proxy, " python dtorder = ", proxy.dt_order)

proxies:
coef trial-function diffop = Id order-dt=1, real
  python dtorder =  1
coef trial-function diffop = Id, real
  python dtorder =  0
coef trial-function diffop = grad, real, dim=2
  python dtorder =  0

proxiesdt0 = [ prox for prox in allproxies if prox.dt_order==0 ]
proxiesdt1 = [ prox for prox in allproxies if prox.dt_order==1 ]

print ("order 0 proxies:")
print (*proxiesdt0)

print ("order 1 proxies:")
print (*proxiesdt1)

order 0 proxies:
coef trial-function diffop = Id, real
 coef trial-function diffop = grad, real, dim=2

order 1 proxies:
coef trial-function diffop = Id order-dt=1, real

The time-dependent equation is

\[ \int \partial_t u v + uv + \nabla u \nabla v - 1 v \; dx = 0 \]

An implicit Euler time-stepping method is to solve

\[ \int \frac{u - u_n}{\tau} v + uv + \nabla u \nabla v - 1 v \; dx = 0 \]

where the unknown \(u\) is the value of the new time-step \(u_{n+1}\).

So we have to replace the

\[ \partial_t u \rightarrow \frac{u - u_n}{\tau} \]
tau = Parameter(0.1)
gfuold = GridFunction(fes)

repl = {}
for prox in allproxies:
    if prox.dt_order==1:
        repl[prox] = 1/tau*(prox.anti_dt - prox.anti_dt.ReplaceFunction(gfuold))

for key,val in repl.items():
    print ("replace:", key)
    print ("by\n", val)

replace: coef trial-function diffop = Id order-dt=1, real

by
 coef binary operation '*', real
  coef binary operation '/', real
    coef 1, real
    coef N5ngfem28ParameterCoefficientFunctionIdEE, real
  coef binary operation '-', real
    coef trial-function diffop = Id, real
    coef N6ngcomp31GridFunctionCoefficientFunctionE, real

ImplEuler_equ = equ.Replace(repl)
print (ImplEuler_equ)

coef binary operation '*', real
  coef binary operation '*', real
    coef binary operation '/', real
      coef 1, real
      coef N5ngfem28ParameterCoefficientFunctionIdEE, real
    coef binary operation '-', real
      coef trial-function diffop = Id, real
      coef N6ngcomp31GridFunctionCoefficientFunctionE, real
  coef test-function diffop = Id, real
 VOL
coef binary operation '*', real
  coef scale 0.001, real
    coef trial-function diffop = Id, real
  coef test-function diffop = Id, real
 VOL
coef innerproduct, fix size = 2, real
  coef trial-function diffop = grad, real, dim=2
  coef test-function diffop = grad, real, dim=2
 VOL
-1*(coef test-function diffop = Id, real
) VOL

bfIE = BilinearForm(ImplEuler_equ)
bfIE.Assemble()
pre = preconditioners.MultiGrid(bfIE)
linsolve = krylovspace.CGSolver(bfIE.mat, pre, maxiter=20)

gfu = GridFunction(fes)
nlsolve = nonlinearsolvers.NewtonSolver(a=bfIE, u=gfu, solver=linsolve)

scene = Draw(gfu, deformation=True, scale=10)
gfuold.vec[:] = 0

from time import sleep
for i in range(20):
    nlsolve.Solve()
    gfuold.vec[:] = gfu.vec
    scene.Redraw()
    sleep(0.2)
Contents