August 2016 - firedrake - imperial.ac.uk

reordering meshes
by Justin Chang 12 Aug '16

12 Aug '16

Hi all, >From what I see in the documentation, reordering of meshes is turned on by default. But when I look at the mesh generation functions like RectangleMesh(...), I see that one of the default options is "reorder=None". Am i actually reordering the mesh if I ran something like: mesh = UnitSquareMesh(10,10) Thanks, Justin

3 3

Re: [firedrake] Mixed system
by David Ham 10 Aug '16

10 Aug '16

The short version is it depends on when I find some time to work on it. It is not simply the case that we have new features lying around waiting to be merged. The current proximate issue is to make the real spaces work properly with MPI communicators. I can't guarantee that will be the last issue. On Tue, 9 Aug 2016 at 17:37 Onno Bokhove <O.Bokhove(a)leeds.ac.uk> wrote: > It is possible to merge the branch for the mixed system relatively soon? > > Anna's results are relevant for the grant application Colin and I aim to > submit in Sept. on further wave-ship modelling. > > Sorry for asking! > ------------------------------ > *From:* firedrake-bounces(a)imperial.ac.uk <firedrake-bounces(a)imperial.ac.uk> > on behalf of David Ham <David.Ham(a)imperial.ac.uk> > *Sent:* Friday, August 5, 2016 11:37:23 AM > *To:* firedrake > *Subject:* Re: [firedrake] Mixed system > > > > On Fri, 5 Aug 2016 at 11:31 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> > wrote: > >> Hi Kramer, >> >> Thank you for the reply. >> >> The only problem I see with the solution you suggest, is that both I and >> lambda will be trial functions, so both terms should be on the bilinear >> part of the variational form. >> > > Not quite. I and lambda together form a Function in a mixed finite element > space. You solve for both of them together, just like pressure and velocity > in a fluid simulation. > > >> >> I could, of course, define a residual and solve everything using >> nonlinear solvers, but I don't know how efficient that would be. >> > > I think you are confused. None of the relevant issues change in the > nonlinear case. After all, the nonlinear solver will differentiate the > residual to create a bilinear form anyway. > > Incidentally, there is no particular reason to believe that writing the > problem in nonlinear form would be inefficient. > > Regards, > > David > > >> >> Best, >> >> Anna. >> >> Dr Anna Kalogirou >> Research Fellow >> School of Mathematics >> University of Leeds >> >> http://www1.maths.leeds.ac.uk/~matak/ >> >> On 04/08/16 17:29, Stephan Kramer wrote: >> > On 27/07/16 11:01, Anna Kalogirou wrote: >> >> Dear all, >> >> >> >> I am trying to solve the system found in the attached pdf as a mixed >> >> system. I want to solve the same >> >> problem as the one found here >> >> < >> https://bitbucket.org/annakalog/buoy2d/src/14334d3c20b9f10ed7c1246cde9e3cb6… >> >, >> >> >> >> but with the use of Schur complements and not linear algebra. >> >> >> >> In a previous discussion about this problem, Lawrence mentioned that >> the >> >> test function for integral(lambda*Theta(x-Lp)dx) needs to be considered >> >> as coming from the real space of constant functions. Could you please >> >> elaborate on this? >> > >> > My guess at what Lawrence means is (it's how I would probably approach >> > it), is that you add an additional equation that says >> > >> > (1) I = integral(lambda*Theta(x-Lp)dx) >> > >> > and then substitute I in your third equation. Equation (1) can >> > be turned into a finite element weak formulation by considering the >> > function space of functions that are constant over the entire domain >> > This function space is just 1-dimensional and therefore equivalent to >> > the space of reals R and hence we denote this function space of >> > constants simply as R. Then we can consider I to be a trial function >> > in R, and with a test function v4 in R, we can write: >> > >> > v4*lambda*Theta(x-Lp)*dx == v4*I*dx/area(domain) >> > >> > This is a proper finite element equation that you could solve in >> > conjuction with the other equations. I believe however that this >> > function space R is not currently implemented in Firedrake. It is >> > available in fenics, and is on the wishlist. >> > >> > Cheers >> > Stephan >> > >> > _______________________________________________ >> > 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 >> >

