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