- firedrake - imperial.ac.uk

NonlinearVariationalSolver problem
by Shipton, Jemma 15 Sep '15

15 Sep '15

Hello Firedrakers! I am attempting to setup a NonlinearVariationalSolver and I am encountering some difficulties which I hope someone can help with... The problem appears to be with the form: F = ( phi*div(v) + inner(w,v) - div(theta*w)*cp*pi + g*w[1] )*dx + theta*cp*inner(w,n)*ds_b I get the message: ufl.log.UFLException: Expecting all integrals in a form to share geometric dimension, got (). The code works if I remove both the g*w[1] and the boundary term. These terms are fine in a LinearVariationalSolver. Example code is attached... and the full error is: Expecting all integrals in a form to share geometric dimension, got (). Traceback (most recent call last): File "nonlinear_hydrostatic_solver.py", line 111, in <module> HSolver.solve() File "<string>", line 2, in solve File "/home/jemma/src/PyOP2/pyop2/profiling.py", line 203, in wrapper return f(*args, **kwargs) File "/home/jemma/src/firedrake/firedrake/variational_solver.py", line 173, in solve self.snes.solve(None, v) File "PETSc/SNES.pyx", line 520, in petsc4py.PETSc.SNES.solve (src/petsc4py.PETSc.c:167382) File "PETSc/petscsnes.pxi", line 286, in petsc4py.PETSc.SNES_Jacobian (src/petsc4py.PETSc.c:33352) File "/home/jemma/src/firedrake/firedrake/solving_utils.py", line 230, in form_jacobian nest=problem._nest) File "/home/jemma/src/firedrake/firedrake/assemble.py", line 68, in assemble inverse=inverse, nest=nest) File "/home/jemma/src/firedrake/firedrake/assemble.py", line 101, in _assemble inverse=inverse) File "/home/jemma/src/firedrake/firedrake/ffc_interface.py", line 267, in compile_form fd = compute_form_data(form) File "/usr/local/lib/python2.7/dist-packages/ufl/algorithms/compute_form_data.py", line 188, in compute_form_data "Expecting all integrals in a form to share geometric dimension, got %s." % str(tuple(sorted(gdims)))) File "/usr/local/lib/python2.7/dist-packages/ufl/assertions.py", line 37, in ufl_assert if not condition: error(*message) File "/usr/local/lib/python2.7/dist-packages/ufl/log.py", line 151, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Expecting all integrals in a form to share geometric dimension, got (). Thanks in advance for any help! Jemma

2 1

Re: [firedrake] results from interpolate not as expected
by Colin Cotter 14 Sep '15

14 Sep '15

