Dear Lawrence, Yes, I have a term where the test function and the basis function are evaluated in separate integrals, this is why I wrote Q*mu*dx. Can I assemble them separately? E.g. Q1 = assemble(v*dx) Q2 = assemble(mu*dx) and then multiply the arrays? This is how I did it in my own FEM code, but I want to understand how to do it using Firedrake as well. Thanks. Best, Anna. On 19/10/15 11:46, Lawrence Mitchell wrote:

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1

Dear Anna,

On 16/10/15 19:03, Anna Kalogirou wrote:

...

Dear all,

Can I use an assembled matrix into a bilinear form?

E.g. I want to write something like

mu1 = Function(V) mu = TrialFunction(V) v = TestFunction(V)

Q = assemble(v*dx) lhs = inner(grad(mu),grad(v))*dx + Q*mu*dx I'm confused, Q*mu*dx is not a bilinear form (it has a trial function but no test function).

Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)

iQEcBAEBAgAGBQJWJMoeAAoJECOc1kQ8PEYvWEQIALKqAik+xEr+64+uzuASvJJM Gl/+kr/hIV8D39a0cYYT8uZb3l0qhW8CQS4ibz0rW4vCjYDTpNhiqnScU+8I9U2F +aOcTZAu48P0DgU3GwoamRJ2VOtvTROFyZARLNwTrQNn7Bokua0CGvxbQB+STF4D kmaL27Cp6kZWTynHVrFR6PZuUBPzYrVsGZL2EASi092sIXQS0Uvc2WE9MpioKyHh 6a2LR+8oadFskaWtAK8paAomMod6josITPzxs+HhejW0fQZoK+v6RIJsJxsrC774 qnbosnq0+MUhi1P9u6ZBy8AHrIzDyseVjNhT9DX5+TI8dFZLc2b07mjwVJZxkaU= =Ss6c -----END PGP SIGNATURE-----

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

I am not familiar with this expression |Q1> <Q2|. What I really want is Q1*Q2^T, where Q1 and Q2 and vectors defined in my previous email. Actually I think it's easier if I just provide you with the actual equation I need to solve, please see attached. Anna. On 19/10/15 14:15, Lawrence Mitchell wrote:

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1

On 19/10/15 14:00, Anna Kalogirou wrote:

...

Dear Lawrence,

Yes, I have a term where the test function and the basis function are evaluated in separate integrals, this is why I wrote Q*mu*dx. Can I assemble them separately? E.g.

Q1 = assemble(v*dx) Q2 = assemble(mu*dx) The resulting values aren't unknowns (either test or trial functions) so I still don't quite understand.

Q1 is a Co-Function (living in the dual space of V) and Q2 is a Function in V. You can certainly point-wise multiply the values together, but I don't think that's what you want, since that still doesn't give you an operator mapping from the trial to the test space.

Do you want the rank-1 operator formed by:

|Q1> <Q2|

?

Cheers,

Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)

iQEcBAEBAgAGBQJWJO0GAAoJECOc1kQ8PEYvtfIIAOsepbX3x0TCCrJYl3PXr3Vz 2AckYTLSiiUgYYFQ5IPhEOXd4yoDlBtzA9DIo6oyotwa9JltE0KpP45ZK/b0oaAO iAB6s73pcNKkZ328GoiT6/Jp5zu5pCQI+2X2Zhjj6EmYzEpu9J5oAGrLSKHy0nbC Z6gyvG72LFDHcTgNauStl+1yC2D6bTD2JOk9e5ph5ySdLbXL1kegTBt+RHYhfdI3 TsiHLaPZzBp7rRUQ4EGq6pDikS+qTU3IEhKz1onEuWOnHqwHCO14zaRSsiFDBDS9 5ALN9kiqz1rW8XHykWrhZec0MZhxjqZ6akDbCUmdMlwdXje6nYTaYjpQrTgvyko= =Vvc4 -----END PGP SIGNATURE-----

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

-- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/

Dear Lawrence, Thank you very much for the detailed reply. So if I understand correctly, I have to define the forms for A_form, Q1_form, Q2_form, rhs_form and use the bits of code you sent me. If I do that, I get the error "Exception: Linear forms must be defined using test functions only: (1, 'k', 4, 4)" for Q2_form = mu*dx where mu is a trial function. So is it not possible to assemble a vector using the trial function only? Best, Anna. On 20/10/15 12:30, Lawrence Mitchell wrote:

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1

