System of PDEs


139
views
0
9 weeks ago by
Nima  
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)

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


0
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 »