Thanks. Remarkably, this misunderstanding of the node ordering didn't affect any of my limiter kernels (mostly as they map to a Taylor basis, which I made myself). cheers --cjc On 14 September 2015 at 20:26, McRae, Andrew <a.mcrae12(a)imperial.ac.uk> wrote: > Yes. This should work: > > remap_kernel =""" > u_out[0][0] = u_in[0][0]; > u_out[1][0] = u_in[1][0]; > u_out[2][0] = u_in[3][0]; > u_out[3][0] = u_in[4][0]; > """ > > On 14 September 2015 at 20:23, Cotter, Colin J < > colin.cotter(a)imperial.ac.uk> wrote: > >> You're talking about the ordering of the vertical? >> >> --cjc >> >> On 14 September 2015 at 20:19, McRae, Andrew <a.mcrae12(a)imperial.ac.uk> >> wrote: >> >>> Ah, my intuition for the local order of nodes was wrong. >>> >>> For example, the following code >>> --- >>> from firedrake import * >>> >>> mesh = UnitIntervalMesh(4) >>> V = FunctionSpace(mesh, "CG", 4) >>> u = Function(V).interpolate(Expression("x[0]")) >>> >>> print_kernel =""" >>> printf("%f %f %f %f %f\\n", u[0][0], u[1][0], u[2][0], u[3][0], u[4][0]); >>> """ >>> >>> par_loop(print_kernel, dx, {"u": (u, RW)}) >>> print u.dat.data >>> --- >>> produces >>> 0.250000 0.000000 0.187500 0.125000 0.062500 >>> 0.500000 0.250000 0.437500 0.375000 0.312500 >>> 0.750000 0.500000 0.687500 0.625000 0.562500 >>> 1.000000 0.750000 0.937500 0.875000 0.812500 >>> >>> (our IntervalMesh cells are 'backwards', just to make things more fun, >>> but this shouldn't be the case in extruded meshes) >>> >>> In other words, it looks like the CG 'ends' come before the middle, so >>> try 0, 1, 3, 4. >>> >>> Thanks for spotting this, because it means some of my own visualisation >>> will be wrong! >>> >>> On 14 September 2015 at 19:57, Cotter, Colin J < >>> colin.cotter(a)imperial.ac.uk> wrote: >>> >>>> Dear friends, >>>> I have been doing some stuff with a DG1 x CG2 field, as we use it for >>>> temperature in the compatible FEM framework. I'm getting some funny >>>> output when trying to interpolate. In the following code: >>>> >>>> https://gist.github.com/colinjcotter/96e23745d40a73303b05 >>>> >>>> I interpolate a function to a DG1 x CG2 field. The function represented >>>> should >>>> be continuous. >>>> >>>> Then, I remap it using a Firedrake kernel to a DG1 x DG1 field. The >>>> function >>>> should still be continuous. >>>> >>>> When I visualise the field in Paraview, I see that it has jumps in the >>>> vertical >>>> direction. >>>> >>>> What went wrong? >>>> >>>> cheers >>>> --cjc >>>> >>>> -- >>>> http://www.imperial.ac.uk/people/colin.cotter >>>> >>>> www.cambridge.org/9781107663916 >>>> >>>> >>>> >>> >> >> >> -- >> http://www.imperial.ac.uk/people/colin.cotter >> >> www.cambridge.org/9781107663916 >> >> >> > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Re: [firedrake] results from interpolate not as expected
by Andrew McRae 14 Sep '15

14 Sep '15

Yes. This should work: remap_kernel =""" u_out[0][0] = u_in[0][0]; u_out[1][0] = u_in[1][0]; u_out[2][0] = u_in[3][0]; u_out[3][0] = u_in[4][0]; """ On 14 September 2015 at 20:23, Cotter, Colin J <colin.cotter(a)imperial.ac.uk> wrote: > You're talking about the ordering of the vertical? > > --cjc > > On 14 September 2015 at 20:19, McRae, Andrew <a.mcrae12(a)imperial.ac.uk> > wrote: > >> Ah, my intuition for the local order of nodes was wrong. >> >> For example, the following code >> --- >> from firedrake import * >> >> mesh = UnitIntervalMesh(4) >> V = FunctionSpace(mesh, "CG", 4) >> u = Function(V).interpolate(Expression("x[0]")) >> >> print_kernel =""" >> printf("%f %f %f %f %f\\n", u[0][0], u[1][0], u[2][0], u[3][0], u[4][0]); >> """ >> >> par_loop(print_kernel, dx, {"u": (u, RW)}) >> print u.dat.data >> --- >> produces >> 0.250000 0.000000 0.187500 0.125000 0.062500 >> 0.500000 0.250000 0.437500 0.375000 0.312500 >> 0.750000 0.500000 0.687500 0.625000 0.562500 >> 1.000000 0.750000 0.937500 0.875000 0.812500 >> >> (our IntervalMesh cells are 'backwards', just to make things more fun, >> but this shouldn't be the case in extruded meshes) >> >> In other words, it looks like the CG 'ends' come before the middle, so >> try 0, 1, 3, 4. >> >> Thanks for spotting this, because it means some of my own visualisation >> will be wrong! >> >> On 14 September 2015 at 19:57, Cotter, Colin J < >> colin.cotter(a)imperial.ac.uk> wrote: >> >>> Dear friends, >>> I have been doing some stuff with a DG1 x CG2 field, as we use it for >>> temperature in the compatible FEM framework. I'm getting some funny >>> output when trying to interpolate. In the following code: >>> >>> https://gist.github.com/colinjcotter/96e23745d40a73303b05 >>> >>> I interpolate a function to a DG1 x CG2 field. The function represented >>> should >>> be continuous. >>> >>> Then, I remap it using a Firedrake kernel to a DG1 x DG1 field. The >>> function >>> should still be continuous. >>> >>> When I visualise the field in Paraview, I see that it has jumps in the >>> vertical >>> direction. >>> >>> What went wrong? >>> >>> cheers >>> --cjc >>> >>> -- >>> http://www.imperial.ac.uk/people/colin.cotter >>> >>> www.cambridge.org/9781107663916 >>> >>> >>> >> > > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > >

