How to Implementing x * dA/dx Terms in FiPy for Coupled PDEs

63 Views Asked by At

I am relatively new to FiPy and have spent some time reading the FAQ and examining examples. My goal is to solve a set of coupled partial differential equations (PDEs), one of which is represented by the following equations:

q1R = TransientTerm(var=A1) == -((g/2)* A4).grad[0] +((g/2)*A3).grad[1]- ((g/2)*A6.grad[0]) - ((g/2)*A5.grad[1]) - ((g/2)*(x-(nx*dx)/2)*A4) + ((g/2)*(y-(ny*dy)/2)*A3) - ((g/2)*(x-(nx*dx)/2)*A6) - ((g/2)*(y-(ny*dy)/2)*A5)- k*((x-(nx*dx)/2)*A1.grad[0]+(y-(ny*dy)/2)*A1.grad[1]) 
eq1I = TransientTerm(var=A2) == ((g/2)* A3).grad[0] +((g/2)*A4).grad[1]+ ((g/2)*A5.grad[0]) - ((g/2)*A6.grad[1]) + ((g/2)*(x-(nx*dx)/2)*A3) + ((g/2)*(y-(ny*dy)/2)*A4) - ((g/2)*(x-(nx*dx)/2)*A5) - ((g/2)*(y-(ny*dy)/2)*A6) - k*((x-(nx*dx)/2)*A2.grad[0]+(y-(ny*dy)/2)*A2.grad[1]) 
eq2R = TransientTerm(var=A3) == -((g/2)* A8).grad[0] +((g/2)*A2).grad[0]- ((g/2)*A7.grad[1]) + ((g/2)*A1.grad[1]) - ((g/2)*(x-(nx*dx)/2)*A2) - ((g/2)*(x-(nx*dx)/2)*A8) - ((g/2)*(y-(ny*dy)/2)*A1) - ((g/2)*(y-(ny*dy)/2)*A7)- k*((x-(nx*dx)/2)*A3.grad[0]+(y-(ny*dy)/2)*A3.grad[1])

eq2I = TransientTerm(var=A4) == -((g/2)* A8).grad[1] +((g/2)*A2).grad[1]+ ((g/2)*A7.grad[0]) - ((g/2)*A1.grad[0]) + ((g/2)*(x-(nx*dx)/2)*A1) - ((g/2)*(y-(ny*dy)/2)*A8) - ((g/2)*(y-(ny*dy)/2)*A2) + ((g/2)*(x-(nx*dx)/2)*A7)- k*((x-(nx*dx)/2)*A4.grad[0]+(y-(ny*dy)/2)*A4.grad[1]) 

eq3R = TransientTerm(var=A5) == -((g/2)*(x-(nx*dx)/2)*A2) - ((g/2)*(x-(nx*dx)/2)*A8) + ((g/2)*(y-(ny*dy)/2)*A1) + ((g/2)*(y-(ny*dy)/2)*A7) + ((g/2)*A2).grad[0] - ((g/2)*A1.grad[1]) - ((g/2)* A8).grad[0] + ((g/2)*A7.grad[1])- k*((x-(nx*dx)/2)*A5.grad[0]+(y-(ny*dy)/2)*A5.grad[1]) 

eq3I = TransientTerm(var=A6) == ((g/2)*(x-(nx*dx)/2)*A1) + ((g/2)*(x-(nx*dx)/2)*A7) + ((g/2)*(y-(ny*dy)/2)*A2) + ((g/2)*(y-(ny*dy)/2)*A8) - ((g/2)*A1).grad[0] - ((g/2)*A2.grad[1]) + ((g/2)* A7).grad[0] + ((g/2)*A8.grad[1])- k*((x-(nx*dx)/2)*A6.grad[0]+(y-(ny*dy)/2)*A6.grad[1]) 

eq4R = TransientTerm(var=A7) == - ((g/2)*(x-(nx*dx)/2)*A4) + ((g/2)*(y-(ny*dy)/2)*A3) - ((g/2)*(x-(nx*dx)/2)*A6) - ((g/2)*(y-(ny*dy)/2)*A5) + ((g/2)* A4).grad[0] - ((g/2)*A3).grad[1] + ((g/2)*A6.grad[0]) + ((g/2)*A5.grad[1])- k*((x-(nx*dx)/2)*A7.grad[0]+(y-(ny*dy)/2)*A7.grad[1])

eq4I = TransientTerm(var=A8) == + ((g/2)*(x-(nx*dx)/2)*A3) + ((g/2)*(y-(ny*dy)/2)*A4) + ((g/2)*(x-(nx*dx)/2)*A5) - ((g/2)*(y-(ny*dy)/2)*A6) - ((g/2)* A3).grad[0] - ((g/2)*A4).grad[1] - ((g/2)*A5.grad[0]) + ((g/2)*A6.grad[1])- k*((x-(nx*dx)/2)*A8.grad[0]+(y-(ny*dy)/2)*A8.grad[1]) 

eqn = eq1R & eq1I & eq2R & eq2I & eq3R & eq3I & eq4R & eq4I

I am encountering difficulty in implementing terms like x * A1/x and x*A2/x . Also,how can I ensure FiPy recognizes that the 'x' in the derivative and the 'x' multiplying the derivative term are the same in the expression: x * A1/x.I have initialized variables A1, A2, A3, and A4 as CellVariables(with gaussian initial conditions, not given here):

A1 = CellVariable(name="A1", mesh=mesh, hasOld=True)
A2 = CellVariable(name="A2", mesh=mesh, hasOld=True)
A3 = CellVariable(name="A3", mesh=mesh, hasOld=True)
A4 = CellVariable(name="A4", mesh=mesh, hasOld=True)
x = mesh.x
y = mesh.y

To solve above equations, I am looping over a certain number of time steps:

for t in range(2000):
    A1.updateOld()
    A2.updateOld()
    A3.updateOld()
    A4.updateOld()

    eqn.solve(dt=1.e-2)
    

However, the solution obtained does not match the expected result (verified with Mathematica). I am having trouble pinpointing the issue and haven't found relevant solutions in my search other than this Stack Overflow post which does not fully resolve my issues. So to summarize:
Q1) How to implement terms like: x *A1/x, x *A2/x, x *A2/y, etc. in an equation of A1?
Q2) How will FiPy know that the 'x' in the derivative and the 'x' multiplying the derivative term above are the same?
Q3) How to define x and y on a 2D mesh?

I would greatly appreciate any feedback, examples, or guidance on solving such equations with FiPy. Thanks in advance!

0

There are 0 best solutions below