Solution of scalar PDE as parameter for Stokes equation

ajaust
Topic Author
Offline
New Member

4 years 9 months ago #2361 by ajaust

Solution of scalar PDE as parameter for Stokes equation was created by ajaust

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

Code:

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

it complains about the following:

Code:

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:

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()

mneunteufel
Offline
Premium Member

4 years 9 months ago #2362 by mneunteufel

Replied by mneunteufel on topic Solution of scalar PDE as parameter for Stokes equation

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)

Code:

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

Best
Michael

The following user(s) said Thank You: ajaust

ajaust
Topic Author
Offline
New Member

4 years 9 months ago #2365 by ajaust

Replied by ajaust on topic Solution of scalar PDE as parameter for Stokes equation

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

Best
Alex

Time to create page: 0.098 seconds