1 0

Re: [firedrake] results from interpolate not as expected
by Colin Cotter 14 Sep '15

14 Sep '15

You're talking about the ordering of the vertical? --cjc On 14 September 2015 at 20:19, McRae, Andrew <a.mcrae12(a)imperial.ac.uk> wrote: > Ah, my intuition for the local order of nodes was wrong. > > For example, the following code > --- > from firedrake import * > > mesh = UnitIntervalMesh(4) > V = FunctionSpace(mesh, "CG", 4) > u = Function(V).interpolate(Expression("x[0]")) > > print_kernel =""" > printf("%f %f %f %f %f\\n", u[0][0], u[1][0], u[2][0], u[3][0], u[4][0]); > """ > > par_loop(print_kernel, dx, {"u": (u, RW)}) > print u.dat.data > --- > produces > 0.250000 0.000000 0.187500 0.125000 0.062500 > 0.500000 0.250000 0.437500 0.375000 0.312500 > 0.750000 0.500000 0.687500 0.625000 0.562500 > 1.000000 0.750000 0.937500 0.875000 0.812500 > > (our IntervalMesh cells are 'backwards', just to make things more fun, but > this shouldn't be the case in extruded meshes) > > In other words, it looks like the CG 'ends' come before the middle, so try > 0, 1, 3, 4. > > Thanks for spotting this, because it means some of my own visualisation > will be wrong! > > On 14 September 2015 at 19:57, Cotter, Colin J < > colin.cotter(a)imperial.ac.uk> wrote: > >> Dear friends, >> I have been doing some stuff with a DG1 x CG2 field, as we use it for >> temperature in the compatible FEM framework. I'm getting some funny >> output when trying to interpolate. In the following code: >> >> https://gist.github.com/colinjcotter/96e23745d40a73303b05 >> >> I interpolate a function to a DG1 x CG2 field. The function represented >> should >> be continuous. >> >> Then, I remap it using a Firedrake kernel to a DG1 x DG1 field. The >> function >> should still be continuous. >> >> When I visualise the field in Paraview, I see that it has jumps in the >> vertical >> direction. >> >> What went wrong? >> >> cheers >> --cjc >> >> -- >> http://www.imperial.ac.uk/people/colin.cotter >> >> www.cambridge.org/9781107663916 >> >> >> > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

Re: [firedrake] results from interpolate not as expected
by Andrew McRae 14 Sep '15

14 Sep '15

Ah, my intuition for the local order of nodes was wrong. For example, the following code --- from firedrake import * mesh = UnitIntervalMesh(4) V = FunctionSpace(mesh, "CG", 4) u = Function(V).interpolate(Expression("x[0]")) print_kernel =""" printf("%f %f %f %f %f\\n", u[0][0], u[1][0], u[2][0], u[3][0], u[4][0]); """ par_loop(print_kernel, dx, {"u": (u, RW)}) print u.dat.data --- produces 0.250000 0.000000 0.187500 0.125000 0.062500 0.500000 0.250000 0.437500 0.375000 0.312500 0.750000 0.500000 0.687500 0.625000 0.562500 1.000000 0.750000 0.937500 0.875000 0.812500 (our IntervalMesh cells are 'backwards', just to make things more fun, but this shouldn't be the case in extruded meshes) In other words, it looks like the CG 'ends' come before the middle, so try 0, 1, 3, 4. Thanks for spotting this, because it means some of my own visualisation will be wrong! On 14 September 2015 at 19:57, Cotter, Colin J <colin.cotter(a)imperial.ac.uk> wrote: > Dear friends, > I have been doing some stuff with a DG1 x CG2 field, as we use it for > temperature in the compatible FEM framework. I'm getting some funny > output when trying to interpolate. In the following code: > > https://gist.github.com/colinjcotter/96e23745d40a73303b05 > > I interpolate a function to a DG1 x CG2 field. The function represented > should > be continuous. > > Then, I remap it using a Firedrake kernel to a DG1 x DG1 field. The > function > should still be continuous. > > When I visualise the field in Paraview, I see that it has jumps in the > vertical > direction. > > What went wrong? > > cheers > --cjc > > -- > http://www.imperial.ac.uk/people/colin.cotter > > www.cambridge.org/9781107663916 > > >

1 0

results from interpolate not as expected
by Colin Cotter 14 Sep '15

14 Sep '15

Dear friends, I have been doing some stuff with a DG1 x CG2 field, as we use it for temperature in the compatible FEM framework. I'm getting some funny output when trying to interpolate. In the following code: https://gist.github.com/colinjcotter/96e23745d40a73303b05 I interpolate a function to a DG1 x CG2 field. The function represented should be continuous. Then, I remap it using a Firedrake kernel to a DG1 x DG1 field. The function should still be continuous. When I visualise the field in Paraview, I see that it has jumps in the vertical direction. What went wrong? cheers --cjc -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0

column-based block-Jacobi preconditioner in firedrake
by Eike Mueller 11 Sep '15

11 Sep '15

Dear firedrakers (and Colin), is there an easy way of building a PETSc-based block-Jacobi preconditioner for a mixed system in firedrake? So if I have the mixed system in velocity, buoyancy, pressure, can I tell PETSc that the dofs in one vertical column belong together, i.e. the full vector has a block-structure, and then use a block-Jacobi preconditioner? Or does that make no sense if I have a horizontally continuous space? Colin, in your full model, do you use the Schur-complement preconditioner and have you tried replacing it with a block-Jacobi)? Do you have a non-linear shallow water code lying around, in which we would test replacing the Schur-complement preconditioner with point-Jacobi? We’re thinking about using this as a short term alternative for the full Schur-complement preconditioner in GungHo. Thanks, Eike

