December 2015 - firedrake - imperial.ac.uk

Can a Function be used in an Expression?
by Justin Chang 08 Dec '15

08 Dec '15

Hi all, As the title states, can a Function solution be used in an Expression? Say I have the following: beta = 0.001 u = Function(Q) alpha = Expression("exp(Beta*solution)",Beta=Betasolution=u) When I try to run this, I get an error saying "ValueError: setting an array element with a sequence". Am I not doing something correctly or is it simply not possible to do include u into an Expression? Thanks, Justin

2 2

Re: [firedrake] Can a Function be used in an Expression?
by Andrew McRae 08 Dec '15

08 Dec '15

Hi Justin, I think the answer is no. I see Expression as a means to execute a simple piece of C code to interpolate into a 'CG' or 'DG' field. Certain variables are always available to you (the coordinate field x[0], x[1], x[2], and the value pi), along with functions from math.h, and you can pass in your own constants, e.g. expr = Expression("cos(2.0*pi*(x[0] - t))", t=t) What is the exact [mathematical] operation you are trying to do? There might be an easy way to do what you want. On 8 December 2015 at 22:39, Justin Chang <jychang48(a)gmail.com> wrote: > sorry, I meant alpha = > Expression("exp(Beta*solution)",Beta=beta,solution=u) > > On Tue, Dec 8, 2015 at 3:38 PM, Justin Chang <jychang48(a)gmail.com> wrote: > >> Hi all, >> >> As the title states, can a Function solution be used in an Expression? >> >> Say I have the following: >> >> beta = 0.001 >> u = Function(Q) >> alpha = Expression("exp(Beta*solution)",Beta=Betasolution=u) >> >> When I try to run this, I get an error saying "ValueError: setting an >> array element with a sequence". Am I not doing something correctly or is it >> simply not possible to do include u into an Expression? >> >> Thanks, >> Justin >> > >

1 0

Re: [firedrake] Preconditioner for Least-squares FEM
by Colin Cotter 08 Dec '15

08 Dec '15

This is a tricky system to precondition in general. Approaches that are known to work well are: 1) Use a block H(div) preconditioner, i.e. precondition by the operator ( 1 - grad(div .) 0 ) ( 0 1 ) This converges fast with direct solvers to obtain the inverse of this operator, but classical iterative methods do not converge well for the top block, and so you need to use Schwartz preconditioners to get a fast parallel implementation. We don't currently have Firedrake support for that. 2) Hybridise the equation. The resulting reduced system behaves well with classical iterative methods or AMG. We have been able to construct the hybridised system using petsc4py. all the best --cjc On 15 October 2015 at 14:18, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: > As usual, the Firedrake person in the best position to answer your > question is Lawrence, but he's currently on holiday. > > I won't try to answer your question, but in an old email > <http://mailman.ic.ac.uk/pipermail/firedrake/2015-July/000953.html>, > Lawrence suggested you try using multigrid on the blocks; e.g. > > 'fieldsplit_0_pc_type': 'gamg', > 'fieldsplit_1_pc_type': 'gamg', > > My linear solver knowledge is pretty dire, so I'll stop talking here :) > > Andrew > > > On 15 October 2015 at 13:58, Justin Chang <jychang48(a)gmail.com> wrote: > >> I did ask them, but it went unanswered for the last 1.5 months... >> >> Let me clarify my question a bit though, I tried using -ksp_type cg >> -pc_type hypre -pc_hypre_type boomeramg, but it is not supported for >> matnest. I added nest=False but it seems to break my solver somehow. >> >> The only thing that works for me is -ksp_type cg -pc_type bjacobi plus >> MixedVectorSpaceBasis(W, ...) but the number of solver iterations increases >> with problem size... >> >> Thanks, >> Justin >> >> On Thu, Oct 15, 2015 at 6:50 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >> wrote: >> >>> I think the PETSc mailing list is a much better place to ask questions >>> like this? >>> >>> On 15 October 2015 at 13:47, Justin Chang <jychang48(a)gmail.com> wrote: >>> >>>> Hi all, >>>> >>>> I am attempting to solve Darcy's equation: >>>> >>>> u + grad[p] = 0 >>>> div[u] = f >>>> >>>> The weak form under the least-squares finite element method (LSFEM) >>>> looks like this: >>>> >>>> (u + grad[p]; v + grad[q]) + div[u]*div[v] = (f; div[v]) >>>> >>>> For H(div) elements like RT0, these options worked nicely: >>>> >>>> -ksp_type gmres >>>> -pc_type fieldsplit >>>> -pc_fieldsplit_type schur >>>> -pc_fieldsplit_schur_precondition selfp >>>> -fieldsplit_0_ksp_type preonly >>>> -fieldsplit_0_pc_type bjacobi >>>> -fieldsplit_0_sub_pc_type ilu >>>> -fieldsplir_1_ksp_type preonly >>>> -fieldsplit_1_pc_type hypre >>>> >>>> but for the above LSFEM, it does not work (get a DIVERGED_MAX_IT >>>> error). I heard multigrid methods are good for these types of problems, so >>>> how do I tweak the above parameters? >>>> >>>> Thanks, >>>> Justin >>>> >>> >>> >>> _______________________________________________ >>> firedrake mailing list >>> firedrake(a)imperial.ac.uk >>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>> >>> >> > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