3 2

Re: [firedrake] Mixed system
by David Ham 09 Aug '16

09 Aug '16

On Fri, 5 Aug 2016 at 11:31 Anna Kalogirou <a.kalogirou(a)leeds.ac.uk> wrote: > Hi Kramer, > > Thank you for the reply. > > The only problem I see with the solution you suggest, is that both I and > lambda will be trial functions, so both terms should be on the bilinear > part of the variational form. > Not quite. I and lambda together form a Function in a mixed finite element space. You solve for both of them together, just like pressure and velocity in a fluid simulation. > > I could, of course, define a residual and solve everything using > nonlinear solvers, but I don't know how efficient that would be. > I think you are confused. None of the relevant issues change in the nonlinear case. After all, the nonlinear solver will differentiate the residual to create a bilinear form anyway. Incidentally, there is no particular reason to believe that writing the problem in nonlinear form would be inefficient. Regards, David > > Best, > > Anna. > > Dr Anna Kalogirou > Research Fellow > School of Mathematics > University of Leeds > > http://www1.maths.leeds.ac.uk/~matak/ > > On 04/08/16 17:29, Stephan Kramer wrote: > > On 27/07/16 11:01, Anna Kalogirou wrote: > >> Dear all, > >> > >> I am trying to solve the system found in the attached pdf as a mixed > >> system. I want to solve the same > >> problem as the one found here > >> < > https://bitbucket.org/annakalog/buoy2d/src/14334d3c20b9f10ed7c1246cde9e3cb6… > >, > >> > >> but with the use of Schur complements and not linear algebra. > >> > >> In a previous discussion about this problem, Lawrence mentioned that the > >> test function for integral(lambda*Theta(x-Lp)dx) needs to be considered > >> as coming from the real space of constant functions. Could you please > >> elaborate on this? > > > > My guess at what Lawrence means is (it's how I would probably approach > > it), is that you add an additional equation that says > > > > (1) I = integral(lambda*Theta(x-Lp)dx) > > > > and then substitute I in your third equation. Equation (1) can > > be turned into a finite element weak formulation by considering the > > function space of functions that are constant over the entire domain > > This function space is just 1-dimensional and therefore equivalent to > > the space of reals R and hence we denote this function space of > > constants simply as R. Then we can consider I to be a trial function > > in R, and with a test function v4 in R, we can write: > > > > v4*lambda*Theta(x-Lp)*dx == v4*I*dx/area(domain) > > > > This is a proper finite element equation that you could solve in > > conjuction with the other equations. I believe however that this > > function space R is not currently implemented in Firedrake. It is > > available in fenics, and is on the wishlist. > > > > Cheers > > Stephan > > > > _______________________________________________ > > 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 >

3 2

Why is firedrake so fast?
by Justin Chang 09 Aug '16

09 Aug '16

Hi all, So I have some personal computers which will give me roughly the same serial performance (in terms of wall-clock time) as a single compute node (via interactive mode e.g., salloc -N 1 -n 16) at one of LANL's machines (Intel Xeon E5 2670). However, once I submit my firedrake program as a job script (via MOAB) to said LANL machine, the performance becomes 10x faster! Certainly not complaining, but I find this unusual. I never encountered such a performance improvement if I were to submit a normal PETSc program. Does this have something to do with the internal framework of firedrake? I can provide more information as needed Thanks, Justin

2 2

Re: [firedrake] Installing firedrake without internet connection
by David Ham 09 Aug '16

09 Aug '16