2 3

Re: [firedrake] Solving this system of nonlinear equations
by David Ham 11 Sep '15

11 Sep '15

Hi Justin, I think the general approach could work. You're not solving exactly the same equation - because you have a Galerkin projection inside your reaction term, but it's not obvious that that's problematic. I'm not sure that it's the cause of your current problem, but you might try a direct solver for your U equation. This eliminates the possibility that jacobi is a bad preconditioner for this problem. Are you sure that the values of the concentrations going into the solve make sense? You might try a test case where you don't solve advection-diffusion at all but just solve the nonlinear problem for concentration based on prescribed initial conditions. You should also try applying Lawrence's collected wisdom, which is to be found at http://firedrakeproject.org/solving-interface.html#debugging-convergence-fa… Let us know how you get on. Regards, David On Fri, 11 Sep 2015 at 00:57 Justin Chang <jychang48(a)gmail.com> wrote: > Hi everyone, > > Attached is a tentative working code I have for an > advective-diffusive-reactive system using operator splitting. I tried to > express the system of equations as a mixed finite element form (not sure if > this is the right approach). > > When I run the attached code, i get an error saying "Nonlinear solver > failed to converge after 0 nonlinear iterations. Reason: > DIVERGED_LINE_SEARCH". > Can you guys have a look and see if this is even a logical approach? The > parts in question begin at line 133, and the solver is at line 217 (if you > comment out line 217, the advection-diffusion part runs well). > > If this doesn't work, I will explore the pyop2 option. > > Thanks, > Justin > > > On Fri, Sep 4, 2015 at 4:35 AM, David Ham <David.Ham(a)imperial.ac.uk> > wrote: > >> I'm not totally sure that's obvious. I agree that it's not the >> theoretically optimal solution, however it is well engineered which may >> make it a convenient solution, and if it's fast enough, then that's just >> fine. >> >> On Fri, 4 Sep 2015 at 11:22 Cotter, Colin J <colin.cotter(a)imperial.ac.uk> >> wrote: >> >>> PETSc shouldn't be used for this, as the problem is solved independently >>> at each node. >>> >>> --cjc >>> >>> On 4 September 2015 at 10:40, Ham, David A <david.ham(a)imperial.ac.uk> >>> wrote: >>> >>>> OK. I think I understand. In essence, assembling PETSc matrices is a >>>> core PyOP2 competence, so it *can* do this, although it'll be a little >>>> clunky because we don't have the full glory of Firedrake on top. >>>> >>>> I think it would help if you try to read some of the PyOP2 >>>> documentation for a bit of background: >>>> >>>> http://op2.github.io/PyOP2/concepts.html >>>> http://op2.github.io/PyOP2/kernels.html >>>> http://op2.github.io/PyOP2/linear_algebra.html >>>> http://op2.github.io/PyOP2/mixed.html >>>> >>>> The last one is important in this case because you're actually solving >>>> a mixed system over the mixed function space (A, B, C). >>>> >>>> Inconveniently, that means we need 3 kernels to assemble the residual >>>> and 7 kernels (one for each nonzero entry) to assemble the Jacobian. >>>> >>>> Once you have the code for assembling, I think it should be >>>> straightforward to provide PETSc4py with the right callbacks so that its >>>> variational inequality solver works. >>>> >>>> (Another possibility would be to bypass PyOP2 and write your own code >>>> accessing PETSc4py directly. That might actually be easier but assembly >>>> will be slower - which you may or may not care about). >>>> >>>> I don't have time right this instant to start writing out the Kernels >>>> (I have to teach on Monday and I'm not ready yet.) I'll try to get back to >>>> this later, or maybe someone else on the list will dive in and save me ;) >>>> >>>> Cheers, >>>> >>>> David >>>> >>> On Fri, 4 Sep 2015 at 09:37 Justin Chang <jychang48(a)gmail.com> wrote: >>>> >>> Hi everyone, >>>>> >>>>> Thanks for the input. So I was leaning more towards employing the >>>>> newton method for solving this. And yes, all the vectors live in the same >>>>> function space. If the system of equations were solved at each node i, the >>>>> residual algebraic vector of size 3x1 would be: >>>>> >>>>> r = [c_A(i)+c_C(i)-psi_A(i); >>>>> c_B(i)+c_C(i)-psi_B(i); >>>>> k1*c_A(i)*c_B(i)-c_C(i)] >>>>> >>>>> and the Jacobian 3x3 matrix would be: >>>>> >>>>> J = [1, 0, 1; >>>>> 0, 1, 1; >>>>> k1*c_B(i), k1*c_A(i), -1] >>>>> >>>>> Ideally I would like to have these expressed as a PETSc matrix and >>>>> vector because I would like to employ PETSc’s variational inequality solver >>>>> to enforce lower/upper bounds for c_A/B/C. >>>>> >>>>> As for the pyop2 option: >>>>> >>>>> void solve_system(double *c_A, double *c_B, double *c_C, double *k1, >>>>>>> double *psi_A, double *psi_B) >>>>>>> { >>>>>>> >>>>>>> // In here, you can dereference these variables to get the >>>>>>> components >>>>>>> // of each variable at the node, e.g. >>>>>>> >>>>>>> // Set first component of c_A to sum of first two components of >>>>>>> psi_A. >>>>>>> c_A[0] = psi_A[0] + psi_A[1]; >>>>>>> } >>>>>>> >>>>>>> I guess you can either write the newton (or fixed point) iteration >>>>>>> out by hand in the kernel (it shouldn't be too difficult), or else make >>>>>>> library calls (getting the correct link lines is possible I think). >>>>>>> >>>>>>> Now in your code you do: >>>>>>> >>>>>>> solve_kernel = op2.Kernel(kernel_code, "solve_system") >>>>>>> >>>>>>> op2.par_loop(solve_kernel, >>>>>>> c_A.dof_dset.node_set, >>>>>>> c_A.dat(op2.WRITE), >>>>>>> c_B.dat(op2.WRITE), >>>>>>> c_C.dat(op2.WRITE), >>>>>>> k1.dat(op2.READ), >>>>>>> psi_A.dat(op2.READ), >>>>>>> psi_B.dat(op2.READ)) >>>>>>> >>>>>>> >>>>>>> Does this make sense? >>>>>>> >>>>>> >>>>> I am guessing void solve_system(…) is a separate C file called >>>>> solve_system.c but how exactly do I get the components of each variable at >>>>> the node? So here’s how I did this with scipy: >>>>> >>>>> ================================ >>>>> >>>>> # Advection-diffusion firedrake code >>>>> # Compute psi_A, psi_B, and k1 ... >>>>> … >>>>> >>>>> from scipy.optimize import solve >>>>> import math >>>>> >>>>> # Solve system of nonlinear reactions >>>>> def equations(p,psi_A,psi_B,k_1): >>>>> c_A,c_B,c_C = p >>>>> return (c_A+c_C-psi_A, c_B+c_C-psi_B, c_C-k_1*c_A*c_B) >>>>> >>>>> # define vectors for primary species >>>>> c_A_vec = np.zeros(len(psi_A.vector())) >>>>> c_B_vec = np.zeros(len(psi_A.vector())) >>>>> >>>>> # define vector for secondary species >>>>> c_C_vec = np.zeros(len(psi_A.vector())) >>>>> >>>>> # Iterate through each node >>>>> for i in range(len(psi_A.vector())): >>>>> c_A_i,c_B_i,c_C_i = >>>>> fsolve(equations,(0.0,0.0,0.0),args=(psi_A.vector()[i],psi_B.vector()[i],k1) >>>>> c_A_vec[i] = c_A_i >>>>> c_B_vec[i] = c_B_i >>>>> c_C_vec[i] = c_C_i >>>>> >>>>> # Output >>>>> c_A = Function(Q) >>>>> c_B = Function(Q) >>>>> c_C = Function(Q) >>>>> c_A.vector()[:] = c_A_vec >>>>> c_B.vector()[:] = c_B_vec >>>>> c_C.vector()[:] = c_C_vec >>>>> File(“Solution_c_A.pvd”) << c_A >>>>> File(“Solution_c_B.pvd”) << c_B >>>>> File(“Solution_c_C.pvd”) << c_C >>>>> >>>>> =================================== >>>>> >>>>> The PyOP2 equivalent of the above would be a good starting point I >>>>> suppose. >>>>> >>>>> Thanks all, >>>>> Justin >>>>> >>>>> >>> -- >>> 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 >> >> >

