System of PDEs

9 weeks ago by

So I am trying to solve a system of pdes. I kept getting a trivial answer, so I reduced it to the following two pdes which are not coupled. I have solved them separately and it gives me a non trivial answer. I thought I should get the same thing using the following code but it does not happen:
from dolfin import *
import numpy as np

Bulk = MeshFunction('size_t' , Domain , 'test_physical_region.xml' ) 
Boundary = MeshFunction('size_t', Domain , 'test_facet_region.xml')

vtkfileu1= File('test/L.pvd')
vtkfileu2= File('test/H.pvd')

ds = Measure('ds', subdomain_data=Boundary)
dx = Measure('dx',subdomain_data= Bulk)

Dh =Constant(3.93)
Dl =Constant(29.89)
kl = Constant(2.35e-4)
kh = Constant(5.29e-6)
L0 = Constant(1e-3)
r0 = Constant(0.26)
H0 = Constant(5e-4)
alphaH = Constant(1)
alphaL =Constant(1)

t , dt= 0.0, 1
k = Constant(dt)

P1 = FiniteElement('P', tetrahedron, 1)
element = MixedElement([P1, P1])
Mixed_Space = FunctionSpace(Domain,element)

v1, v2 =TestFunction(Mixed_Space)
U = Function(Mixed_Space)
U_n = Function(Mixed_Space)
L, H= split(U)
L_n, H_n= split(U_n)


inf_norm = 1
while (inf_norm >= 1e-10):
    t += dt
    solve(F == 0, U)
    inf_norm = project(U-U_n,Mixed_Space).vector().max()
    L_,H_ = U_n.split()

    vtkfileu1<< (L_,t)
    vtkfileu2<< (H_,t)

    print t
    print inf_norm​
and here are my .xml files:
File attached: test.xml (1.36 MB)

File attached: test_physical_region.xml (509.85 KB)

File attached: test_facet_region.xml (1.1 MB)

Could someone please tell me what is so different about solving an uncoupled system as opposed to solving each equation separately that I get such different answers?

Thank you in advance

The following is a code that solves just for L:

from dolfin import *
import numpy as np

Bulk = MeshFunction('size_t' , Domain , 'general_boundary_wenrui+5_physical_region.xml' ) 
Boundary = MeshFunction('size_t', Domain , 'general_boundary_wenrui+5_facet_region.xml')

ds = Measure('ds', subdomain_data=Boundary)
dx = Measure('dx',subdomain_data= Bulk)

Dl =Constant(29.89)
kl = Constant(2.35e-4)
L0 = Constant(1e-3)
r0 = Constant(0.26)
alphaL =Constant(1)

t , dt= 0.0, 0.01
k = Constant(dt)

Single_Space = FunctionSpace(Domain,'P', 1)

L = Function(Single_Space)
L_n = Function(Single_Space)
v1 = TestFunction(Single_Space)


inf_norm = 1
while (inf_norm >= 1e-15):
    t += dt
    solve(F == 0, L)
    inf_norm = project(L-L_n,Single_Space).vector().max()
    print t
    print inf_norm

Which gives a completely different answer.

Community: FEniCS Project
In principal, it should give the same result, since the discrete linear system in the monolithic case should look like this:

K_L & 0 \\
0 & K_H
\end{bmatrix}=\begin{bmatrix} R_L\\R_H \end{bmatrix}
whereas in the case of solving separately you should solve the systems
where, if the ordering of dofs is the same, the matrices \(K_L\) and \(K_H\) and the right hand side vectors \(R_L\) and \(R_H\) are identical.

Can you maybe also post the code where you solve the PDEs separately? Perhaps it has something to do with boundary conditions.
written 9 weeks ago by klunkean  
Thanks for the reply.
That is what I thought, until I did this problem. I added a code for solving L separately at the end of my post which uses the same .xml files for mesh.
written 9 weeks ago by Nima  
If your PDE is linear, then doing newtons method is just essentially solving a linear system of equations. So newtons method just simplifies to a linear system (i.e GMRES, CG) when the root finding problem is an affine function.
written 9 weeks ago by Ryan Vogt  
But do you think that is going to change the results?
I don't think it should.
written 9 weeks ago by Nima  

1 Answer

9 weeks ago by
So I just ran both your codes. The result for L look identical with sufficient precision. The image I attached shows the absolute value of the difference of L computed with the monolithic code and L computed without H:

Maybe you got confused because in the coupled solution you write the solution from the last time step to file (U_n and not U).
Please login to add an answer/comment or follow this question.

Similar posts:
Search »