On 19/10/15 19:43, Anna Kalogirou wrote:

...

I am not familiar with this expression |Q1> <Q2|. What I really want is Q1*Q2^T, where Q1 and Q2 and vectors defined in my previous email. Actually I think it's easier if I just provide you with the actual equation I need to solve, please see attached. OK, so we're on the same page.

Note that the operator Q1*Q2^T is dense, so you don't want to form it. However, its action is cheap, since if you have the assembled vectors you just do:

Q1*Q2^T x == Q1 dot(Q2, x)

For the purposes of solving the system you want, I think what you'll want to do is use firedrake to assemble the pieces and then define a matrix-free operator to pass to PETSc (via petsc4py) to actually solve the problem. Finally, you'll want to think about how to precondition the system: B == A + Q1*Q2^T.

Fortunately, the perturbuation is rank-1 and so the Sherman-Morrison-Woodbury formula tells you how to compute a good preconditioner for B, given one for A:

B^{-1} = A^{-1} + (A^{-1} Q1 * Q2^T A^{-1}) / (1 + Q2^T A^{-1} Q1)

To hook this up, you'll need to know some petsc4py, as well as firedrake.

Assume the forms are already defined elsewhere.

First I'll build the operator B:

from firedrake.petsc import PETSc

class MatrixFreeB(object): def __init__(self, A, Q1, Q2): self.A = A self.Q1 = Q1 self.Q2 = Q2

def mult(self, mat, x, y): "Compute y = B*x" # y = Ax self.A.mult(x, y) # alpha = Q2^T x alpha = self.Q2.dot(x) # y = alpha Q1 + y y.axpy(alpha, self.Q1)

# Get reference to PETSc Matrix A = assemble(A_form).M.handle Q1f = assemble(Q1_form) Q2f = assemble(Q2_form)

# Get reference to PETSc Vectors for Q1 and Q2 # I take a copy here, so if Q1 and Q2 change, # you'll need to do something difference with Q1f.dat.vec_ro as v: Q1 = v.copy() A^{- with Q2f.dat.vec_ro as v: Q2 = v.copy()

B = PETSc.Mat().create()

B.setSizes(*A.getSizes()) B.setType(B.Type.PYTHON) # Indicate that B should use the MatrixFreeB class to # compute mult. B.setPythonContext(MatrixFreeB(A, Q1, Q2)) B.setUp()

# For small problems, we don't need a preconditioner at all, so let's # check this works:

solver = PETSc.KSP().create()

solver.setOperators(B)

opts = PETSc.Options() opts["pc_type"] = "none"

solver.setUp() solver.setFromOptions()

# Now let's solve the system.

rhs = assemble(rhs_form) solution = Function(result_space)

with rhs.dat.vec_ro as b: with solution.dat.vec as x: solver.solve(b, x)

Let's try and get this working on a small problem first and then move on to constructing the preconditioner as well. Note that if you have boundary conditions in your operator, we'll have to do a little more work (but not much).

Cheers,

Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)

iQEcBAEBAgAGBQJWJiXEAAoJECOc1kQ8PEYvSGsH/2mDSfh9VkHb+JwJpMt6IEhs +Wx3SfqyfaztbEKay1QPfd50qU7Ojfx1Gv2j7bKFxNh1ppHBqRWo/awZmgbsm2wI aPbVxUL9v5Pj6294DiWVoB8MKFxdQvDcMxm5FNXKRof9M8aLHOOFuAX2fK+aTae7 5a4u/a7yduOh/2qZ19tylaCvDlDx1pBMMU8/mGI28ecmsBaGsBmLR5ObvOnn/LPH HSe/uvBPwautxD5VssgmveF6C5Nwqa/LkYRh5a8hKZhAE3r1PqQhVCiyLbFtEw34 +RdQhldLKKwGenFhF5fuii5v1F4hzctfrotoCkAojuf75UIFXlmdsOKD3eP7lPM= =fus0 -----END PGP SIGNATURE-----

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

-- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/

Dear Lawrence, I believe the code (attached) sets up the matrices/vectrors/operators correctly. However, I have one question: Due to the use of a heavyside function, the operator B on the LHS is nonzero in a block NxN, say, in the bottom right corner. For that reason, when I solved the linear system in Matlab I only considered these nonzero values in the block, because otherwise the LHS would not be invertible. How come and it doesn't complain here? Best, Anna. On 20/10/15 17:03, Lawrence Mitchell wrote:

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1