5 11

Re: [firedrake] implementing these FEniCS functions in Firedrake
by David Ham 08 Dec '15

08 Dec '15

Yes. I think that's right. On Tue, 8 Dec 2015 at 17:54 Justin Chang <jychang48(a)gmail.com> wrote: > Lawrence, > > Okay yeah I got it to work, at least for structured grids. > > Now if I had unstructured grids, element_meas is not immediately known > and varies among elements. In this case, would I use h = CellSize(mesh) and > define element_meas = 1/h/h (I am using quadrilaterals for this)? > > Thanks, > Justin > > On Tue, Dec 8, 2015 at 4:17 AM, Lawrence Mitchell < > lawrence.mitchell(a)imperial.ac.uk> wrote: > >> Hi Justin, >> >> On 08/12/15 09:19, Justin Chang wrote: >> > Hi all, >> > >> > So our research group has this augmented lagrangian solver for the >> > LSFEM in FEniCS. I have attached the code. >> > >> > While the code works and does what we want, it is extremely slow! I >> > suspect it is because of these lines of code: >> > >> > 1) >> > >> > a_LMB, L_LMB = {}, {} >> > element_meas = 0.5 / (XSeed * YSeed) >> > for i in range(mesh.num_cells()): >> > a_LMB[i] = (element_meas) * (1./element_meas**2) * div(v) * div(w) * >> > dx(i) >> > L_LMB[i] = (1./element_meas) * div(w) * dx(i) >> >> So you mark each cell as its own subdomain. Then you build a weak >> form integrating over each cell. >> >> ... >> >> > 3) >> > >> > a, L = aStandard, LStandard >> > for i in range(mesh.num_cells()): >> > a += gamma * a_LMB[i] >> > L += lmbda[i] * L_LMB[i] >> >> Here you build the global weak form: >> >> a = aStandard + \sum_i a_LMB[i] >> >> surely this is just equivalent to saying: >> >> a = a_Standard + gamma * (1/element_meas)*div(v)*div(w)*dx >> >> and similarly: >> >> L = LStandard + (1/element_meas) * lmbda * div(w) * dx >> >> Where lmbda is just a DG0 Function. >> >> > What the FFC just-in-time compiler does is that it compiles each and >> > every single one of these element-wise forms (e.g., if my mesh has >> > over 1000 elements, the JIT compiler will be called 1000 times). This >> > takes a stupidly long time to run. If I wanted to do something like >> > this in Firedrake, how would I go about making this more efficient? >> >> I think you don't have to define any of these tiny forms (unless I'm >> really missing something). Moreover, since they're all symbolic, they >> don't need to be redefined inside the loop over and over either. >> >> >> > For 2) I believe the solution is to do something like >> > "Error.assign(assemble(div(v)*e*dx)/assemble(e*dx))" as suggested in a >> > previous thread. >> >> Yes, it looks like this is just the projection of div(mysol.sub(0)) >> into a DG0 space. >> >> Does that all make sense? >> >> Cheers, >> >> Lawrence >> >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >

1 0

Re: [firedrake] Preconditioner for Least-squares FEM
by Colin Cotter 08 Dec '15

08 Dec '15

I'm not too sure that the Hdiv preconditioner is actually good for the least squares formulation. Worth a try though. all the best --cjc On 8 December 2015 at 17:58, Justin Chang <jychang48(a)gmail.com> wrote: > Hi Colin (or anyone), > > Coming back to this > > 1) Could you elaborate a little more on implementing this block H(div) > operator in firedrake/petsc: > > ( 1 - grad(div .) 0 ) > ( 0 1 ) > > My guess at least for that bottom right entry would be to use the > fieldsplit, where fieldsplit_1_pc would be bjacobi/ilu? > > also > > 2) What exactly do you guys mean by "hybridise the equation"? Like how > would it look in the context of LSFEM? > > Thanks, > Justin > > > On Thu, Oct 15, 2015 at 7:32 AM, Colin Cotter <colin.cotter(a)imperial.ac.uk > > wrote: > >> This is a tricky system to precondition in general. Approaches that are >> known to work well are: >> >> 1) Use a block H(div) preconditioner, i.e. precondition by the operator >> >> ( 1 - grad(div .) 0 ) >> ( 0 1 ) >> >> This converges fast with direct solvers to obtain the inverse of this >> operator, but classical iterative methods do not converge well for the top >> block, and so you need to use Schwartz preconditioners to get a fast >> parallel implementation. We don't currently have Firedrake support for that. >> >> 2) Hybridise the equation. The resulting reduced system behaves well with >> classical iterative methods or AMG. We have been able to construct the >> hybridised system using petsc4py. >> >> all the best >> --cjc >> >> On 15 October 2015 at 14:18, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: >> >>> As usual, the Firedrake person in the best position to answer your >>> question is Lawrence, but he's currently on holiday. >>> >>> I won't try to answer your question, but in an old email >>> <http://mailman.ic.ac.uk/pipermail/firedrake/2015-July/000953.html>, >>> Lawrence suggested you try using multigrid on the blocks; e.g. >>> >>> 'fieldsplit_0_pc_type': 'gamg', >>> 'fieldsplit_1_pc_type': 'gamg', >>> >>> My linear solver knowledge is pretty dire, so I'll stop talking here :) >>> >>> Andrew >>> >>> >>> On 15 October 2015 at 13:58, Justin Chang <jychang48(a)gmail.com> wrote: >>> >>>> I did ask them, but it went unanswered for the last 1.5 months... >>>> >>>> Let me clarify my question a bit though, I tried using -ksp_type cg >>>> -pc_type hypre -pc_hypre_type boomeramg, but it is not supported for >>>> matnest. I added nest=False but it seems to break my solver somehow. >>>> >>>> The only thing that works for me is -ksp_type cg -pc_type bjacobi plus >>>> MixedVectorSpaceBasis(W, ...) but the number of solver iterations increases >>>> with problem size... >>>> >>>> Thanks, >>>> Justin >>>> >>>> On Thu, Oct 15, 2015 at 6:50 AM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> >>>> wrote: >>>> >>>>> I think the PETSc mailing list is a much better place to ask questions >>>>> like this? >>>>> >>>>> On 15 October 2015 at 13:47, Justin Chang <jychang48(a)gmail.com> wrote: >>>>> >>>>>> Hi all, >>>>>> >>>>>> I am attempting to solve Darcy's equation: >>>>>> >>>>>> u + grad[p] = 0 >>>>>> div[u] = f >>>>>> >>>>>> The weak form under the least-squares finite element method (LSFEM) >>>>>> looks like this: >>>>>> >>>>>> (u + grad[p]; v + grad[q]) + div[u]*div[v] = (f; div[v]) >>>>>> >>>>>> For H(div) elements like RT0, these options worked nicely: >>>>>> >>>>>> -ksp_type gmres >>>>>> -pc_type fieldsplit >>>>>> -pc_fieldsplit_type schur >>>>>> -pc_fieldsplit_schur_precondition selfp >>>>>> -fieldsplit_0_ksp_type preonly >>>>>> -fieldsplit_0_pc_type bjacobi >>>>>> -fieldsplit_0_sub_pc_type ilu >>>>>> -fieldsplir_1_ksp_type preonly >>>>>> -fieldsplit_1_pc_type hypre >>>>>> >>>>>> but for the above LSFEM, it does not work (get a DIVERGED_MAX_IT >>>>>> error). I heard multigrid methods are good for these types of problems, so >>>>>> how do I tweak the above parameters? >>>>>> >>>>>> Thanks, >>>>>> Justin >>>>>> >>>>> >>>>> >>>>> _______________________________________________ >>>>> firedrake mailing list >>>>> firedrake(a)imperial.ac.uk >>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >>>>> >>>>> >>>> >>> >> >> >> -- >> http://www.imperial.ac.uk/people/colin.cotter >> >> www.cambridge.org/9781107663916 >> >> >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

