Question on residual update costs

Guosheng Fu
Topic Author
Offline
Elite Member

4 years 7 months ago #2548 by Guosheng Fu

Question on residual update costs was created by Guosheng Fu

Hello,

I am trying to speed up my HDG Naver Stokes solver and need some help on technical details.
In my 2D lid driven cavity problem, I am using a sparse cholesky factorization to solve the linear system,
and all integration terms are already appended with the command

Code:

.Compile(True,True)

.

I found that for a particular problem (k=2 on 128 X 128 grid), the cost of residual update is about twice more expensive than the linear system solver. So I am looking into alternative ways to speed up the residual updates. The residual term contains the nonlinear convection and linear viscosity operators.
I tried two ways for the viscosity operator updates:

Code:

visc = LinearForm(V)
visc += Grad(uh)*Grad(v)*dx
...
visc.Assemble()

and

Code:

visc =BinearForm(V, nonassemble=True)
visc += Grad(u)*Grad(v)*dx
...
res = visc.mat*uh.vec

Surprisingly, the second BilinearForm approach is about 4 times faster than the LinearForm approach. Why is this the case? I thought the two approaches shall be the same (is it faster to take derivatives of proxi functions than grid functions??).

Also, attached is the snippet for my residuals:

Code:

#### diffusion update: very expensive                                        
visc = BilinearForm(fes, nonassemble=True)                                     
visc += (-nu*InnerProduct(gradu,gradv)).Compile(True,True)*dx                  
visc += (nu*(gradu*n*tang(v-vhat)+gradv*n*tang(u-uhat)                         
     -penalty*tang(u-uhat)*tang(v-vhat))).Compile(True,True)*dx(element_boundary=True)
                                                                                
conv = BilinearForm(fes, nonassemble=True)                                      
conv += (-InnerProduct(Grad(u)*u, v)).Compile(True,True)*dx                     
vminus = IfPos(u*n, v.Other(), v)                                     
conv += u2*n*(u-u.Other())*vminus.Compile(True,True)*dx(skeleton=True)

....
res.data = visc.mat*gfu.vec + conv.mat*gfu.vec     

Are there immediate ways to further speed up this step? Thank you.

Best,
Guosheng

schruste
Offline
Administrator

4 years 7 months ago #2552 by schruste

Replied by schruste on topic Question on residual update costs

Hi Guosheng,

Why don't you assemble the viscosity bilinear form and apply the matrix? If you do so, take the nonsymmetric storage of the matrix (it has a faster matrix-vector multiplication than the storage-saving symmetric storage).

For the convection operator there are different tricks around. If you have an Hdiv space in your HDG discretization it is typically more efficient to project into an L2-space, apply the convection there and project the functional back. I recommend having a look at the two model templates for NavierStokes which you can find here: github.com/ngsolve/modeltemplates

Best,
Christoph

Guosheng Fu
Topic Author
Offline
Elite Member

4 years 7 months ago #2554 by Guosheng Fu

Replied by Guosheng Fu on topic Question on residual update costs

Thanks, Christoph. This is helpful. I thought the matrix version shall be slower than assembly on the fly... anyway, I ended up not updating the viscous part on the right hand side.

Also, I (re)found a bug in my code when applying Dirichlet BC conditions for static condensed linear system and found out the remedy in bvp.py (just put it here in case I forgot about it later):
we shall first set zero all internal degrees of freedom before applying the boundary conditions which totally make sense.

# zero local dofs (THIS IS IMPORTANT FOR CONDENSATION!!!!)
innerbits = fes.FreeDofs(False) & ~fes.FreeDofs(True)
Projector(innerbits, False).Project(gfu.vec)
rhs.data -= a.mat*gfu.vec # add dirichlet BC

Time to create page: 0.104 seconds