# Solution of scalar PDE as parameter for Stokes equation

Hi,

I am trying to solve a modified Stokes equation. In the first step I solve a heat equation for a scalar quantity phi. This quantity should be used as a parameter in the Stokes equation afterwards. I have followed the Tutorials such that I was able to solve the PDE for phi, but I could not make it work for the Stokes equation yet. In line 93 of my code

``  stokes = nu*InnerProduct( grad(gf_phi*u), grad(gf_phi*v) ) ``

``AttributeError: 'ngsolve.fem.CoefficientFunction' object has no attribute 'derivname'``

I am not sure how to solve this. I am allowed to take the gradient of gf_phi and u alone, but I am not allowed to take the gradient of the product gf_phi*u. Does anyone understand why this does not work and how to fix it?

I have added a code example below.

Alex

Code:

``````from ngsolve import *
from netgen.geom2d import *

# Grain
grain_x_origin = 0.5
grain_y_origin = 0.5

# viscosity
nu = 1.00
# domain properties
edge_length = 1.0

# forcing term
forcing = CoefficientFunction(( 1.0, 0. ))

# mesh
max_mesh_size=0.025

# Phase-field
lam = 0.02
gamma = 0.01

n = 10.
K = 25.

# Time stepping parameters
tau = 1e-5
tend = 1e-5

def create_periodic_mesh():
periodic = SplineGeometry()
pnts = [ (0,0), (edge_length,0.0), (edge_length,edge_length), (0.0,edge_length) ]
pnums = [periodic.AppendPoint(*p) for p in pnts]

# Make it periodic in y-direction
lbottom = periodic.Append( ["line", pnums[0], pnums[1]], bc="periodic" )
periodic.Append( ["line", pnums[3], pnums[2]], leftdomain=0, rightdomain=1, copy=lbottom, bc="periodic" )
# Make it periodic in x-direction
right = periodic.Append( ["line", pnums[1], pnums[2]], bc="periodic" )
periodic.Append( ["line", pnums[0], pnums[3]], leftdomain=0, rightdomain=1, copy=right, bc="periodic" )

return Mesh( periodic.GenerateMesh(maxh=max_mesh_size))

def main():

mesh = create_periodic_mesh()

u_initial = CoefficientFunction( ( 1. / (1. + exp( -4. / lam * ( sqrt( (x-grain_x_origin)**2+(y-grain_y_origin)**2) - radius) ) ) ) )

Vphi = Periodic(H1(mesh,order=3))
phi = Vphi.TrialFunction()
v = Vphi.TestFunction()

gf_phi = GridFunction(Vphi)
gf_phi.Set( u_initial )
gf_phi_old = GridFunction(Vphi)
gf_phi_old.Set( u_initial )

a = BilinearForm(Vphi)
phasefield = phi * v + tau * grad(phi) * grad(v) - gf_phi_old * v
phasefield *= dx

a += phasefield
a.Assemble()

Draw(u_initial, mesh, "u_initial")
Draw(gf_phi-u_initial,mesh,"diff")
Draw(gf_phi,mesh,"phi")

# Solve heat equation for some time steps to make it more diffuse
t = 0.
while t < tend:
print ("t={}".format( t ))

solvers.Newton(a, gf_phi, maxerr=1e-12, maxit=10, printing=True)
gf_phi_old.vec.data = gf_phi.vec

t = t + tau
Redraw()
Redraw()

# Solve Stokes (taken from examples)
V = Periodic(VectorH1(mesh,order=3))
Q = Periodic(H1(mesh,order=2))

X = FESpace([V,Q])

u,p = X.TrialFunction()
v,q = X.TestFunction()

# Bilinear form
# Alternative with phi that does not work either :(
stokes -= div(gf_phi*u)*q-div(gf_phi*v)*p
stokes += K / lam * (1.-gf_phi)*n / (gf_phi+n) *u * v
stokes *= dx
a = BilinearForm(X)
a += stokes
a.Assemble()

# forcing term
f = LinearForm(X)
f += SymbolicLFI(gf_phi * forcing * v)
f.Assemble()

# gridfunction for the solution
gf_phi_stokes = GridFunction(X)
velocity = gf_phi_stokes.components[0]
Draw(velocity,mesh,"u",sd=3)
Draw(gf_phi.components[1],mesh,"p",sd=3)
from ngsolve.internal import visoptions
visoptions.scalfunction = "u:0"

# solve Stokes problem for initial conditions:
inv_stokes = a.mat.Inverse(X.FreeDofs())

res = f.vec.CreateVector()
res.data = f.vec - a.mat*gf_phi.vec
gf_phi.vec.data += inv_stokes * res

main()``````

Hi ajaust,

you can use the grad operator for GridFunctions, Trial- and Testfunctions (if they provide this operator). When multiplying a Trialfunction and GridFunction the result is a CoefficientFunction, which in general does not support the grad operator (a CoefficientFunction is a quite generic object).

Therefore, you have to compute the product rule by hand and write the result into your bilinearform, which should be something like (maybe a transpose is missing)

``nu*InnerProduct( OuterProduct(u,grad(gf_phi)) + gf_phi*grad(u),  OuterProduct(v,grad(gf_phi)) + gf_phi*grad(v) )``

Best
Michael

Thanks! This solved the problem. It does not seem as any transpose was necessary.

Best
Alex