# Re: 2 coupled convection equations - how to apply coundary conditions?

Welcome, Tanja.

Your declaration of boundary conditions looks fine to me, other than being on
the opposite faces (e.g., x=0 is facesLeft and x=L is facesRight) from your
problem definition.

FiPy is not well suited to solving purely hyperbolic PDEs and the combination
of implicit solvers and steady state solutions means that I'm not particularly
surprised that the "solution" doesn't see your boundary conditions, but I'm
also not really able to explain it.

I've modified your script below to include TransientTerms to hopefully allow
relaxing toward your solution. The B(0,z) =  z boundary condition definitely
has an impact, but there is a severe checkerboard instability. Reducing the
mesh spacing and the time step does not make this instability go away.

I suggest you search our mailing list archives for "hyperbolic" for previous
discussions about solving equations like this, but I don't hold out much hope.

https://www.mail-archive.com/search?q=hyperbolic&l=fipy%40nist.gov

Daniel, Thibault, or James may have more insight than I.

You also might have better luck with a package like CLAWPACK which is designed
for hyperbolic equations.

- Jon

----

from fipy import *
# Define the domain and make a mesh
H = 50.0
L = H
nx = 20
nz = nx
dx = L / nx
dz = H / nz
mesh = Grid2D(dx=dx, dy=dz, nx=nx, ny=nz)
# Create two CellVariable and initialize it to zero:
A = CellVariable(name = "normal stress",mesh = mesh,value = 0.)
B = CellVariable(name = "shear stress",mesh = mesh,value = 0.)
# apply boundary conditions
X, Z = mesh.faceCenters
B.constrain(Z,mesh.facesLeft)
B.constrain(0.,mesh.facesBottom)
B.constrain(0.,mesh.facesRight)
A.constrain(0.,mesh.facesBottom)
A.constrain(0.,mesh.facesRight)

viewer = Viewer(vars=(A, B))
viewer.plot()

# define the system of equations
eq1 = (CentralDifferenceConvectionTerm(var=A,coeff=(1.,0.,)) +
CentralDifferenceConvectionTerm(var=B,coeff=(0.,1.,))) == TransientTerm(var=B)
eq2 = (CentralDifferenceConvectionTerm(var=B,coeff=(1.,0.,)) +
CentralDifferenceConvectionTerm(var=A,coeff=(0.,-1.,))) == TransientTerm(var=A)
coupled = eq1 & eq2
# solve
for step in range(1000):
coupled.solve(dt=.1)

# view solution

viewer.plot()

> On Oct 4, 2017, at 6:02 AM, Tanja Schlemm <schl...@pik-potsdam.de> wrote:
>
> Hello!
>
> I'm new to FiPy and I'm trying to solve a 2-dimensional problem of 2
> coupled equations with 2 variables. I'm having troubles with the
> boundary conditions. I checked the examples, but all examples containing
> convection are 1d.
> I'd be very thankful for some tips on how I can make this work!
>
> The equations I want to solve:
> \partial_x A + \partial_z B = 0
> \partial_x B - \partial_z A = 0
> on a box with 0<x<L and 0<z<H. Boundary conditions are
> A(x,0) = B(x,0) = 0
> A(L,z) = B(L,z) = 0
> B(0,z) =  z
> I tried rewriting my problem in terms of convection terms. Introducing
> the vectors
> u = (1,0),v = (0,1)
> the equations above can be written as convection terms:
> \nabla \cdot (u A) + \nabla \cdot (v B) = 0
> \nabla \cdot (u B) - \nabla \cdot (v A) = 0
>
> My attempt to solve this with FiPy:
> # Define the domain and make a mesh
> H = 50.0
> L = H
> nx = 20
> nz = nx
> dx = L / nx
> dz = H / nz
> mesh = Grid2D(dx=dx, dy=dz, nx=nx, ny=nz)
> # Create two CellVariable and initialize it to zero:
> A = CellVariable(name = "normal stress",mesh = mesh,value = 0.)
> B = CellVariable(name = "shear stress",mesh = mesh,value = 0.)
> # apply boundary conditions
> X, Z = mesh.faceCenters
> B.constrain(Z,mesh.facesRight)
> B.constrain(0.,mesh.facesTop)
> B.constrain(0.,mesh.facesLeft)
> A.constrain(0.,mesh.facesTop)
> A.constrain(0.,mesh.facesLeft)
> # define the system of equations
> eq1 = (ConvectionTerm(var=A,coeff=(1.,0.)) +
> ConvectionTerm(var=B,coeff=(0.,1.))) == 0.
> eq2 = (ConvectionTerm(var=B,coeff=(1.,0.)) +
> ConvectionTerm(var=A,coeff=(0.,-1.))) == 0.
> coupled = eq1 & eq2
> # solve
> coupled.solve()
> # view solution
> viewer = Viewer(vars=A)
> viewer.plot()
> viewer = Viewer(vars=B)
> viewer.plot()
>
> The solution I get is zero everywhere. It doesn't even fulfill the
> boundary conditions... I experimented with some other boundary
> conditions (setting some to 0. and some to 1.), the results were
> somewhat different, but the resulting solution never quite fulfilled the
> boundary conditions.
> So I feel like I'm missing something vital about how boundary conditions
> work with convection terms. Can someone give me a hint how it works and
> what I'm doign wrong?
>
> Kind regards
> Tanja
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]