Hi Floriane, I'm afraid that what firedrake-install does is downloads and builds firedrake and all its components so it's not going to work without internet access. Also note that firedrake-update works by downloading updates to firedrake so even if you manage to install once with no network, updating is not going to work. Now, what you in principle can do is manually set up all the dependencies of firedrake on the cluster (assuming you can find a way of copying them to the cluster) and then build everything in a virtualenv yourself, but this is not going to be a lot of fun. You should also be aware that cluster environments are somewhat different from the usual PC environment. You in particular need to make sure you pick up the officially blessed MPI implementation in order to access the high speed network correctly. It may also be necessary to load environment modules or take other steps in order to access the right compilers, Python or even git. Regards, David On Mon, 8 Aug 2016 at 16:06 Floriane Gidel [RPG] <mmfg(a)leeds.ac.uk> wrote: > Dear all, > > > Is it possible to install Firedrake without internet connection ? I > managed to get the firedrake-install file, and I get the following error > when running it: > > > [fgidel@marclus3login2 ~]$ python firedrake-install --no-package-manager > apt-get not found or disabled. Proceeding on the rash assumption that your > compiled dependencies are in place. > If this is not the case, please install the following and try again: > * A C and C++ compiler (for example gcc/g++ or clang), GNU make > * A Fortran compiler (for PETSc) > * MPI > * Blas and Lapack > * Git, Mercurial > * Python version 2.7 > * pip and the Python headers > * autoconf, automake, libtool > * CMake > * zlib > Creating firedrake virtualenv. > /home/fgidel/firedrake/bin/pip install -U pip > > Requirement already up-to-date: pip in > ./firedrake/lib/python2.7/site-packages > Cloning firedrake > > Failed to clone firedrake using ssh, falling back to https. > Traceback (most recent call last): > File "firedrake-install", line 857, in <module> > git_clone("git+https://github.com/firedrakeproject/firedrake.git") > File "firedrake-install", line 285, in git_clone > check_call(["git", "clone", "-q", "-b", branch, git_url(plain_url, > "https")]) > File "firedrake-install", line 205, in check_call > subprocess.check_call(arguments, env=env) > File "/usr/lib64/python2.7/subprocess.py", line 537, in check_call > retcode = call(*popenargs, **kwargs) > File "/usr/lib64/python2.7/subprocess.py", line 524, in call > return Popen(*popenargs, **kwargs).wait() > File "/usr/lib64/python2.7/subprocess.py", line 711, in __init__ > errread, errwrite) > File "/usr/lib64/python2.7/subprocess.py", line 1327, in _execute_child > raise child_exception > OSError: [Errno 2] No such file or directory > > > I know that the cluster on which I am trying to install it has no access > to internet. Is it the reason for the error ? How could I manage to install > firedrake despite that ? > > > Thanks and best regards, > > Floriane >

2 1

Installing firedrake without internet connection
by Floriane Gidel [RPG] 08 Aug '16

08 Aug '16

Dear all, Is it possible to install Firedrake without internet connection ? I managed to get the firedrake-install file, and I get the following error when running it: [fgidel@marclus3login2 ~]$ python firedrake-install --no-package-manager apt-get not found or disabled. Proceeding on the rash assumption that your compiled dependencies are in place. If this is not the case, please install the following and try again: * A C and C++ compiler (for example gcc/g++ or clang), GNU make * A Fortran compiler (for PETSc) * MPI * Blas and Lapack * Git, Mercurial * Python version 2.7 * pip and the Python headers * autoconf, automake, libtool * CMake * zlib Creating firedrake virtualenv. /home/fgidel/firedrake/bin/pip install -U pip Requirement already up-to-date: pip in ./firedrake/lib/python2.7/site-packages Cloning firedrake Failed to clone firedrake using ssh, falling back to https. Traceback (most recent call last): File "firedrake-install", line 857, in <module> git_clone("git+https://github.com/firedrakeproject/firedrake.git") File "firedrake-install", line 285, in git_clone check_call(["git", "clone", "-q", "-b", branch, git_url(plain_url, "https")]) File "firedrake-install", line 205, in check_call subprocess.check_call(arguments, env=env) File "/usr/lib64/python2.7/subprocess.py", line 537, in check_call retcode = call(*popenargs, **kwargs) File "/usr/lib64/python2.7/subprocess.py", line 524, in call return Popen(*popenargs, **kwargs).wait() File "/usr/lib64/python2.7/subprocess.py", line 711, in __init__ errread, errwrite) File "/usr/lib64/python2.7/subprocess.py", line 1327, in _execute_child raise child_exception OSError: [Errno 2] No such file or directory I know that the cluster on which I am trying to install it has no access to internet. Is it the reason for the error ? How could I manage to install firedrake despite that ? Thanks and best regards, Floriane

