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 ]