1 0

Re: [firedrake] The wonderful world of firedrake-install
by David Ham 07 Sep '15

07 Sep '15

Hi Justin, The prospective fix is in install-cython-workaround . Can you try that and tell me if it works for you? It will be merged to master as soon as Travis gives the green light. Cheers, David On Fri, 4 Sep 2015 at 09:42 Ham, David A <david.ham(a)imperial.ac.uk> wrote: > Hi Justin, > > Ah yes. I saw something like this when attempting to set up travis. The > problem appears to be that the dependency analysis in pip is broken so it > doesn't realise it has to build and install Cython before it builds h5py. I > think the workaround is to force Cython installation in the install script > before looking at the requirements file. I will push a branch presently and > merge it as soon as it passes tests. > > Cheers, > > David > > On Fri, 4 Sep 2015 at 08:17 Justin Chang <jychang48(a)gmail.com> wrote: > >> Hi David >> >> Thanks for the updated script. So I am using the system provided >> openmpi/1.8.3 module but built my own Python 2.7.10 and pip from source. I >> am getting an error as shown in the attached firedrake-install.log >> >> I also got a similar issue when I built my own openmpi library, although >> in that case I didn’t have the cython related errors and only had the ones >> related to hdf5/h5py >> >> Any idea what this is? >> >> Thanks, >> Justin >> >> >> On Sep 3, 2015, at 8:54 AM, David Ham <David.Ham(a)imperial.ac.uk> wrote: >> >> The branch is merged now so I believe this is now fixed. >> >> David >> >> On Thu, 3 Sep 2015 at 12:09 David Ham <David.Ham(a)imperial.ac.uk> wrote: >> >> It turns out that this is for the usual reason: I am an idiot. A fix for >> this is in the further-install-fixes branch and I will merge the branch as >> soon as tests pass. >> >> (What's actually going on is that I am abusing CalledProcessError when I >> should really raise a special exception for this case. I have now created a >> new exception and raised that instead). >> >> Cheers, >> >> David >> >> On Thu, 3 Sep 2015 at 01:56 Justin Chang <jychang48(a)gmail.com> wrote: >> >> Thanks David, >> >> I actually forgot about the --no_package_manager option. But when I tried >> installing with it, i get this error: >> >> $ python firedrake-install --no_package_manager >> Traceback (most recent call last): >> File "firedrake-install", line 265, in <module> >> raise subprocess.CalledProcessError >> TypeError: __init__() takes at least 3 arguments (1 given) >> >> Even if I add the --disable_ssh and --developer tags I get the above >> error. Am I suppose to do something else with this? >> >> Thanks, >> Justin >> >> On Tue, Sep 1, 2015 at 3:57 AM, David Ham <David.Ham(a)imperial.ac.uk> >> wrote: >> >> Thanks, >> >> On Sun, 30 Aug 2015 at 08:15 Justin Chang <jychang48(a)gmail.com> wrote: >> >> I agree that hype and a sparse direct solver like mumps is necessary. >> The latter needs --download-metis --download-parmetis >> --download-scalapack. >> >> >> Added and currently testing. >> >> >> More comments: >> >> 4) If mumps is to be included, the install script has to ensure that >> one has cmake >=2.5 when configuring metis, otherwise PETSc configure >> will return an error. Perhaps include "brew install cmake" and "sudo >> apt-get install cmake" in the script? >> >> >> Also added. Thanks for catching this. >> >> >> 5) Should --download-exodusii (and it's dependencies --download-netcdf >> --download-hdf5) also be included in the PETSC_CONFIGURE_OPTIONS? I >> imagine that if one wants to solve problems on large-scale >> unstructured grids, .exo files would be more feasible to use over .msh >> files. Or is exodusii already accounted for in firedrake/PyOP2/etc? >> >> >> And this. >> >> >> 6) If I run the script on a pristine MacOSX, I get an error saying >> "OSError: Unable to find SWIG installation. Please install SWIG >> version 2.0.0 or higher." I am guessing "brew install swig" was >> somehow missed in the firedrake-install script? >> >> Ah yes. Fecking swig. I've added the dependency to the new branch. >> Hopefully the swig dependency is going away very soon. >> >> >> >> 7) Perhaps not as important, but do y'all think it's possible to make >> this script more non-Ubuntu HPC system friendly? For isntance, either >> a) require the user to install his or her own local OpenMPI and Python >> libraries from source, b) let script download the tarballs/git >> repositories remotely and install it for you (kind of like how PETSc >> handles external packages through --download-<package>), or c) enable >> the user to point to the system provided OpenMPI and Python libraries. >> Same would need to be done for CMake, SWIG, and PCRE as these packages >> cannot simply be obtained from pip. >> >> >> I think (c) is the case right now. If you run the script on a non-Ubuntu >> system or pass the --no_package_manager option then it will just assume >> that you have installed the dependencies somewhere visible and will get on >> with pip. >> >> I think (b) is [a] hard and [b] in the case of MPI probably a bad idea. >> Supercomputers are such a diverse bunch of hacked up distributions that I >> think that installing the compiled dependencies automatically in a way >> which is portable would require a lot of install script work. In the case >> of MPI, you really need to use the MPI build which talks to the fast >> interconnect on your supercomputer, so automatically building vanilla >> OpenMPI is a Bad Thing. >> >> Cheers, >> >> David >> >> >> >> >> Thanks, >> Justin >> >> On Fri, Aug 28, 2015 at 7:20 AM, Lawrence Mitchell >> <lawrence.mitchell(a)imperial.ac.uk> wrote: >> > >> >> On 28 Aug 2015, at 14:01, David Ham <David.Ham(a)imperial.ac.uk> wrote: >> >> >> >> Hi All, >> >> >> >> I have a branch going through testing now which will handle >> PETSC_CONFIGURE_OPTIONS in a smarter way. Basically we'll just make sure >> that the things which are required are there, and we'll also honour >> anything the user has already set. I've also updated docs to say that >> that's what happens. >> >> >> >> I am also open to the suggestion that we could add more configuration >> options to the default set. Would anyone like to suggest what they should >> be? >> > >> > I think --download-hypre=1 at least. Maybe also a sparse direct solver >> (mumps?) which I think needs --download-metis --download-parmetis >> --download-mumps (maybe some others?) >> > >> > 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 >> >> >> _______________________________________________ >> 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