1 0

firedrake-install
by Onno Bokhove 08 Aug '16

08 Aug '16

We are installing from scratch and get the following; access to brew is in order: ... /usr/local/Cellar/libtool/2.4.6_1/share/man/ (2 files) Creating firedrake virtualenv in 'firedrake'. ERROR: The executable firedrake/bin/python2.7 is not functioning ERROR: It thinks sys.prefix is u'/usr/local/Cellar/python/2.7.10_2/Frameworks/Python.framework/Versions/2.7' (should be u'/Users/amtob/werk/firedrake') ERROR: virtualenv is not compatible with this system or executable

2 1

Re: [firedrake] 3D SUPG on extruded mesh not working
by Andrew McRae 05 Aug '16

05 Aug '16

If you're planning to use the tet code (UnitCubeMesh) then you should do the same. I guess you snuck just under the 1000-point error check on that, due to small differences in degree estimation. Without a check like that, users would be none-the-wiser and just complain "Firedrake is slow" :) On 5 August 2016 at 23:09, Justin Chang <jychang48(a)gmail.com> wrote: > That did the trick, thanks Andrew > > On Fri, Aug 5, 2016 at 5:01 PM, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> > wrote: > >> If you want your a quadrature rule which is exact for *form* polynomial >> degree 5 in the horizontal and 3 in the vertical, use dx(degree=(5, 3)). >> >> UFL symbolically analyses the form to obtain the polynomial degree, and >> by default uses quadrature which is exact for this. When non-polynomial >> operations are present, it uses a heuristic, and when lots are present, it >> tends to be far too conservative. >> >> E.g. deg(a/b) is estimated as deg(a) + deg(b), and deg(sqrt(f)) is >> estimated as 2 + deg(f). These are reasonable on their own, but now your >> classic unit vector v/sqrt(dot(v, v)) has degree >> deg(v) + deg(sqrt(dot(v, v))) >> = deg(v) + 2 + deg(dot(v, v)) >> = deg(v) + 2 + deg(v) + deg(v) >> If deg(v) = 1, this is 5(!) >> >> In your case, I guess it's gone for a quadrature rule with 12 points in >> each direction, so it's estimated the total polynomial degree as 22 or >> 23... probably far more than you want. >> >> On 5 August 2016 at 22:52, Justin Chang <jychang48(a)gmail.com> wrote: >> >>> Andrew, where would I do that? >>> >>> >>> On Friday, August 5, 2016, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: >>> >>>> Set a manual quadrature degree, you probably don't want 1728 points per >>>> cell. >>>> >>>> On 5 August 2016 at 22:47, Justin Chang <jychang48(a)gmail.com> wrote: >>>> >>>>> Hi all, >>>>> >>>>> I have the following 3D advection diffusion problem with SUPG >>>>> formulation: >>>>> >>>>> # Mesh >>>>> meshbase = UnitSquareMesh(seed,seed,quadrilateral=True) >>>>> mesh = ExtrudedMesh(meshbase,seed) >>>>> V = FunctionSpace(mesh, 'CG', 1) >>>>> Q = VectorFunctionSpace(mesh, 'CG', 1) >>>>> u = TrialFunction(V) >>>>> v = TestFunction(V) >>>>> sol = Function(V, name="solution") >>>>> >>>>> # Velocity >>>>> velocity = interpolate(Expression(( ... )),Q) >>>>> >>>>> # Dispersion tensor >>>>> aT = 0.0001 >>>>> aL = 1 >>>>> aD = 0.0001 >>>>> alpha_t = Constant(aT) >>>>> alpha_L = Constant(aL) >>>>> alpha_D = Constant(aD) >>>>> normv = sqrt(dot(velocity,velocity)) >>>>> Id = Identity(mesh.geometric_dimension()) >>>>> D = (alpha_D + alpha_t*normv)*Id+(alpha_L-alp >>>>> ha_t)*outer(velocity,velocity)/normv >>>>> >>>>> # Forcing function >>>>> f = Constant(0.0) >>>>> >>>>> # Weak form >>>>> h = CellSize(V.mesh()) >>>>> Pe = h/(2*normv)*dot(velocity,grad(v)) >>>>> a_r = Pe*(dot(velocity,grad(u)) - div(D*grad(u)))*dx >>>>> a = a_r + v*dot(velocity,grad(u))*dx + dot(grad(v),D*grad(u))*dx >>>>> L = (v + Pe)*f*dx >>>>> >>>>> ... >>>>> >>>>> But when I try to run this code i get the following error: >>>>> >>>>> Traceback (most recent call last): >>>>> File "3D_advection_diffusion_ex1.py", line 122, in <module> >>>>> execfile('driver/AD_solver.py') >>>>> File "driver/AD_solver.py", line 36, in <module> >>>>> solver_parameters=solver_parameters, options_prefix="linear_") >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/variational_solver.py", >>>>> line 261, in __init__ >>>>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/variational_solver.py", >>>>> line 127, in __init__ >>>>> ctx = solving_utils._SNESContext(problem) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", >>>>> line 107, in __init__ >>>>> for problem in problems) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", >>>>> line 107, in <genexpr> >>>>> for problem in problems) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>>>> line 67, in assemble >>>>> inverse=inverse, nest=nest) >>>>> File "<decorator-gen-295>", line 2, in _assemble >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", >>>>> line 62, in wrapper >>>>> return f(*args, **kwargs) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>>>> line 101, in _assemble >>>>> inverse=inverse) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>>>> line 288, in compile_form >>>>> number_map).kernels >>>>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>>>> line 198, in __new__ >>>>> obj = make_obj() >>>>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>>>> line 188, in make_obj >>>>> obj.__init__(*args, **kwargs) >>>>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>>>> line 210, in __init__ >>>>> tree = tsfc_compile_form(form, prefix=name, parameters=parameters) >>>>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", >>>>> line 54, in compile_form >>>>> kernels.append(compile_integral(integral_data, fd, prefix, >>>>> parameters)) >>>>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", >>>>> line 198, in compile_integral >>>>> quadrature_degree)) >>>>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/fiatinterface.py", >>>>> line 108, in create_quadrature >>>>> raise RuntimeError("Requested a quadrature rule with %d points" % >>>>> len(rule.get_points())) >>>>> RuntimeError: Requested a quadrature rule with 1728 points >>>>> >>>>> If I had used a UnitCubeMesh I had no problems. What is going on here >>>>> and how do I fix/circumvent this issue? >>>>> >>>>> Thanks, >>>>> Justin >>>>> >>>> >>>> >> >> _______________________________________________ >> firedrake mailing list >> firedrake(a)imperial.ac.uk >> https://mailman.ic.ac.uk/mailman/listinfo/firedrake >> >> >

