How to assemble separately the mass and stiffness matrix

360
views
0
8 months ago by
Dear all,
I have been working with some nonlinear equations and wanted to try to implement a predictor corrector scheme, however in order to do so I would need to access both the mass and the stiffness matrix arrays after I assemble them.
I was playing a little bit with FEniCS and I think that I am stuck.

I provide the following MWE to show what I mean:

from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import sympy as sy
from sympy import init_printing
init_printing()
%matplotlib inline

dt = 0.1
mesh = IntervalMesh(2, 0, 1)
P1 = FiniteElement('CG', interval, 2)
V = FunctionSpace(mesh, MixedElement([P1, P1]))

v_2, v_1 = TestFunctions(V)
T, P = TrialFunctions(V)

bcl_T = 'near(x[0],0)'
bc = DirichletBC(V.sub(0), 300, bcl_T)

T_0 = Constant('30')
P_0 = Constant('2900')
T_n = interpolate(T_0, V.sub(0).collapse())
P_n = interpolate(P_0, V.sub(1).collapse())

T_en = Constant('30')
P_en = Constant('2900')

F = (P + T)*v_1*dx - (P_n + T_n)*v_1*dx - dt*dot(grad(P), grad(v_1))*dx - \
dt*(P - P_en)*v_1 * ds(1) + \
(P + T)*v_2*dx - (P_n + T_n)*v_2*dx - dt*0.1*dot(grad(T), grad(v_2))*dx - \

a_u_v = lhs(F)
A = assemble(a_u_v)
L_v = rhs(F)
b = assemble(L_v)
bc.apply(A)
bc.apply(b)

# Separeted Version

M_wf = (P + T)*v_1*dx + (P + T)*v_2*dx - dt*P*v_1*ds(1) - dt*T*v_2*ds(2)
F_wf = - (P_n + T_n)*v_1*dx - (P_n + T_n)*v_2*dx + dt*P_en*v_1*ds(1) + dt*T_en*v_2*ds(2)
M = assemble(lhs(M_wf))
K = assemble(lhs(K_wf))
F_ = assemble(rhs(F_wf))
bc.apply(M+K)
bc.apply(F_)

u = Function(V)
w = Function(V)
solve(A, u.vector(), b)
solve((M+K), w.vector(), F_)

u.compute_vertex_values(),w.compute_vertex_values()


I expected that I would have the same results regardless the way I set up the problem, however the output of such code is:

(array([  300.,   300.,   300.,  2630.,  2630.,  2630.]),
array([ -4.31058821e+17,  -4.31058821e+17,  -4.31058821e+17,
4.31058821e+17,   4.31058821e+17,   4.31058821e+17]))​
I suspect that I am missing something, if anyone could please point me a direction I would appretiate a lot.

All the best, Murilo.
Community: FEniCS Project

3
8 months ago by
Change the last few lines to this

M = assemble(lhs(M_wf)+ lhs(K_wf))
F_ = assemble(rhs(F_wf))
bc.apply(M)
bc.apply(F_)

u = Function(V)
w = Function(V)
solve(A, u.vector(), b)
solve(M, w.vector(), F_)

print u.compute_vertex_values()
print w.compute_vertex_values()​

and it seems to work

fenics@eb194e04ecb5:~/shared\$ python test.py
[  300.   300.   300.  2630.  2630.  2630.]
[  300.   300.   300.  2630.  2630.  2630.]

or like this also works

M = assemble(lhs(M_wf))
K = assemble(lhs(K_wf))
F_ = assemble(rhs(F_wf))
A1 = M + K
bc.apply(A1)
bc.apply(F_)

u = Function(V)
w = Function(V)
solve(A, u.vector(), b)
solve(A1, w.vector(), F_)

print u.compute_vertex_values()
print w.compute_vertex_values()

Dear Praveen C,
Thank you very much for you answer.
It works indeed, however, I am sorry for not being so clear about what I would like to have.
I wanted to have the two matrices, for the mass and stiffness, separately as I would need to perform the calculations on the matrix themselves.
I was only solving for the vector to check if I would have the same outcome if I assembled them separately or if I assembled in the more usual way.
With your given answer the two solutions are the same, however I still have only one matrix with both contributions from the mass and stiffness matrices. I would like to have two matrix. Do you think this would be possible?
Anyway thank you very much for your kind help!
All the best, Murilo.

written 8 months ago by Murilo Moreira
1
I added a second way with M and K. Does that not work for you ?
written 8 months ago by Praveen C
Wow! It works on my simple MWE! Thank you very much. However, for my more complex problem it gives still different results.
I will need to test a little bit more.
Thank you very much for your help!
All the best, Murilo.
written 8 months ago by Murilo Moreira
And do you have an explanation of why it works this way and not the others?
Thank you.
written 8 months ago by Murilo Moreira
1
I think if you do bc.apply(M+K) it creates temporary matrix containing M+K on which bc is applied. Your original matrices are not changed.
written 8 months ago by Praveen C
It makes sense! Thank you very much!
written 8 months ago by Murilo Moreira