dolfin team mailing list archive

[Vincenzo C <725208@xxxxxxxxxxxxxxxxxx>]

Public bug reported:

Hi, I'm trying to solve a periodic channel (with stokes eq.) with periodic boundary condition, 
but unfortunately I have an oscillating field of pressure (that I suppose) due 
to the assignment of the periodicBC.. how can I solve the problem?
I attach the code.. 


#Stokes taylor-hood periodic channel

from dolfin import *
xf=5
h=2
nx=90
ny=60
mesh=Rectangle(0,0,xf,h,nx,ny,'left')

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q

class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool((abs(x[1]) < DOLFIN_EPS or abs(x[1]-h) <(DOLFIN_EPS)) and on_boundary)

def by1(x):
      return x[0] >(xf-DOLFIN_EPS)

class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(abs(x[0]) < DOLFIN_EPS and on_boundary)
    def map(self, x, y):
        y[1]=x[1]
    def map(self, x, y):
        y[0]=by1(x)


pbc = PeriodicBoundary()
pbc_x=PeriodicBC(V.sub(0), pbc)
pbc_y=PeriodicBC(V.sub(1), pbc)


u00=Constant((0,0))
dbc = DirichletBoundary()
bc0 = DirichletBC(V,u00, dbc)

bcs = [bc0, pbc_x, pbc_y]

(v, q) = TestFunctions(W)
(u, p) = TrialFunctions(W)
f = Constant((0.0019, 0.0))
a = (inner(grad(v), grad(u)) - div(v)*p + q*div(u))*dx
L = inner(v, f)*dx

problem = VariationalProblem(a, L, bcs)
U = problem.solve()

** Affects: dolfin
     Importance: Undecided
         Status: New

-- 
You received this bug notification because you are a member of DOLFIN
Team, which is subscribed to DOLFIN.
https://bugs.launchpad.net/bugs/725208

Title:
  periodic Boundary condition

Status in DOLFIN:
  New

Bug description:
  Hi, I'm trying to solve a periodic channel (with stokes eq.) with periodic boundary condition, 
  but unfortunately I have an oscillating field of pressure (that I suppose) due 
  to the assignment of the periodicBC.. how can I solve the problem?
  I attach the code.. 

  
  #Stokes taylor-hood periodic channel

  from dolfin import *
  xf=5
  h=2
  nx=90
  ny=60
  mesh=Rectangle(0,0,xf,h,nx,ny,'left')

  V = VectorFunctionSpace(mesh, "CG", 2)
  Q = FunctionSpace(mesh, "CG", 1)
  W = V * Q

  class DirichletBoundary(SubDomain):
      def inside(self, x, on_boundary):
          return bool((abs(x[1]) < DOLFIN_EPS or abs(x[1]-h) <(DOLFIN_EPS)) and on_boundary)

  def by1(x):
        return x[0] >(xf-DOLFIN_EPS)

  class PeriodicBoundary(SubDomain):
      def inside(self, x, on_boundary):
          return bool(abs(x[0]) < DOLFIN_EPS and on_boundary)
      def map(self, x, y):
          y[1]=x[1]
      def map(self, x, y):
          y[0]=by1(x)


  
  pbc = PeriodicBoundary()
  pbc_x=PeriodicBC(V.sub(0), pbc)
  pbc_y=PeriodicBC(V.sub(1), pbc)

  
  u00=Constant((0,0))
  dbc = DirichletBoundary()
  bc0 = DirichletBC(V,u00, dbc)

  bcs = [bc0, pbc_x, pbc_y]

  (v, q) = TestFunctions(W)
  (u, p) = TrialFunctions(W)
  f = Constant((0.0019, 0.0))
  a = (inner(grad(v), grad(u)) - div(v)*p + q*div(u))*dx
  L = inner(v, f)*dx

  problem = VariationalProblem(a, L, bcs)
  U = problem.solve()