1 0

Re: [firedrake] 3D SUPG on extruded mesh not working
by Andrew McRae 05 Aug '16

05 Aug '16

If you want your a quadrature rule which is exact for *form* polynomial degree 5 in the horizontal and 3 in the vertical, use dx(degree=(5, 3)). UFL symbolically analyses the form to obtain the polynomial degree, and by default uses quadrature which is exact for this. When non-polynomial operations are present, it uses a heuristic, and when lots are present, it tends to be far too conservative. E.g. deg(a/b) is estimated as deg(a) + deg(b), and deg(sqrt(f)) is estimated as 2 + deg(f). These are reasonable on their own, but now your classic unit vector v/sqrt(dot(v, v)) has degree deg(v) + deg(sqrt(dot(v, v))) = deg(v) + 2 + deg(dot(v, v)) = deg(v) + 2 + deg(v) + deg(v) If deg(v) = 1, this is 5(!) In your case, I guess it's gone for a quadrature rule with 12 points in each direction, so it's estimated the total polynomial degree as 22 or 23... probably far more than you want. On 5 August 2016 at 22:52, Justin Chang <jychang48(a)gmail.com> wrote: > Andrew, where would I do that? > > > On Friday, August 5, 2016, Andrew McRae <A.T.T.McRae(a)bath.ac.uk> wrote: > >> Set a manual quadrature degree, you probably don't want 1728 points per >> cell. >> >> On 5 August 2016 at 22:47, Justin Chang <jychang48(a)gmail.com> wrote: >> >>> Hi all, >>> >>> I have the following 3D advection diffusion problem with SUPG >>> formulation: >>> >>> # Mesh >>> meshbase = UnitSquareMesh(seed,seed,quadrilateral=True) >>> mesh = ExtrudedMesh(meshbase,seed) >>> V = FunctionSpace(mesh, 'CG', 1) >>> Q = VectorFunctionSpace(mesh, 'CG', 1) >>> u = TrialFunction(V) >>> v = TestFunction(V) >>> sol = Function(V, name="solution") >>> >>> # Velocity >>> velocity = interpolate(Expression(( ... )),Q) >>> >>> # Dispersion tensor >>> aT = 0.0001 >>> aL = 1 >>> aD = 0.0001 >>> alpha_t = Constant(aT) >>> alpha_L = Constant(aL) >>> alpha_D = Constant(aD) >>> normv = sqrt(dot(velocity,velocity)) >>> Id = Identity(mesh.geometric_dimension()) >>> D = (alpha_D + alpha_t*normv)*Id+(alpha_L-alp >>> ha_t)*outer(velocity,velocity)/normv >>> >>> # Forcing function >>> f = Constant(0.0) >>> >>> # Weak form >>> h = CellSize(V.mesh()) >>> Pe = h/(2*normv)*dot(velocity,grad(v)) >>> a_r = Pe*(dot(velocity,grad(u)) - div(D*grad(u)))*dx >>> a = a_r + v*dot(velocity,grad(u))*dx + dot(grad(v),D*grad(u))*dx >>> L = (v + Pe)*f*dx >>> >>> ... >>> >>> But when I try to run this code i get the following error: >>> >>> Traceback (most recent call last): >>> File "3D_advection_diffusion_ex1.py", line 122, in <module> >>> execfile('driver/AD_solver.py') >>> File "driver/AD_solver.py", line 36, in <module> >>> solver_parameters=solver_parameters, options_prefix="linear_") >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/variational_solver.py", >>> line 261, in __init__ >>> super(LinearVariationalSolver, self).__init__(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/variational_solver.py", >>> line 127, in __init__ >>> ctx = solving_utils._SNESContext(problem) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", >>> line 107, in __init__ >>> for problem in problems) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", >>> line 107, in <genexpr> >>> for problem in problems) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>> line 67, in assemble >>> inverse=inverse, nest=nest) >>> File "<decorator-gen-295>", line 2, in _assemble >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", >>> line 62, in wrapper >>> return f(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", >>> line 101, in _assemble >>> inverse=inverse) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>> line 288, in compile_form >>> number_map).kernels >>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>> line 198, in __new__ >>> obj = make_obj() >>> File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", >>> line 188, in make_obj >>> obj.__init__(*args, **kwargs) >>> File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", >>> line 210, in __init__ >>> tree = tsfc_compile_form(form, prefix=name, parameters=parameters) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line >>> 54, in compile_form >>> kernels.append(compile_integral(integral_data, fd, prefix, >>> parameters)) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line >>> 198, in compile_integral >>> quadrature_degree)) >>> File "/home/justin/Software/firedrake/src/tsfc/tsfc/fiatinterface.py", >>> line 108, in create_quadrature >>> raise RuntimeError("Requested a quadrature rule with %d points" % >>> len(rule.get_points())) >>> RuntimeError: Requested a quadrature rule with 1728 points >>> >>> If I had used a UnitCubeMesh I had no problems. What is going on here >>> and how do I fix/circumvent this issue? >>> >>> Thanks, >>> Justin >>> >> >>

