Explicit assembly and external solve


326
views
1
10 months ago by
AJ  
Hi,

For the validation purpose, I am trying to export FEniCS assembly (Poisson's equation) matrix and vector to Matlab and solve using (u = A\b). However, when I solve externally the solution vector looks different (same values but different sparsity) compared to the solution vector extracted from internal FEniCS solve(A, U, b).

For example, take a look at the following Snippet:

A = assemble(a)
b = assemble(L)
bc.apply(A, b)
sio.savemat('Ab.mat', {'A':A.array(), 'b':b.array()})
u = Function(V)
U = u.vector()
solve(A, U, b)

U = [0
0
0
0
0
0
0
0
-0.479166666666669
-0.276041666666667
-0.244791666666667
-0.244791666666667
-0.276041666666667]


Matlab solve (u = A\b)
u=[0
0
-0.244791666666667
0
-0.479166666666669
-0.276041666666667
0
-0.244791666666667
0
-0.276041666666667
0
0
0].

Does anyone have a clue what am I missing here?
Thanks,
-A







Community: FEniCS Project
I think a good place to start is to verify that A and b are identical between software.
written 10 months ago by pf4d  
Hi pf4d, Thanks for the comment. 

Both A and b are assembled in FEniCS and after solving for U in FEniCS I export them to Matlab and then solve again. So A and b in both places are one and the same. 

I wonder if FEniCS changes something inside during the cal of solve(A, U, b). 

-A
written 10 months ago by AJ  
The "export" process is what I refer to; have you verified post export within MATLAB that these tensor elements are the same?

That is, if you were to manually create a MATLAB matrix and vector from the result of printing the Python tenors "A.array()" and "b.array()", and then solve, what do you get then?
written 10 months ago by pf4d  
Thanks for the comment. I get the same solution vector if I copied values from Python tenors "A.array()" and "b.array()" and solve in Matlab. 

I posed the question wrongly (my apologies).

My concern is when we print solution vector after solving it is:
U = [ 0.  0.  -0.24479167  0. -0.47916667   -0.27604167   0.  -0.24479167  0. -0.27604167  0. 0. 0.]

And when I open "poisson000000.vtu" it is:
U = [ 0. 0. 0. 0. 0. 0. 0. 0. -0.47916667  -0.27604167  -0.24479167  -0.24479167  -0.27604167]

How did the non-zero positions change between two?

If I solve the FEniCS assembled system externally how can I do the same so that I can plot the solution vector using mesh loaded externally?

 
written 10 months ago by AJ  
Oh, I see.  I believe that the internal representation of tensor components has been optimized for speed, and when exporting they change the ordering to something more convenient for post-processing tools...  perhaps.  This is indeed a different question entirely.

Maybe this will help you:

https://github.com/mikaem/fenicstools/wiki/DofMap-plotter
written 10 months ago by pf4d  
That link can surely help to figure out the reordering. However, if you're aware, is there a way to get the assembly matrix and vector without node reordering, even if it is slow?

I am more interested in using assembly routines from FEniCS and then use the in-house DDM solver to solve the system. To do that, I need the Assembly matrix and vector in its original node order.

Thanks for all your comments, I appreciate your help.  
-A

written 10 months ago by AJ  
Adding "parameters['reorder_dofs_serial'] = False" seems to do the tric. 

-A
written 10 months ago by AJ  

1 Answer


0
10 months ago by
pf4d  
I'm not sure what MATLAB is doing, but here is a MWE:

from fenics import *

mesh = UnitIntervalMesh(3)
Q    = FunctionSpace(mesh, 'CG', 1)

u    = TrialFunction(Q)
v    = TestFunction(Q)
x    = Function(Q)
f    = Constant(1.0)

def boundary(x, on_boundary):
  return on_boundary

bc = DirichletBC(Q, 0.0, boundary)

A = assemble( u.dx(0) * v.dx(0) * dx)
b = assemble( f * v * dx)

bc.apply(A,b)

solve(A, x.vector(), b)

print "A:\n", A.array(),          "\n"
print "B:\n", b.array(),          "\n"
print "x:\n", x.vector().array(), "\n"

from scipy.linalg import solve as sci_solve

x_sci = sci_solve(A.array(), b.array())

print "x_sci:\n", x_sci, "\n"
print "|| x - x_sci ||_oo:", \
      (x.vector() - x_sci).max()

It is easy to verify that the tensors and solution are correct for this function space, see for example how to build the tensor from

https://www.amazon.com/J-N-Reddy-Introduction-Element/dp/B008VS0BPC/ref=sr_1_1?s=books&ie=UTF8&qid=1507387554&sr=1-1&keywords=An+Introduction+to+the+Finite+Element+Method+2nd

Hope that helps...

Thank you pf4d! That is helpful. 

If I want to visualize the solution after solving externally, let's say using "pdesurf in Matlab". If I use the same mesh and provide "pdesurf(points, elements, U)". I don't get the same plot as I get from FEniCS.
Because the non-zero solution values are associated with different nodal points.

As I posted in the comment:
Printed soution after solving:
U = [ 0.  0.  -0.24479167  0. -0.47916667   -0.27604167   0.  -0.24479167  0. -0.27604167  0. 0. 0.]

And when I open "poisson000000.vtu" it is:
U = [ 0. 0. 0. 0. 0. 0. 0. 0. -0.47916667  -0.27604167  -0.24479167  -0.24479167  -0.27604167]

-A 
written 10 months ago by AJ  
.vtu files are used for viewing with visit or paraview; save your solution with numpy.savetxt() or scipy.savemat().
written 10 months ago by pf4d  
Or, if you have HDF5 support in Matlab, save in parallel with that format:

https://www.mathworks.com/help/matlab/hdf5-files.html

see my comment here:

https://www.allanswered.com/post/aqox/#pzew

hope that helps!
written 10 months ago by pf4d  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »