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) ) it complains about the following:
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.
Thanks in advance!
Alex
Code:
from ngsolve import *
from netgen.geom2d import *
# Grain
radius = 0.3
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
stokes = nu*InnerProduct( grad(gf_phi*u), grad(gf_phi*v) )
# Alternative with phi that does not work either :(
# stokes = nu*InnerProduct( grad(phi*u), grad(phi*v) )
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()