2 1

Re: [firedrake] 3D SUPG on extruded mesh not working
by Andrew McRae 05 Aug '16

05 Aug '16

Set a manual quadrature degree, you probably don't want 1728 points per cell. On 5 August 2016 at 22:47, Justin Chang <jychang48(a)gmail.com> wrote: > Hi all, > > I have the following 3D advection diffusion problem with SUPG formulation: > > # Mesh > meshbase = UnitSquareMesh(seed,seed,quadrilateral=True) > mesh = ExtrudedMesh(meshbase,seed) > V = FunctionSpace(mesh, 'CG', 1) > Q = VectorFunctionSpace(mesh, 'CG', 1) > u = TrialFunction(V) > v = TestFunction(V) > sol = Function(V, name="solution") > > # Velocity > velocity = interpolate(Expression(( ... )),Q) > > # Dispersion tensor > aT = 0.0001 > aL = 1 > aD = 0.0001 > alpha_t = Constant(aT) > alpha_L = Constant(aL) > alpha_D = Constant(aD) > normv = sqrt(dot(velocity,velocity)) > Id = Identity(mesh.geometric_dimension()) > D = (alpha_D + alpha_t*normv)*Id+(alpha_L-alpha_t)*outer(velocity, > velocity)/normv > > # Forcing function > f = Constant(0.0) > > # Weak form > h = CellSize(V.mesh()) > Pe = h/(2*normv)*dot(velocity,grad(v)) > a_r = Pe*(dot(velocity,grad(u)) - div(D*grad(u)))*dx > a = a_r + v*dot(velocity,grad(u))*dx + dot(grad(v),D*grad(u))*dx > L = (v + Pe)*f*dx > > ... > > But when I try to run this code i get the following error: > > Traceback (most recent call last): > File "3D_advection_diffusion_ex1.py", line 122, in <module> > execfile('driver/AD_solver.py') > File "driver/AD_solver.py", line 36, in <module> > solver_parameters=solver_parameters, options_prefix="linear_") > File "/home/justin/Software/firedrake/src/firedrake/ > firedrake/variational_solver.py", line 261, in __init__ > super(LinearVariationalSolver, self).__init__(*args, **kwargs) > File "/home/justin/Software/firedrake/src/firedrake/ > firedrake/variational_solver.py", line 127, in __init__ > ctx = solving_utils._SNESContext(problem) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", > line 107, in __init__ > for problem in problems) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", > line 107, in <genexpr> > for problem in problems) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", > line 67, in assemble > inverse=inverse, nest=nest) > File "<decorator-gen-295>", line 2, in _assemble > File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", > line 62, in wrapper > return f(*args, **kwargs) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", > line 101, in _assemble > inverse=inverse) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", > line 288, in compile_form > number_map).kernels > File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line > 198, in __new__ > obj = make_obj() > File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line > 188, in make_obj > obj.__init__(*args, **kwargs) > File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", > line 210, in __init__ > tree = tsfc_compile_form(form, prefix=name, parameters=parameters) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line > 54, in compile_form > kernels.append(compile_integral(integral_data, fd, prefix, > parameters)) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line > 198, in compile_integral > quadrature_degree)) > File "/home/justin/Software/firedrake/src/tsfc/tsfc/fiatinterface.py", > line 108, in create_quadrature > raise RuntimeError("Requested a quadrature rule with %d points" % > len(rule.get_points())) > RuntimeError: Requested a quadrature rule with 1728 points > > If I had used a UnitCubeMesh I had no problems. What is going on here and > how do I fix/circumvent this issue? > > Thanks, > Justin >

2 1