On 20/10/15 16:05, Anna Kalogirou wrote:

...

Dear Lawrence,

Thank you very much for the detailed reply. So if I understand correctly, I have to define the forms for A_form, Q1_form, Q2_form, rhs_form and use the bits of code you sent me. I hope so, note that I just typed it into the email so there are probably many other issues other than the one below.

...

If I do that, I get the error "Exception: Linear forms must be defined using test functions only: (1, 'k', 4, 4)" for Q2_form = mu*dx

where mu is a trial function. So is it not possible to assemble a vector using the trial function only? I think you can just use a test function from the appropriate space instead (the result will turn out the same I believe).

Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)

iQEcBAEBAgAGBQJWJmW7AAoJECOc1kQ8PEYvR7oH/Rhb9exp1TpOiDOH5QMO2EZT WyYexbWTUrPok3pA8PFvtg3kefpr9ao7j1UZpt9kiXpWCz8V4nL29tmPOTwCnXrp A3OzjCf1N1guwgk48rLatIH8bOuYvOP7yABA9LEOO5Gjzs0it4BdYozBcs0/of47 6Ue4QuKYFkJekgTgU27A8f3ubxJQtuWLaCM3vYaCLh2ecQHNxAhROk2H5zWChz4G E6HblBCY9EvvTbGLGOSUCmBfSra5DXh5pRQSBy+57c0rCM1f/ZYaZp9uxkdUIr4M rS29MAD5VCR6vD4rWDP3uM6qhpzCscPAc6bczhhWDFqhG0TIc2dp6z6o1Rb5RiA= =DF/m -----END PGP SIGNATURE-----

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

-- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/

But the matrix A has the same form as C, namely it has nonzero variables in a block submatrix only. On 21/10/15 14:26, Lawrence Mitchell wrote:

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1

On 21/10/15 14:14, Anna Kalogirou wrote:

...

Dear Lawrence,

I believe the code (attached) sets up the matrices/vectrors/operators correctly. However, I have one question: Due to the use of a heavyside function, the operator B on the LHS is nonzero in a block NxN, say, in the bottom right corner. For that reason, when I solved the linear system in Matlab I only considered these nonzero values in the block, because otherwise the LHS would not be invertible. How come and it doesn't complain here? So it's clear that the matrix C = Q1*Q2^T is not invertible, however, presumably you're solving for:

B = A + C

where A is full rank and invertible. So I don't immediately see why there would be a problem in trying to invert B which is what the solver does.

Lawrence

-----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)

iQEcBAEBAgAGBQJWJ5KFAAoJECOc1kQ8PEYv7F8IANoAlaJs8bPi72Hi/YdveiG+ y02EBN/gzM8JQtJAV/4IPyldt//koM4BEp6BAUG0T4GlZrXPoJY8wz9lG9SJnVer XARcV8C7EmQV5uW9vgK2pGMbt7ZHrIEClVZCAaGYT/rmP+l0XKolMxDXH7tKipIq enjIEoS1mrvVj2iSj58T3NZLM3547r9ChrdS/+Mdkn15rXoR2KS0bapEHITkVppi 6qFSU8OVRbcx/HA8GypX7K2ya1X3MTSwmBhdTO0SAdsOnTgPPBb11y9ziQ3onWPw hqaGpasc3ACSB06VDisGlkZYFNg/iS3nbF3mIsqP/2+8VPaoyxkQmnngKLgGegw= =0uCV -----END PGP SIGNATURE-----

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

-- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/

Hi Lawrence, In the email below you explained to me how to define a matrix free operator B. Is there a way to inspect the values of the matrix behind this operator? Best, Anna. On 20/10/15 12:30, Lawrence Mitchell wrote:

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1

On 19/10/15 19:43, Anna Kalogirou wrote:

...

I am not familiar with this expression |Q1> <Q2|. What I really want is Q1*Q2^T, where Q1 and Q2 and vectors defined in my previous email. Actually I think it's easier if I just provide you with the actual equation I need to solve, please see attached. OK, so we're on the same page.

Note that the operator Q1*Q2^T is dense, so you don't want to form it. However, its action is cheap, since if you have the assembled vectors you just do:

Q1*Q2^T x == Q1 dot(Q2, x)

For the purposes of solving the system you want, I think what you'll want to do is use firedrake to assemble the pieces and then define a matrix-free operator to pass to PETSc (via petsc4py) to actually solve the problem. Finally, you'll want to think about how to precondition the system: B == A + Q1*Q2^T.

