search

firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
508 stars 159 forks source link

Solving a PDE in complex branch using hybridisation PC fails #1279

Open colinjcotter opened 6 years ago

colinjcotter commented 6 years ago

The hybridisation preconditioner produces an error message on the complex branch

I think this is because the surface terms in the PC need to be changed to use inner instead of *.

from firedrake import *

n = 20
mesh = UnitSquareMesh(n,n)

V = FunctionSpace(mesh,"BDM",1)
Q = FunctionSpace(mesh, "DG", 0)

M = MixedFunctionSpace((V,Q))

w, q = TestFunctions(M)
u, h = TrialFunctions(M)

perp = lambda u: as_vector([-u[1], u[0]])
dt = Constant(1.0)
g = Constant(1.0)
f = Constant(1.0)
H = Constant(1.0)
#alpha = Constant(1+1j)                                                                                        
alpha = Constant(1+1j)

a = (
    inner(alpha*u,w)
    - dt*f*inner(perp(u), w)
    + dt*inner(h,div(w))
    + inner(alpha*h,q)
    - dt*H*inner(div(u),q)
)*dx

assemble(a)

x, y = SpatialCoordinate(mesh)
f = Function(Q).interpolate(cos(2*pi*x))
L = inner(f,q)*dx

assemble(L)

m0 = Function(M)

lu_parameters = {'mat_type':'aij',
                 'ksp_type':'preonly',
                 'pc_type':'lu',
                 'pc_factor_mat_solver_type': 'mumps'}

hyb_parameters = {'mat_type': 'matfree',
                  'ksp_type': 'preonly',
                  'pc_type': 'python',
        'pc_python_type': 'firedrake.HybridizationPC',
                  'hybridization': {'ksp_type': 'preonly',
                                    'pc_type': 'lu',
                                    'pc_factor_mat_solver_type': 'mumps'}}

solve(a==L, m0, solver_parameters=lu_parameters)

u0 = Function(V)
q0 = Function(Q)

u1, q1 = m0.split()
u1.assign(u0)
q1.assign(q0)
File('lu.pvd').write(u1,q1)

solve(a==L, m0, solver_parameters=hyb_parameters)
u0 = Function(V)
q0 = Function(Q)

u1, q1 = m0.split()
u1.assign(u0)
q1.assign(q0)
File('hyb.pvd').write(u1,q1)
thomasgibson commented 6 years ago

Can you be a bit more specific with the error message? Is it failing at the UFL level? It might have something to do (as you said) with the surface terms the preconditioner introduces:

Kform = gammar('+') * ufl.jump(sigma, n=n) * ufl.dS

I need to look at the current complex branch and documentation to see how this should be rewritten. Would inner(gammar('+'), jump(sigma, n=n))*dS work? Or do I need to always conj the second argument?

rckirby commented 6 years ago

Inner inserts a conj, dot and * do not

Sent from my iPhone

On Sep 4, 2018, at 10:38 AM, Thomas H. Gibson notifications@github.com wrote:

Can you be a bit more specific with the error message? Is it failing at the UFL level? It might have something to do (as you said) with the surface terms the preconditioner introduces:

Kform = gammar('+') ufl.jump(sigma, n=n) ufl.dS I need to look at the current complex branch and documentation to see how this should be rewritten. Would inner(gammar('+'), jump(sigma, n=n))*dS work? Or do I need to always conj the second argument?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

colinjcotter commented 6 years ago

I pushed a fix for the forms to the complex-hybridised branch. Now I get compile errors from Slate:

/home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:67:0: warning: ignoring #pragma coffee expression [-Wunknown-pragmas]

pragma coffee expression

^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:115:0: warning: ignoring #pragma coffee expression [-Wunknown-pragmas]

pragma coffee expression

^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:169:0: warning: ignoring #pragma coffee expression [-Wunknown-pragmas]

pragma coffee expression

^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:204:0: warning: ignoring #pragma coffee expression [-Wunknown-pragmas]

