This all seems a bit odd. Inflow boundary conditions for the advection equation should be handled weakly, then you can make an upwind flag to do this as normal. cheers --cjc On 26 May 2016 at 09:19, Mitchell, Lawrence < lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,
On 26 May 2016, at 08:43, Justin Chang <jychang48@gmail.com> wrote:
Hi all,
I am attempting to simulate a chaotic advection equation where the velocity vector field is based on the Arnold–Beltrami–Childress (ABC) flow equation:
v_x = A*sin(z) + C*cos(y) v_x = B*sin(x) + A*cos(z) v_x = C*sin(y) + B*cos(x)
Right now I have the following:
from firedrake import * mesh = UnitCubeMesh(20,20,20) A = 0.08 B = 0.1 C = 0.15 V = VectorFunctionSpace(mesh, "CG", 1) velocity = Function(V) velocity.interpolate(Expression(( "A*sin(2*pi*x[2])+C*cos(2*pi*x[1])", "B*sin(2*pi*x[0])+A*cos(2*pi*x[2])", "C*sin(8*pi*x[1])+B*cos(8*pi*x[0])"), A=A,B=B,C=C))
FWIW, I recommend for interpolating expressions to utilise the new UFL expression interpolation. http://firedrakeproject.org/interpolation.html#ufl-expressions This gets you much better type checking and more useful error messages in the case of typos.
Where A,B,C are user-defined constants (will add such constants inside all sin and cos functions later). I want to have DirichletBC's only applied to faces and regions where the velocity points into the cube (i.e., the "upstream" location).
This is doable, you need to subclass DirichletBC and provide your own property that computes the *nodes* at which the boundary condition should apply.
In the following, I assume you want to apply a bc to the same space that the velocity lives in. I would do the following.
import operator import numpy as np
class UpstreamBC(DirichletBC): def __init__(self, V, value, sub_domain, method="topological", velocity=None): super(UpstreamBC, self).__init__(V, value, sub_domain, method=method) self._velocity = velocity self._all_nodes = super(UpstreamBC, self).nodes
@property def nodes(self): face_velocity = self._velocity.dat.data_ro_with_halos[self._all_nodes] component, comparator = {1: (0, operator.gt), # x == 0 2: (0, operator.lt), # x == 1 3: (1, operator.gt), # y == 0 4: (1, operator.lt), # y == 1 5: (2, operator.gt), # z == 0 6: (2, operator.lt), # z == 1 }[self.sub_domain] return self._all_nodes[comparator(face_velocity[:, component], 0.0)]
@property def node_set(self): return op2.Subset(self._function_space.node_set, self.nodes)
Untested, but you can see the idea. I provide, in addition the the normal arguments an additional velocity field. In the nodes property I use the value of the velocity field to select the appropriate subset of the nodes.
Note that to make this selection simple, I require that the bc object be made with a single subdomain value (that way I can simply determine which component of the velocity field I should be checking). So you need one BC object for each side of the cube. This is possibly slightly inefficient, so we could improve it, but let's try and get everything working first.
The region locations and sizes will vary depending on the constants chosen. Yes I could create separate GMSH files and manually determine where these "inflow" regions will be for every set of A/B/C/etc constants, but I was wondering if there's a way to automate this the Firedrake way, particularly if these constants will vary as a function of time.
The setup I have above is already time-dependent, so you don't need to do anything else.
Hope this works!
Lawrence
-- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916