implementing these FEniCS functions in Firedrake
by Justin Chang 08 Dec '15

08 Dec '15

Hi all, So our research group has this augmented lagrangian solver for the LSFEM in FEniCS. I have attached the code. While the code works and does what we want, it is extremely slow! I suspect it is because of these lines of code: 1) a_LMB, L_LMB = {}, {} element_meas = 0.5 / (XSeed * YSeed) for i in range(mesh.num_cells()): a_LMB[i] = (element_meas) * (1./element_meas**2) * div(v) * div(w) * dx(i) L_LMB[i] = (1./element_meas) * div(w) * dx(i) 2) ErrorSpace = FunctionSpace(mesh,"DG",0) Error = Function(ErrorSpace) ... for i in range(mesh.num_cells()): Error.vector()[i] = (1./element_meas) * assemble(div(mysol.sub(0))*dx(i)) 3) a, L = aStandard, LStandard for i in range(mesh.num_cells()): a += gamma * a_LMB[i] L += lmbda[i] * L_LMB[i] What the FFC just-in-time compiler does is that it compiles each and every single one of these element-wise forms (e.g., if my mesh has over 1000 elements, the JIT compiler will be called 1000 times). This takes a stupidly long time to run. If I wanted to do something like this in Firedrake, how would I go about making this more efficient? For 2) I believe the solution is to do something like "Error.assign(assemble(div(v)*e*dx)/assemble(e*dx))" as suggested in a previous thread. How would I more efficiently implement 1) and 3) in Firedrake? Thanks, Justin