pragma coffee expression

^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:287:0: warning: ignoring #pragma coffee expression [-Wunknown-pragmas]

pragma coffee expression

^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:13:130: error: expected ‘,’ or ‘...’ before ‘’ token static inline void subkernel0_cell_to_00_cell_integralotherwise (const Eigen::MatrixBase & A , co nst double complex restrict coords , const double complex restrict w_0 , const double complex restrict w_1 , const double complex restrict w_2 , const double complex restrict w_3 )

               ^

/home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp: In function ‘v oid subkernel0_cell_to_00_cell_integral_otherwise(const Eigen::MatrixBase&, double)’: /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:18:21: error: e xpected initializer before ‘t0’ double complex t0 = ((w_0[0]) (w_1[0])); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:19:21: error: e xpected initializer before ‘t1’ double complex t1 = ((-1) (coords[0])); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:20:21: error: e xpected initializer before ‘t2’ double complex t2 = ((t1) + (coords[2])); /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:21:21: error: e xpected initializer before ‘t3’ double complex t3 = ((-1) (coords[1])); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:22:21: error: e xpected initializer before ‘t4’ double complex t4 = ((t3) + (coords[5])); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:23:21: error: e xpected initializer before ‘t5’ double complex t5 = ((t1) + (coords[4])); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:24:21: error: e xpected initializer before ‘t6’ double complex t6 = ((t3) + (coords[3])); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:25:21: error: expected initializer before ‘t7’ double complex t7 = ((((t2) (t4))) + (((-1) (((t5) (t6)))))); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:26:21: error: expected initializer before ‘t8’ double complex t8 = cabs(t7); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:27:34: error: expected initializer before ‘t9’ static const double complex t9[4] = {0.166666666666667, 0.166666666666667, 0.166666666666667}; ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:28:21: error: expected initializer before ‘t10’ double complex t10 = ((1) / (t7)); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:29:21: error: expected initializer before ‘t11’ double complex t11 = ((t5) (t10)); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:30:21: error: expected initializer before ‘t12’ double complex t12 = ((t2) (t10)); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:31:21: error: expected initializer before ‘t13’ double complex t13 = ((t4) * (t10)); ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:32:34: error: expected initializer before ‘t14’ static const double complex t14[3][6] = {{-0.166666666666667, 0.333333333333333, 0.166666666666667, -0.333333333333333, -1.16666666666667, 0.333333333333334}, ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:410:1: error: expected ‘}’ at end of input } ^ /home/cjc1/firedrake-complex/firedrake/.cache/pyop2/f11650969a16677cd0d7edbd038edf17_p18751.cpp:410:1: error: expected ‘}’ at end of input

thomasgibson commented 6 years ago

Looks like something funny is happening in the compiler on the complex branch. Can you send me the files:

colinjcotter commented 6 years ago

Hi Thomas,

Probably for us to make progress on this you will need to build a complex version of Firedrake - please can you do this? You can do this in a separate directory and different virtual environment. The trick is to get the install script from the complex branch, and then use the invocation in the Jenkins file on that branch, then finally change to the complex branch after installation.

all the best

--Colin

From: Thomas H. Gibson notifications@github.com Sent: 05 September 2018 12:37:20 To: firedrakeproject/firedrake Cc: Cotter, Colin J; Author Subject: Re: [firedrakeproject/firedrake] Solving a PDE in complex branch using hybridisation PC fails (#1279)

Looks like something funny is happening in the compiler on the complex branch. Can you send me the files:

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/firedrakeproject/firedrake/issues/1279#issuecomment-418699042, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AF0q8shN5wpLSp00-z8z_UlwLixwBrh4ks5uX7dwgaJpZM4WY6ng.

miklos1 commented 6 years ago

My guess just by looking at the error message: SLATE triggers code to be compiled as C++, but C++ does not understand double complex, even though C does. See: https://en.wikipedia.org/wiki/Compatibility_of_C_and_C++ and grep for "Complex arithmetic."