dolfin team mailing list archive

[Melanie Jahny <question147529@xxxxxxxxxxxxxxxxxxxxx>]

New question #147529 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/147529

I'm trying to solve the advection diffusion equation in a rectangle.
I'm using the velocity function computed by stokes brinkmann equation. 
The velocity field describes a flow streaming from the right and the left side into the channel and getting out on
the lower boundary in the middle of the channel
The ligand for the advection diffusion equation flows into the channel by a neumann-boundary condition.
I want it to stream in only from the left side of the boundary, but fenics always computes concentrations
of the ligand on both sides of the channel (already in the first time step). 
How can this happen??? How can I avoid it???

If I use a velocity field flowing only from the left side to the right side this does not occur, but in my
model I have to use a velocity field streaming from both sides. 

Here is the code:
mesh = Rectangle(0,0,154.0,7.8,154,20)

V = FunctionSpace(mesh, "CG", 2)
conc = TrialFunction(V) # solution
conc_prev = Function(V) # solution from previous time step
eta = TestFunction(V)

# Read velocity Function calculated by Stokes-Brinkman-Equation
W = VectorFunctionSpace(mesh, "CG", 2)
velocity = Function(W, "./velocity.xml") 

F_time = inner(eta,(conc-conc_prev))*dx 
F_Diff = dt*(D_const*inner(grad(eta), grad(conc))*dx)  
F_Adv = dt*(inner(eta,inner(velocity, grad(conc)))*dx)

boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)

class NeumannBoundary(SubDomain):
    def inside(self, x, on_boundary):        
         return ((x[0] >= 1.0 and x[0] <= 3.0) and x[1] == y_max and on_boundary)
       
NeumBound = NeumannBoundary()
NeumBound.mark(boundary_parts,0)

bound_flow = Expression("0.87")
bound_nm = dt*bound_flow*eta*ds(0) 

F = F_time + F_Diff + F_Adv - bound_nm 

a = lhs(F) 
L = rhs(F)

A = assemble(a, exterior_facet_domains=boundary_parts) 

solver = cpp.UmfpackLUSolver(A) 
solver.factorize() 

while(t<T):
      
   b = assemble(L, exterior_facet_domains=boundary_parts) 
   solver.solve_factorized(conc.vector(), b)        
   conc_prev.assign(conc)
   t += dt   
      
      
I would appreciate it if anyone could help me. Thank you!





-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.