2 2

Re: [firedrake] FD on mac
by David Ham 05 Dec '15

05 Dec '15

On Mon, 30 Nov 2015 at 18:26 Mitchell, Lawrence < lawrence.mitchell(a)imperial.ac.uk> wrote: > > > On 30 Nov 2015, at 18:01, David Ham <David.Ham(a)imperial.ac.uk> wrote: > > > > This is quite possible. Unfortunately we have no access to Mavericks > machines which makes it difficult to work fix this problem. > > > > The script does currently work for me on El Capitan! > > The proximate cause of Onno's problem appears to be the lack of fortran > compiler. In the past, this came via the brew package "gfortran" is this > still the case? Those of us for whom everything works have had it all in > place for a while, no? > gfortran is currently provided by the plain gcc package. ls -l $(which gfortran) lrwxr-xr-x 1 dham admin 32 23 Jul 12:44 /usr/local/bin/gfortran -> ../Cellar/gcc/5.2.0/bin/gfortran > > Lawrence > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

5 23

Re: [firedrake] Help visualization
by Andrew McRae 03 Dec '15

03 Dec '15

Look for a box that says 'Solid Color' and change that to 'function_5' or similar. Best, Andrew On 3 December 2015 at 16:34, Angwenyi David <kipkoej(a)gmail.com> wrote: > The visual output of the helmholtz, burgers and mixed poisson in the > tutorials are rectangular surfaces. This is when I open their respective > .vtu files. What could be the problem? > > > On 03 Dec 2015, at 17:07, Lawrence Mitchell < > lawrence.mitchell(a)imperial.ac.uk> wrote: > > > > > > > > On 03/12/15 16:05, Angwenyi David wrote: > >> I am actually trying to visualize the output in the tutorials. We > already have for example: File(“poisson_mixed.pvd”) << u > >> I have done: output = File(“poisson_mixed.pvd”) > >> output << u > >> del output > >> > >> And after that, I cannot locate the file poisson_mixed.pvd > > > > This file is created in the directory you started the python process > > in. You can find that from within the python process with: > > > > import os > > > > print os.path.abspath(os.path.curdir) > > > > Cheers, > > > > Lawrence > > > > _______________________________________________ > > firedrake mailing list > > firedrake(a)imperial.ac.uk > > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > _______________________________________________ > firedrake mailing list > firedrake(a)imperial.ac.uk > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >

2 1

Help visualization
by Angwenyi David 03 Dec '15

03 Dec '15

So far so good. However, I need help on visualization of results in paraview. I have it installed and running but I don't know on how to visualize .pvd files using paraview. Thanks in advance. David

2 8