### System of PDEs

139
views
0
9 weeks ago by
Hi,

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

Domain=Mesh('test.xml')
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()
L_.rename('L','L')
H_.rename('H','H')

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

U_n.assign(U)
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?

The following is a code that solves just for L:

from dolfin import *
import numpy as np

Domain=Mesh('general_boundary_wenrui+5.xml')
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()
L_n.assign(L)
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:

$\begin{bmatrix} K_L & 0 \\ 0 & K_H \end{bmatrix} \begin{bmatrix} L\\ H \end{bmatrix}=\begin{bmatrix} R_L\\R_H \end{bmatrix}$
whereas in the case of solving separately you should solve the systems
\begin{align*}
K_LL&=R_L\\
K_HH&=R_H
\end{align*}
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
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