I am solving elastodynamic+elastostatic equations for waves and fracture using a spectral element method in 2D.

One of the steps involves solving a large system of equations. I am currently using a preconditioned conjugate gradient with jacobi preconditioner. This step is very slow and does not converge for higher mesh resolution. I tried looking into Petsc but I am not sure if Petsc.jl works with julia 1.0.

I am trying to use iterativesolvers.jl to just play around with simple conjugate gradient method for a smaller problem. The issue is that my code has operators defined to construct the vectors from element level matrices and therefore I do not have a ‘matrix-form’ for `A`

in the equation `Ax = b`

.

The problem is defined as:

`K1u1 = K2u2`

I need to solve for `u1`

, given `K1, K2, u2`

. Right now, I use an initial guess for `u1`

and assemble `K1u1`

as well as `K2u2`

from their local parts defined for each element in the finite element mesh. Given the function to assemble these vectors, how would I solve this using iterativesolvers? I have read that iterativesolvers gives the option to define custom operators, but I couldn’t find an example anywhere.

A second question is whether I can use the standard Petsc library written in c and call that from julia?

Any other input about solving this is also appreciated. @mohamed82008