Re: [firedrake] Solving this system of nonlinear equations
by Colin Cotter 04 Sep '15

04 Sep '15

PETSc shouldn't be used for this, as the problem is solved independently at each node. --cjc On 4 September 2015 at 10:40, Ham, David A <david.ham(a)imperial.ac.uk> wrote: > OK. I think I understand. In essence, assembling PETSc matrices is a core > PyOP2 competence, so it *can* do this, although it'll be a little clunky > because we don't have the full glory of Firedrake on top. > > I think it would help if you try to read some of the PyOP2 documentation > for a bit of background: > > http://op2.github.io/PyOP2/concepts.html > http://op2.github.io/PyOP2/kernels.html > http://op2.github.io/PyOP2/linear_algebra.html > http://op2.github.io/PyOP2/mixed.html > > The last one is important in this case because you're actually solving a > mixed system over the mixed function space (A, B, C). > > Inconveniently, that means we need 3 kernels to assemble the residual and > 7 kernels (one for each nonzero entry) to assemble the Jacobian. > > Once you have the code for assembling, I think it should be > straightforward to provide PETSc4py with the right callbacks so that its > variational inequality solver works. > > (Another possibility would be to bypass PyOP2 and write your own code > accessing PETSc4py directly. That might actually be easier but assembly > will be slower - which you may or may not care about). > > I don't have time right this instant to start writing out the Kernels (I > have to teach on Monday and I'm not ready yet.) I'll try to get back to > this later, or maybe someone else on the list will dive in and save me ;) > > Cheers, > > David > > On Fri, 4 Sep 2015 at 09:37 Justin Chang <jychang48(a)gmail.com> wrote: > >> Hi everyone, >> >> Thanks for the input. So I was leaning more towards employing the newton >> method for solving this. And yes, all the vectors live in the same function >> space. If the system of equations were solved at each node i, the residual >> algebraic vector of size 3x1 would be: >> >> r = [c_A(i)+c_C(i)-psi_A(i); >> c_B(i)+c_C(i)-psi_B(i); >> k1*c_A(i)*c_B(i)-c_C(i)] >> >> and the Jacobian 3x3 matrix would be: >> >> J = [1, 0, 1; >> 0, 1, 1; >> k1*c_B(i), k1*c_A(i), -1] >> >> Ideally I would like to have these expressed as a PETSc matrix and vector >> because I would like to employ PETSc’s variational inequality solver to >> enforce lower/upper bounds for c_A/B/C. >> >> As for the pyop2 option: >> >> void solve_system(double *c_A, double *c_B, double *c_C, double *k1, >>>> double *psi_A, double *psi_B) >>>> { >>>> >>>> // In here, you can dereference these variables to get the components >>>> // of each variable at the node, e.g. >>>> >>>> // Set first component of c_A to sum of first two components of psi_A. >>>> c_A[0] = psi_A[0] + psi_A[1]; >>>> } >>>> >>>> I guess you can either write the newton (or fixed point) iteration out >>>> by hand in the kernel (it shouldn't be too difficult), or else make library >>>> calls (getting the correct link lines is possible I think). >>>> >>>> Now in your code you do: >>>> >>>> solve_kernel = op2.Kernel(kernel_code, "solve_system") >>>> >>>> op2.par_loop(solve_kernel, >>>> c_A.dof_dset.node_set, >>>> c_A.dat(op2.WRITE), >>>> c_B.dat(op2.WRITE), >>>> c_C.dat(op2.WRITE), >>>> k1.dat(op2.READ), >>>> psi_A.dat(op2.READ), >>>> psi_B.dat(op2.READ)) >>>> >>>> >>>> Does this make sense? >>>> >>> >> I am guessing void solve_system(…) is a separate C file called >> solve_system.c but how exactly do I get the components of each variable at >> the node? So here’s how I did this with scipy: >> >> ================================ >> >> # Advection-diffusion firedrake code >> # Compute psi_A, psi_B, and k1 ... >> … >> >> from scipy.optimize import solve >> import math >> >> # Solve system of nonlinear reactions >> def equations(p,psi_A,psi_B,k_1): >> c_A,c_B,c_C = p >> return (c_A+c_C-psi_A, c_B+c_C-psi_B, c_C-k_1*c_A*c_B) >> >> # define vectors for primary species >> c_A_vec = np.zeros(len(psi_A.vector())) >> c_B_vec = np.zeros(len(psi_A.vector())) >> >> # define vector for secondary species >> c_C_vec = np.zeros(len(psi_A.vector())) >> >> # Iterate through each node >> for i in range(len(psi_A.vector())): >> c_A_i,c_B_i,c_C_i = >> fsolve(equations,(0.0,0.0,0.0),args=(psi_A.vector()[i],psi_B.vector()[i],k1) >> c_A_vec[i] = c_A_i >> c_B_vec[i] = c_B_i >> c_C_vec[i] = c_C_i >> >> # Output >> c_A = Function(Q) >> c_B = Function(Q) >> c_C = Function(Q) >> c_A.vector()[:] = c_A_vec >> c_B.vector()[:] = c_B_vec >> c_C.vector()[:] = c_C_vec >> File(“Solution_c_A.pvd”) << c_A >> File(“Solution_c_B.pvd”) << c_B >> File(“Solution_c_C.pvd”) << c_C >> >> =================================== >> >> The PyOP2 equivalent of the above would be a good starting point I >> suppose. >> >> Thanks all, >> Justin >> >> -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916

1 0