Returning eigenvectors of a solution to a PDE at each node
Hello, I am currently using firedrake to solve a PDE. The solution to this PDE is vector-valued, and I have no problem finding the solution. However, I'd like to go a step further. I have named the solution "q_soln" and this is a 2-dimesional vector. I have used this vector to construct a 2x2 matrix at the node [0.5,0.5]. I then extracted an eigenvector from this matrix, and then I printed it: import numpy as np from numpy import linalg as la E1 = (1 / np.sqrt(2)) * np.array([[1,0],[0,-1]]) E2 = (1 / np.sqrt(2)) * np.array([[0,1],[1,0]]) # as an example, we evaluate the solution at the point [0.5,0.5], and pick the first eigenvector q_soln_atpoint = q_soln.at([0.5,0.5]) Q_soln = q_soln_atpoint[0] * E1 + q_soln_atpoint[1] * E2 evals, evecs = la.eig(Q_soln) # evec is a 2x2 array with the columns being the e-vectors q_soln_evec = np.zeros((1,2)) q_soln_evec += evecs[0,:] print(q_soln_evec) My question is this: how can I get firedrake to do this automatically for every node in my mesh, and then store it in as a Function "q_soln_evec"? The only reason I used numpy to get this one point is because I'm unfamiliar with firedrake. I bet firedrake has some built-in way of doing what I need here. Thanks, Andrew Hicks
Dear Andrew, I'm afraid that I don't know the most efficient way to do this sort of many small eigensolves, but you could get your solution in `np.ndarray` using `q_soln.dat.data` (`q_soln.dat` returns `PyOP2.Dat` object: https://op2.github.io/PyOP2/user.html#pyop2.Dat). You could also get the mesh data calling `mesh.coordinates.dat.data` . For instance:
from firedrake import * mesh = UnitSquareMesh(2,2) V = VectorFunctionSpace(mesh, "CG", 1) f = Function(V) f.dat.data array([[0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.]]) mesh.coordinates.dat.data array([[0. , 0. ], [0. , 0.5], [0.5, 0. ], [0.5, 0.5], [1. , 0. ], [0. , 1. ], [1. , 0.5], [0.5, 1. ], [1. , 1. ]])
I think you can use these data in your context. Taking a glance at the documentation, `numpy.linalg.eig` seems to take multiple matrices, so I think you can solve the node-wise eigen problems all at once. Thank you, Koki ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Andrew Hicks <ahick17@lsu.edu> Sent: Friday, March 20, 2020 4:23 AM To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] Returning eigenvectors of a solution to a PDE at each node Hello, I am currently using firedrake to solve a PDE. The solution to this PDE is vector-valued, and I have no problem finding the solution. However, I’d like to go a step further. I have named the solution “q_soln” and this is a 2-dimesional vector. I have used this vector to construct a 2x2 matrix at the node [0.5,0.5]. I then extracted an eigenvector from this matrix, and then I printed it: import numpy as np from numpy import linalg as la E1 = (1 / np.sqrt(2)) * np.array([[1,0],[0,-1]]) E2 = (1 / np.sqrt(2)) * np.array([[0,1],[1,0]]) # as an example, we evaluate the solution at the point [0.5,0.5], and pick the first eigenvector q_soln_atpoint = q_soln.at([0.5,0.5]) Q_soln = q_soln_atpoint[0] * E1 + q_soln_atpoint[1] * E2 evals, evecs = la.eig(Q_soln) # evec is a 2x2 array with the columns being the e-vectors q_soln_evec = np.zeros((1,2)) q_soln_evec += evecs[0,:] print(q_soln_evec) My question is this: how can I get firedrake to do this automatically for every node in my mesh, and then store it in as a Function “q_soln_evec”? The only reason I used numpy to get this one point is because I’m unfamiliar with firedrake. I bet firedrake has some built-in way of doing what I need here. Thanks, Andrew Hicks
Hi Andrew, Koki, To solve local eigenproblems I currently call Eigen within par_loops. After specifying a mesh, a tensor field M (which you seek to compute the eigendecompositon of) and your PETSC_ARCH, the following should work. Extension to 3d just means changing the kernel_str appropriately. Cheers, Joe PETSC_ARCH = ... mesh = ... P1_ten = TensorFunctionSpace(mesh, "CG", 1) M = interpolate(..., P1_ten) # Tensor field to compute eigendecomposition of evec = Function(P1_ten) # Tensor field to hold eigenvectors eval = Function(P1_ten) # Tensor field to hold eigenvalues kernel_str = """ #include <Eigen/Dense> using namespace Eigen; void get_eigendecomposition(double EVecs_[4], double EVals_[2], const double * M_) { Map<Matrix<double, 2, 2, RowMajor> > EVecs((double *)EVecs_); Map<Vector2d> EVals((double *)EVals_); Map<Matrix<double, 2, 2, RowMajor> > M((double *)M_); SelfAdjointEigenSolver<Matrix<double, 2, 2, RowMajor>> eigensolver(M); Matrix<double, 2, 2, RowMajor> Q = eigensolver.eigenvectors(); Vector2d D = eigensolver.eigenvalues(); EVecs = Q; EVals = D; } """ kernel = op2.Kernel(kernel_str, 'get_eigendecomposition', cpp=True, include_dirs= ["%s/include/eigen3" % PETSC_ARCH]) op2.par_loop(kernel, P1_ten.node_set, evec.dat(op2.RW), eval.dat(op2.RW), M.dat(op2.READ)) ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Sagiyama, Koki <k.sagiyama@imperial.ac.uk> Sent: 20 March 2020 10:52 To: Andrew Hicks <ahick17@lsu.edu>; firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] Returning eigenvectors of a solution to a PDE at each node Dear Andrew, I'm afraid that I don't know the most efficient way to do this sort of many small eigensolves, but you could get your solution in `np.ndarray` using `q_soln.dat.data` (`q_soln.dat` returns `PyOP2.Dat` object: https://op2.github.io/PyOP2/user.html#pyop2.Dat). You could also get the mesh data calling `mesh.coordinates.dat.data` . For instance:
from firedrake import * mesh = UnitSquareMesh(2,2) V = VectorFunctionSpace(mesh, "CG", 1) f = Function(V) f.dat.data array([[0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.], [0., 0.]]) mesh.coordinates.dat.data array([[0. , 0. ], [0. , 0.5], [0.5, 0. ], [0.5, 0.5], [1. , 0. ], [0. , 1. ], [1. , 0.5], [0.5, 1. ], [1. , 1. ]])
I think you can use these data in your context. Taking a glance at the documentation, `numpy.linalg.eig` seems to take multiple matrices, so I think you can solve the node-wise eigen problems all at once. Thank you, Koki ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Andrew Hicks <ahick17@lsu.edu> Sent: Friday, March 20, 2020 4:23 AM To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] Returning eigenvectors of a solution to a PDE at each node Hello, I am currently using firedrake to solve a PDE. The solution to this PDE is vector-valued, and I have no problem finding the solution. However, I’d like to go a step further. I have named the solution “q_soln” and this is a 2-dimesional vector. I have used this vector to construct a 2x2 matrix at the node [0.5,0.5]. I then extracted an eigenvector from this matrix, and then I printed it: import numpy as np from numpy import linalg as la E1 = (1 / np.sqrt(2)) * np.array([[1,0],[0,-1]]) E2 = (1 / np.sqrt(2)) * np.array([[0,1],[1,0]]) # as an example, we evaluate the solution at the point [0.5,0.5], and pick the first eigenvector q_soln_atpoint = q_soln.at([0.5,0.5]) Q_soln = q_soln_atpoint[0] * E1 + q_soln_atpoint[1] * E2 evals, evecs = la.eig(Q_soln) # evec is a 2x2 array with the columns being the e-vectors q_soln_evec = np.zeros((1,2)) q_soln_evec += evecs[0,:] print(q_soln_evec) My question is this: how can I get firedrake to do this automatically for every node in my mesh, and then store it in as a Function “q_soln_evec”? The only reason I used numpy to get this one point is because I’m unfamiliar with firedrake. I bet firedrake has some built-in way of doing what I need here. Thanks, Andrew Hicks
participants (3)
-
Andrew Hicks
-
Sagiyama, Koki
-
Wallwork, Joe