### System of PDEs

139

views

0

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:

File attached: test.xml (1.36 MB)
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.

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)
F=((L-L_n)/k)*v1*dx+kl*r0*L*v1*dx-Dl*alphaL*(L0-L)*v1*ds(1)+Dl*dot(grad(L),grad(v1))*dx\
+((H-H_n)/k)*v2*dx+kh*r0*H*v2*dx-Dh*alphaH*(H0-H)*v2*ds(1)+Dh*dot(grad(H),grad(v2))*dx
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?

Thank you in advance

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)
F=((L-L_n)/k)*v1*dx+kl*r0*L*v1*dx-Dl*alphaL*(L0-L)*v1*ds(1)+Dl*dot(grad(L),grad(v1))*dx
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

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.

I don't think it should.

written
9 weeks ago by
Nima

### 1 Answer

0

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).

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.

\[

\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.