Fortunately, the perturbuation is rank-1 and so the Sherman-Morrison-Woodbury formula tells you how to compute a good preconditioner for B, given one for A:

B^{-1} = A^{-1} + (A^{-1} Q1 * Q2^T A^{-1}) / (1 + Q2^T A^{-1} Q1)

To hook this up, you'll need to know some petsc4py, as well as firedrake.

Assume the forms are already defined elsewhere.

First I'll build the operator B:

from firedrake.petsc import PETSc

class MatrixFreeB(object): def __init__(self, A, Q1, Q2): self.A = A self.Q1 = Q1 self.Q2 = Q2

def mult(self, mat, x, y): "Compute y = B*x" # y = Ax self.A.mult(x, y) # alpha = Q2^T x alpha = self.Q2.dot(x) # y = alpha Q1 + y y.axpy(alpha, self.Q1)

# Get reference to PETSc Matrix A = assemble(A_form).M.handle Q1f = assemble(Q1_form) Q2f = assemble(Q2_form)

# Get reference to PETSc Vectors for Q1 and Q2 # I take a copy here, so if Q1 and Q2 change, # you'll need to do something difference with Q1f.dat.vec_ro as v: Q1 = v.copy() A^{- with Q2f.dat.vec_ro as v: Q2 = v.copy()

B = PETSc.Mat().create()

B.setSizes(*A.getSizes()) B.setType(B.Type.PYTHON) # Indicate that B should use the MatrixFreeB class to # compute mult. B.setPythonContext(MatrixFreeB(A, Q1, Q2)) B.setUp()

# For small problems, we don't need a preconditioner at all, so let's # check this works:

solver = PETSc.KSP().create()

solver.setOperators(B)

opts = PETSc.Options() opts["pc_type"] = "none"

solver.setUp() solver.setFromOptions()

# Now let's solve the system.

rhs = assemble(rhs_form) solution = Function(result_space)

with rhs.dat.vec_ro as b: with solution.dat.vec as x: solver.solve(b, x)

Let's try and get this working on a small problem first and then move on to constructing the preconditioner as well. Note that if you have boundary conditions in your operator, we'll have to do a little more work (but not much).

Cheers,

Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)

iQEcBAEBAgAGBQJWJiXEAAoJECOc1kQ8PEYvSGsH/2mDSfh9VkHb+JwJpMt6IEhs +Wx3SfqyfaztbEKay1QPfd50qU7Ojfx1Gv2j7bKFxNh1ppHBqRWo/awZmgbsm2wI aPbVxUL9v5Pj6294DiWVoB8MKFxdQvDcMxm5FNXKRof9M8aLHOOFuAX2fK+aTae7 5a4u/a7yduOh/2qZ19tylaCvDlDx1pBMMU8/mGI28ecmsBaGsBmLR5ObvOnn/LPH HSe/uvBPwautxD5VssgmveF6C5Nwqa/LkYRh5a8hKZhAE3r1PqQhVCiyLbFtEw34 +RdQhldLKKwGenFhF5fuii5v1F4hzctfrotoCkAojuf75UIFXlmdsOKD3eP7lPM= =fus0 -----END PGP SIGNATURE-----

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake

-- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/

Dear Lawrence, I have the following term in my RHS form (all quantities in there are known functions) Hb*step_b*inner(grad(0.5*dt*g*eta0-phi0),grad(v))*dx which is similar to the "local" stiffness matrix Ab. I want to replace this with something of the following form: Mb*Minv*A*(0.5*dt*g*eta0_k-phi0_k) where eta0_k, phi0_k are the coefficients in the approximations for eta0, phi0, respectively. Is this possible? Thanks, Anna. On 12/11/15 14:28, Lawrence Mitchell wrote:

Hi Anna,

On 12/11/15 12:38, Anna Kalogirou wrote:

...

Hi Lawrence,

In the email below you explained to me how to define a matrix free operator B. Is there a way to inspect the values of the matrix behind this operator? Assuming that you are not running in parallel, and you do not have two many degrees of freedom in your problem (I guess you want to look at the values for a small system), you can do the following:

ctx = B.getPythonContext()

dense_A = ctx.A[:, :] Q1_array = ctx.Q1.array_r Q2_array = ctx.Q2.array_r

import numpy as np

dense_B = dense_A + np.outer(Q1_array, Q2_array)

print dense_B

Cheers,

Lawrence

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake