Hi Lawrence, Thanks for the suggestions. Indeed, it doesn’t give any segmentation faults when running it from my mac laptop. On my linux machine, I ran firedrake-update, firedrake-clean and make alltest without any problems, but it still complains about those segmentation faults. Do you have any idea what the problem on linux might be? Now regarding the code itself, I fixed a few mistakes and updated it on bitbucket. It is now running but it’s blowing up. I suspect is due to bad preconditioning (or lack of any). How can I set up the preconditioner for the mixed system via setting the correct PETSc options? Thanks, Anna. On 3 Jun 2016, at 10:22, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
Hi Anna,
On 3 Jun 2016, at 09:59, Anna Kalogirou <a.kalogirou@leeds.ac.uk> wrote:
Good morning,
The code can be found here: https://bitbucket.org/annakalog/buoy2d/src/d05fb55ab1b3b6d75486e5c4b9a268933...
Running the buoy-swe.py in this directory fails for me with compilation errors in your expression for mu_bar:
diff --git a/Inequality constraint/parameters.py b/Inequality constraint/parameters.py index b38c0f6..f35b94b 100755 --- a/Inequality constraint/parameters.py +++ b/Inequality constraint/parameters.py @@ -85,7 +85,7 @@ def def_Hb(Hb, x, L, Lp, H0, rho, Mass, a, Nx): def def_lam_mu_bar(lambda_bar, mu_bar, x, L, Lp, H0, g, rho, Mass, a):
lambda_bar.interpolate(Expression("-0.5*(1.0+copysign(1.0,x[0]-Lp))*g*(sqrt(2*M*tan(a)/rho)+tan(a)*(x[0]-L))", Lp=Lp, g=g, M=Mass, a=a, rho=rho, L=L)) - mu_bar.interpolate(Expression("0.5*(1.0-copysign(1.0,x[0]-Lp))*(sqrt(tan(a)*(Lp-x[0]))", Lp=Lp, a=a)) + mu_bar.interpolate(Expression("0.5*(1.0-copysign(1.0,x[0]-Lp))*(sqrt(tan(a)*(Lp-x[0])))", Lp=Lp, a=a))
return (lambda_bar, mu_bar);
Then I get a problem because making the mixed matrix in solver_mu doesn't work right. This seems to fix things for me:
diff --git a/Inequality constraint/solvers.py b/Inequality constraint/solvers.py index 789be96..8987bca 100755 --- a/Inequality constraint/solvers.py +++ b/Inequality constraint/solvers.py @@ -58,22 +58,24 @@ def solver_mu(eta0, phi0, Z0, W0, LM_0_5, mu, v, v1, dt, Hb, lam_bar, mu_bar, g,
# Build the matrix free operator B class MatrixFreeB(object): - def __init__(self, A, Q1, Q2): + def __init__(self, A, Q1, Q2, scale): self.A = A self.Q1 = Q1 self.Q2 = Q2 + self.scale = scale
# Compute y = alpha Q1 + y = Q2^T*Q1 + A*x def mult(self, mat, x, y): self.A.mult(x, y) alpha = self.Q2.dot(x) y.axpy(alpha, self.Q1) + y.scale(self.scale)
# Define MatrixFreeB class operator B B = PETSc.Mat().create() B.setSizes(*A.getSizes()) B.setType(B.Type.PYTHON) - B.setPythonContext(MatrixFreeB(A, rho/Mass*Q, Q)) + B.setPythonContext(MatrixFreeB(A, rho/Mass*Q, Q, (dt/2))) B.setUp() ''' # Print LHS 'matrix' B to inspect its values @@ -84,7 +86,7 @@ def solver_mu(eta0, phi0, Z0, W0, LM_0_5, mu, v, v1, dt, Hb, lam_bar, mu_bar, g, print 'LHS matrix A + Q*Q^T \n------------\n' + str( dense_B ) + ' \n' ''' # Define combined matrix - mixed_matrix = PETSc.Mat().createNest([[(dt/2)*B, (2/dt)*B_mu], [B_mu, C_lam]]) + mixed_matrix = PETSc.Mat().createNest([[B, (2/dt)*B_mu], [B_mu, C_lam]])
# For small problems, we don't need a preconditioner at all mu_solver = PETSc.KSP().create()
Then one of your forms in buoy-swe.py is wrong, since one of the terms contains a test function but the other doesn't (the RHS just before the solvers_SV call). I don't know what this should be though.
After fixing that you'll probably find that the mu_solver now doesn't work, because the RHS that you feed in comes from V but the solver works on the space W.
But I didn't get any of the segfaults you described! Assuming you have run firedrake-update successfully, can you run firedrake-clean and see if that removes those problems?
Cheers,
Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake