Hi all,
I’d like to solve this system of PDE :
I ended up with this kind of formulation :
from fenics import *
mesh = UnitIntervalMesh(10)
V = VectorFunctionSpace(mesh, 'P', 1, 2)
u = Function(V)
u_n = Function(V)
v = TestFunction(V)
sol = list(split(u))
prev_sol = list(split(u_n))
test = list(split(v))
F = (sol[0] - prev_sol[0])/2*test[0]*dx - 1 * test[0]*dx
F += (sol[0]*sol[1] - prev_sol[0]*prev_sol[1])/2*test[1]*dx - 1 * test[1]*dx
solve(F == 0, u, [])
But this doesn’t converge (residual is nan). I reckon there’s something I’m not doing right here. Thanks for the help !
Rem
Think about what the left-hand side matrix looks like in the first iteration. Since sol[0] is zero, the mass matrix that you get from linearizing the sol[1]*test[1] term is multiplied by a zero coefficient, and all the rows corresponding to test[1] are zero.
Ah that makes sense. So I suppose I should developp the time derivative term above as :
(sol[0]*(sol[1] - prev_sol[1]) + sol[1]*(sol[0] - prev_sol[0]))/2*test[1]*dx
is that right ?
That still has the problem of the sol[0]–test[1] block in the LHS matrix being scaled by zero. There is now a sol[0]–test[1] block, but it is scaled by sol[1], which is also zero, so this reformulation still crashes in the same way.
Another option might be to have soln[1] represent C_1C_2, then recover C_2 as a post-processing step (although that will become undefined when C_1=0). This would completely decouple the system mentioned in the original post, but I assume it’s really an MWE for something more complicated, where the change-of-variables might not make sense.
It seems that the following is working:
from fenics import *
mesh = UnitIntervalMesh(10)
V = VectorFunctionSpace(mesh, 'P', 1, 2)
u = Function(V)
u_n = Function(V)
v = TestFunction(V)
sol = list(split(u))
prev_sol = list(split(u_n))
test = list(split(v))
F = (sol[0] - prev_sol[0])/2*test[0]*dx - 1 * test[0]*dx
# F += (sol[0]*sol[1] - prev_sol[0]*prev_sol[1])/2*test[1]*dx - 1 * test[1]*dx
F += sol[0]*(sol[1] - prev_sol[1])/2*test[1]*dx - 1 * test[1]*dx + 1*sol[1]*test[1]*dx
solve(F == 0, u, [])
what could be the issue then ?
The 1*sol[1]*test[1]*dx is giving you a nonzero mass matrix in the 11 block of the LHS matrix, although I’m not seeing how it follows from applying the product rule to the second equation.
it is equal to the time derivative of sol[0]