"Adding expressions with non-matching form arguments {0} vs {1}." It looks like you declared one of the forms incorrectly, e.g. used a trial function instead of a test function, or tried to add a form with a test function to a form with no test function, etc. On 12 March 2017 at 12:49, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
I am trying to build on the QG SLEPc demo that I put together a while ago to solve for the linear instability of a 1D Shallow Water Jet. I've done this in the past using spectral methods and have had much success, but I often get spurious modes, which requires very high resolution to remove. I have been told by a reliable source that the FEM should be less prone to these spurious modes, which would be great.
I have a notebook and a python script at the following links:
https://github.com/francispoulin/firedrakeQG/blob/master/Linear% 20Stability%20Analysis%20SW%20Jet%20(Firedrake).ipynb
https://github.com/francispoulin/firedrakeQG/blob/master/LSA_SW_Jet_ Firedrake.py
The error comes from the assemble part. I have done this previously with QG and SW in 2D and managed to get it working. This is only a 1D problem, so a little bit different, but should be easier. Unfortunately, when I do try it I get the following error, see below. I have tried different combinations but am not having any success.
One issue is that there are three variables: u, v, eta (height). I define the solution space to be V*V*V where V is a CG space. This is the major difference compared to what I tried before.
I am pretty sure that I am getting the syntax wrong here. Can anyone offer any advice?
Cheers, Francis
(firedrake) fpoulin@fpoulin-Gazelle:~/Research/2015-Ben_Storer/Firedrake$ python LSA_SW_Jet_Firedrake.py ('Ro = ', 0.25) ('Bu = ', 0.0004) ('profile = ', 'bickley') Traceback (most recent call last): File "LSA_SW_Jet_Firedrake.py", line 103, in <module> petsc_a = assemble(a, mat_type='aij', bcs=bc).M.handle File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 102, in assemble allocate_only=allocate_only) File "<decorator-gen-279>", line 2, in _assemble File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/driver.py", line 44, in compile_form fd = ufl_utils.compute_form_data(form) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/home/fpoulin/software/firedrake/src/ufl/ufl/ algorithms/compute_form_data.py", line 387, in compute_form_data check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 152, in check_form_arity check_integrand_arity(itg.integrand(), arguments) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 145, in check_integrand_arity args = map_expr_dag(rules, expr, compress=False) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 37, in map_expr_dag result, = map_expr_dags(function, [expression], compress=compress) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 86, in map_expr_dags r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands]) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 42, in sum raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(a, b)) ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments (Argument(WithGeometry( FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None), Argument(WithGeometry( FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 1, None)) vs (Argument(WithGeometry( FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None),).
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 <(519)%20888-4567>
You are computing Eigenvalues from a problem which involves Dirichlet boundary conditions but don't shift them out of your spectrum. So your if your Eigenvalues are around 0 real part, this leads to spurious Eigenmodes which arise from Dirichlet values on the boundary. On Sun, Mar 12, 2017 at 2:22 PM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
"Adding expressions with non-matching form arguments {0} vs {1}."
It looks like you declared one of the forms incorrectly, e.g. used a trial function instead of a test function, or tried to add a form with a test function to a form with no test function, etc.
On 12 March 2017 at 12:49, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
I am trying to build on the QG SLEPc demo that I put together a while ago to solve for the linear instability of a 1D Shallow Water Jet. I've done this in the past using spectral methods and have had much success, but I often get spurious modes, which requires very high resolution to remove. I have been told by a reliable source that the FEM should be less prone to these spurious modes, which would be great.
I have a notebook and a python script at the following links:
https://github.com/francispoulin/firedrakeQG/blob/master/Linear%20Stability%...
https://github.com/francispoulin/firedrakeQG/blob/master/LSA_SW_Jet_Firedrak...
The error comes from the assemble part. I have done this previously with QG and SW in 2D and managed to get it working. This is only a 1D problem, so a little bit different, but should be easier. Unfortunately, when I do try it I get the following error, see below. I have tried different combinations but am not having any success.
One issue is that there are three variables: u, v, eta (height). I define the solution space to be V*V*V where V is a CG space. This is the major difference compared to what I tried before.
I am pretty sure that I am getting the syntax wrong here. Can anyone offer any advice?
Cheers, Francis
(firedrake) fpoulin@fpoulin-Gazelle:~/Research/2015-Ben_Storer/Firedrake$ python LSA_SW_Jet_Firedrake.py ('Ro = ', 0.25) ('Bu = ', 0.0004) ('profile = ', 'bickley') Traceback (most recent call last): File "LSA_SW_Jet_Firedrake.py", line 103, in <module> petsc_a = assemble(a, mat_type='aij', bcs=bc).M.handle File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 102, in assemble allocate_only=allocate_only) File "<decorator-gen-279>", line 2, in _assemble File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/driver.py", line 44, in compile_form fd = ufl_utils.compute_form_data(form) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 387, in compute_form_data check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 152, in check_form_arity check_integrand_arity(itg.integrand(), arguments) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 145, in check_integrand_arity args = map_expr_dag(rules, expr, compress=False) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 37, in map_expr_dag result, = map_expr_dags(function, [expression], compress=compress) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 86, in map_expr_dags r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands]) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 42, in sum raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(a, b)) ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None), Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 1, None)) vs (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None),).
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Thanks Andrew an Julian for your helpful feedback. @Andrew: Yes, you were right that I had improperly defined the form. I have since fixed that and don't get an error. I am now working on getting the correct result. @Julian: This is a very interesting idea that you mention, which I admit I don't quite understand. What do you mean by shifting the boundary conditions out of my spectrum? When I use a matrix method, say spectral collocation on a chebyschev grid, I neglect the rows and columns corresponding to the end points since those are zero. Do you mean something should/can be done with the FEM as well? To make things a bit more complicated, I am only imposing zero Dirichlet BCs on the across channel velocity. If you had a reference you could point me to I will eagerly dig it up and try and better understand your concerns. I certainly want to do this right and avoid spurious eigenmodes if at all possible. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Julian Andrej [juan@tf.uni-kiel.de] Sent: Sunday, March 12, 2017 10:08 AM To: Firedrake Project Subject: Re: [firedrake] Problems with SLEPc You are computing Eigenvalues from a problem which involves Dirichlet boundary conditions but don't shift them out of your spectrum. So your if your Eigenvalues are around 0 real part, this leads to spurious Eigenmodes which arise from Dirichlet values on the boundary. On Sun, Mar 12, 2017 at 2:22 PM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
"Adding expressions with non-matching form arguments {0} vs {1}."
It looks like you declared one of the forms incorrectly, e.g. used a trial function instead of a test function, or tried to add a form with a test function to a form with no test function, etc.
On 12 March 2017 at 12:49, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
I am trying to build on the QG SLEPc demo that I put together a while ago to solve for the linear instability of a 1D Shallow Water Jet. I've done this in the past using spectral methods and have had much success, but I often get spurious modes, which requires very high resolution to remove. I have been told by a reliable source that the FEM should be less prone to these spurious modes, which would be great.
I have a notebook and a python script at the following links:
https://github.com/francispoulin/firedrakeQG/blob/master/Linear%20Stability%...
https://github.com/francispoulin/firedrakeQG/blob/master/LSA_SW_Jet_Firedrak...
The error comes from the assemble part. I have done this previously with QG and SW in 2D and managed to get it working. This is only a 1D problem, so a little bit different, but should be easier. Unfortunately, when I do try it I get the following error, see below. I have tried different combinations but am not having any success.
One issue is that there are three variables: u, v, eta (height). I define the solution space to be V*V*V where V is a CG space. This is the major difference compared to what I tried before.
I am pretty sure that I am getting the syntax wrong here. Can anyone offer any advice?
Cheers, Francis
(firedrake) fpoulin@fpoulin-Gazelle:~/Research/2015-Ben_Storer/Firedrake$ python LSA_SW_Jet_Firedrake.py ('Ro = ', 0.25) ('Bu = ', 0.0004) ('profile = ', 'bickley') Traceback (most recent call last): File "LSA_SW_Jet_Firedrake.py", line 103, in <module> petsc_a = assemble(a, mat_type='aij', bcs=bc).M.handle File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 102, in assemble allocate_only=allocate_only) File "<decorator-gen-279>", line 2, in _assemble File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/driver.py", line 44, in compile_form fd = ufl_utils.compute_form_data(form) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 387, in compute_form_data check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 152, in check_form_arity check_integrand_arity(itg.integrand(), arguments) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 145, in check_integrand_arity args = map_expr_dag(rules, expr, compress=False) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 37, in map_expr_dag result, = map_expr_dags(function, [expression], compress=compress) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 86, in map_expr_dags r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands]) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 42, in sum raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(a, b)) ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None), Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 1, None)) vs (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None),).
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On Sun, Mar 12, 2017 at 3:35 PM, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Thanks Andrew an Julian for your helpful feedback.
@Andrew: Yes, you were right that I had improperly defined the form. I have since fixed that and don't get an error. I am now working on getting the correct result.
@Julian: This is a very interesting idea that you mention, which I admit I don't quite understand. What do you mean by shifting the boundary conditions out of my spectrum?
When I use a matrix method, say spectral collocation on a chebyschev grid, I neglect the rows and columns corresponding to the end points since those are zero. Do you mean something should/can be done with the FEM as well?
Yes, you would have to do something similar. Basically its the same approach since you end up with just a 1 on the diagonal of your matrix where the boundary condition is applied (rows/cols are zero otherwise). This causes SLEPc to find an eigenvalue at 1 (which is obvious and correct numerically, but not physically). So the tasks should be 1. Find the entry where you apply the Dirichlet BC (so the corresponding row/col pair for the node) 2. If you are looking for the smallest eigenvalues add a huge number to the diagonal entry (thats what i meant by shifting out of the spectrum) If you know you only have eigenvalues with large real part, the Dirichlet entries are not a problem imho.
To make things a bit more complicated, I am only imposing zero Dirichlet BCs on the across channel velocity.
If you had a reference you could point me to I will eagerly dig it up and try and better understand your concerns. I certainly want to do this right and avoid spurious eigenmodes if at all possible.
Cheers, Francis
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Julian Andrej [juan@tf.uni-kiel.de] Sent: Sunday, March 12, 2017 10:08 AM To: Firedrake Project Subject: Re: [firedrake] Problems with SLEPc
You are computing Eigenvalues from a problem which involves Dirichlet boundary conditions but don't shift them out of your spectrum. So your if your Eigenvalues are around 0 real part, this leads to spurious Eigenmodes which arise from Dirichlet values on the boundary.
On Sun, Mar 12, 2017 at 2:22 PM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
"Adding expressions with non-matching form arguments {0} vs {1}."
It looks like you declared one of the forms incorrectly, e.g. used a trial function instead of a test function, or tried to add a form with a test function to a form with no test function, etc.
On 12 March 2017 at 12:49, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
I am trying to build on the QG SLEPc demo that I put together a while ago to solve for the linear instability of a 1D Shallow Water Jet. I've done this in the past using spectral methods and have had much success, but I often get spurious modes, which requires very high resolution to remove. I have been told by a reliable source that the FEM should be less prone to these spurious modes, which would be great.
I have a notebook and a python script at the following links:
https://github.com/francispoulin/firedrakeQG/blob/master/Linear%20Stability%...
https://github.com/francispoulin/firedrakeQG/blob/master/LSA_SW_Jet_Firedrak...
The error comes from the assemble part. I have done this previously with QG and SW in 2D and managed to get it working. This is only a 1D problem, so a little bit different, but should be easier. Unfortunately, when I do try it I get the following error, see below. I have tried different combinations but am not having any success.
One issue is that there are three variables: u, v, eta (height). I define the solution space to be V*V*V where V is a CG space. This is the major difference compared to what I tried before.
I am pretty sure that I am getting the syntax wrong here. Can anyone offer any advice?
Cheers, Francis
(firedrake) fpoulin@fpoulin-Gazelle:~/Research/2015-Ben_Storer/Firedrake$ python LSA_SW_Jet_Firedrake.py ('Ro = ', 0.25) ('Bu = ', 0.0004) ('profile = ', 'bickley') Traceback (most recent call last): File "LSA_SW_Jet_Firedrake.py", line 103, in <module> petsc_a = assemble(a, mat_type='aij', bcs=bc).M.handle File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 102, in assemble allocate_only=allocate_only) File "<decorator-gen-279>", line 2, in _assemble File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/driver.py", line 44, in compile_form fd = ufl_utils.compute_form_data(form) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 387, in compute_form_data check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 152, in check_form_arity check_integrand_arity(itg.integrand(), arguments) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 145, in check_integrand_arity args = map_expr_dag(rules, expr, compress=False) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 37, in map_expr_dag result, = map_expr_dags(function, [expression], compress=compress) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 86, in map_expr_dags r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands]) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 42, in sum raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(a, b)) ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None), Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 1, None)) vs (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None),).
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hello everyone, I am happy to say that I have the weak form working and, using lapack, I get virtually the same results as using a spectral method. I believe the FEM converges faster but at this point I'm just happy to get something working. Updated python scripts and notebooks can be found at the following: https://github.com/francispoulin/firedrakeQG/blob/master/LSA_SW_Jet_Firedrak... https://github.com/francispoulin/firedrakeQG/blob/master/Linear%20Stability%... The problem that I'm having is that I want to use an indirect solver it converges to zero, which might be an eigenvalue but the one I want is the eigenvalue with the largest imaginary part. I know that SLEPc can be a bit finicky but does anyone have any advice as to how I might get the direct solver converting to the eigenvalues with a non-zero imaginary party? I remember a few months ago I discovered that cayley seemed to converge better than shift_invert for some problems but I haven't had any success with this problem. opts.setValue("st_type", "cayley") @ Julian: Thanks for the advice. Sadly I'm not sure where to start looking to find the entry for the DirichletBC. I will have to think about that some more. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Julian Andrej [juan@tf.uni-kiel.de] Sent: Sunday, March 12, 2017 10:53 AM To: Firedrake Project Subject: Re: [firedrake] Problems with SLEPc On Sun, Mar 12, 2017 at 3:35 PM, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Thanks Andrew an Julian for your helpful feedback.
@Andrew: Yes, you were right that I had improperly defined the form. I have since fixed that and don't get an error. I am now working on getting the correct result.
@Julian: This is a very interesting idea that you mention, which I admit I don't quite understand. What do you mean by shifting the boundary conditions out of my spectrum?
When I use a matrix method, say spectral collocation on a chebyschev grid, I neglect the rows and columns corresponding to the end points since those are zero. Do you mean something should/can be done with the FEM as well?
Yes, you would have to do something similar. Basically its the same approach since you end up with just a 1 on the diagonal of your matrix where the boundary condition is applied (rows/cols are zero otherwise). This causes SLEPc to find an eigenvalue at 1 (which is obvious and correct numerically, but not physically). So the tasks should be 1. Find the entry where you apply the Dirichlet BC (so the corresponding row/col pair for the node) 2. If you are looking for the smallest eigenvalues add a huge number to the diagonal entry (thats what i meant by shifting out of the spectrum) If you know you only have eigenvalues with large real part, the Dirichlet entries are not a problem imho.
To make things a bit more complicated, I am only imposing zero Dirichlet BCs on the across channel velocity.
If you had a reference you could point me to I will eagerly dig it up and try and better understand your concerns. I certainly want to do this right and avoid spurious eigenmodes if at all possible.
Cheers, Francis
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Julian Andrej [juan@tf.uni-kiel.de] Sent: Sunday, March 12, 2017 10:08 AM To: Firedrake Project Subject: Re: [firedrake] Problems with SLEPc
You are computing Eigenvalues from a problem which involves Dirichlet boundary conditions but don't shift them out of your spectrum. So your if your Eigenvalues are around 0 real part, this leads to spurious Eigenmodes which arise from Dirichlet values on the boundary.
On Sun, Mar 12, 2017 at 2:22 PM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
"Adding expressions with non-matching form arguments {0} vs {1}."
It looks like you declared one of the forms incorrectly, e.g. used a trial function instead of a test function, or tried to add a form with a test function to a form with no test function, etc.
On 12 March 2017 at 12:49, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
I am trying to build on the QG SLEPc demo that I put together a while ago to solve for the linear instability of a 1D Shallow Water Jet. I've done this in the past using spectral methods and have had much success, but I often get spurious modes, which requires very high resolution to remove. I have been told by a reliable source that the FEM should be less prone to these spurious modes, which would be great.
I have a notebook and a python script at the following links:
https://github.com/francispoulin/firedrakeQG/blob/master/Linear%20Stability%...
https://github.com/francispoulin/firedrakeQG/blob/master/LSA_SW_Jet_Firedrak...
The error comes from the assemble part. I have done this previously with QG and SW in 2D and managed to get it working. This is only a 1D problem, so a little bit different, but should be easier. Unfortunately, when I do try it I get the following error, see below. I have tried different combinations but am not having any success.
One issue is that there are three variables: u, v, eta (height). I define the solution space to be V*V*V where V is a CG space. This is the major difference compared to what I tried before.
I am pretty sure that I am getting the syntax wrong here. Can anyone offer any advice?
Cheers, Francis
(firedrake) fpoulin@fpoulin-Gazelle:~/Research/2015-Ben_Storer/Firedrake$ python LSA_SW_Jet_Firedrake.py ('Ro = ', 0.25) ('Bu = ', 0.0004) ('profile = ', 'bickley') Traceback (most recent call last): File "LSA_SW_Jet_Firedrake.py", line 103, in <module> petsc_a = assemble(a, mat_type='aij', bcs=bc).M.handle File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 102, in assemble allocate_only=allocate_only) File "<decorator-gen-279>", line 2, in _assemble File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/driver.py", line 44, in compile_form fd = ufl_utils.compute_form_data(form) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 387, in compute_form_data check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 152, in check_form_arity check_integrand_arity(itg.integrand(), arguments) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 145, in check_integrand_arity args = map_expr_dag(rules, expr, compress=False) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 37, in map_expr_dag result, = map_expr_dags(function, [expression], compress=compress) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 86, in map_expr_dags r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands]) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 42, in sum raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(a, b)) ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None), Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 1, None)) vs (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None),).
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hello again, I am trying to play with the SLEPc options to solve the eigenvalue problem that I asked about yesterday. The direct method works is when I have the following: opts = PETSc.Options() opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") opts.setValue("eps_type", "lapack") opts.setValue("eps_largest_imaginary", None) opts.setValue("eps_tol", 1e-10) I am trying to set up an indirect method and knowing that the eigenvalue for k=0.45 has an imaginary part of 0.11, I would like to give that as a seed. Ideally, I would like to find it without knowing the answer but I am happy to start with a known result for testing purposes. When I try the following, opts = PETSc.Options() opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") #opts.setValue("eps_type", "lapack") opts.setValue("eps_type", "krylovschur") #opts.setValue("eps_largest_imaginary", None) opts.setValue("eps_target_imaginary", None) opts.setValue("eps_target", 0.11) opts.setValue("st_type", "cayley") opts.setValue("eps_tol", 1e-10) I get the the error below. When I ask for largest_imaginary I get zero. Does someone have some advice as to how I might be able to ask for an imaginary target? I know this is similar to a question I asked a few months ago. There the solution was to ask for a magnitude. That worked ok when the spectrum was purely real or imaginary but if it is complex, as it is in this case, it seems less obvious how to do that. Any advice would be greatly appreciated. Cheers, Francis Traceback (most recent call last): File "LSA_SW_Jet_Firedrake.py", line 123, in <module> es.setFromOptions() File "SLEPc/EPS.pyx", line 376, in slepc4py.SLEPc.EPS.setFromOptions (src/slepc4py.SLEPc.c:23442) petsc4py.PETSc.Error: error code 63 [0] EPSSetFromOptions() line 220 in /tmp/pip-0gbT1h-build/src/eps/interface/epsopts.c [0] EPSSetWhichEigenpairs() line 550 in /tmp/pip-0gbT1h-build/src/eps/interface/epsopts.c [0] Argument out of range [0] Invalid 'which' value ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Francis Poulin [fpoulin@uwaterloo.ca] Sent: Monday, March 13, 2017 11:46 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Problems with SLEPc Hello everyone, I am happy to say that I have the weak form working and, using lapack, I get virtually the same results as using a spectral method. I believe the FEM converges faster but at this point I'm just happy to get something working. Updated python scripts and notebooks can be found at the following: https://github.com/francispoulin/firedrakeQG/blob/master/LSA_SW_Jet_Firedrak... https://github.com/francispoulin/firedrakeQG/blob/master/Linear%20Stability%... The problem that I'm having is that I want to use an indirect solver it converges to zero, which might be an eigenvalue but the one I want is the eigenvalue with the largest imaginary part. I know that SLEPc can be a bit finicky but does anyone have any advice as to how I might get the direct solver converting to the eigenvalues with a non-zero imaginary party? I remember a few months ago I discovered that cayley seemed to converge better than shift_invert for some problems but I haven't had any success with this problem. opts.setValue("st_type", "cayley") @ Julian: Thanks for the advice. Sadly I'm not sure where to start looking to find the entry for the DirichletBC. I will have to think about that some more. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Julian Andrej [juan@tf.uni-kiel.de] Sent: Sunday, March 12, 2017 10:53 AM To: Firedrake Project Subject: Re: [firedrake] Problems with SLEPc On Sun, Mar 12, 2017 at 3:35 PM, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Thanks Andrew an Julian for your helpful feedback.
@Andrew: Yes, you were right that I had improperly defined the form. I have since fixed that and don't get an error. I am now working on getting the correct result.
@Julian: This is a very interesting idea that you mention, which I admit I don't quite understand. What do you mean by shifting the boundary conditions out of my spectrum?
When I use a matrix method, say spectral collocation on a chebyschev grid, I neglect the rows and columns corresponding to the end points since those are zero. Do you mean something should/can be done with the FEM as well?
Yes, you would have to do something similar. Basically its the same approach since you end up with just a 1 on the diagonal of your matrix where the boundary condition is applied (rows/cols are zero otherwise). This causes SLEPc to find an eigenvalue at 1 (which is obvious and correct numerically, but not physically). So the tasks should be 1. Find the entry where you apply the Dirichlet BC (so the corresponding row/col pair for the node) 2. If you are looking for the smallest eigenvalues add a huge number to the diagonal entry (thats what i meant by shifting out of the spectrum) If you know you only have eigenvalues with large real part, the Dirichlet entries are not a problem imho.
To make things a bit more complicated, I am only imposing zero Dirichlet BCs on the across channel velocity.
If you had a reference you could point me to I will eagerly dig it up and try and better understand your concerns. I certainly want to do this right and avoid spurious eigenmodes if at all possible.
Cheers, Francis
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Julian Andrej [juan@tf.uni-kiel.de] Sent: Sunday, March 12, 2017 10:08 AM To: Firedrake Project Subject: Re: [firedrake] Problems with SLEPc
You are computing Eigenvalues from a problem which involves Dirichlet boundary conditions but don't shift them out of your spectrum. So your if your Eigenvalues are around 0 real part, this leads to spurious Eigenmodes which arise from Dirichlet values on the boundary.
On Sun, Mar 12, 2017 at 2:22 PM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
"Adding expressions with non-matching form arguments {0} vs {1}."
It looks like you declared one of the forms incorrectly, e.g. used a trial function instead of a test function, or tried to add a form with a test function to a form with no test function, etc.
On 12 March 2017 at 12:49, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
I am trying to build on the QG SLEPc demo that I put together a while ago to solve for the linear instability of a 1D Shallow Water Jet. I've done this in the past using spectral methods and have had much success, but I often get spurious modes, which requires very high resolution to remove. I have been told by a reliable source that the FEM should be less prone to these spurious modes, which would be great.
I have a notebook and a python script at the following links:
https://github.com/francispoulin/firedrakeQG/blob/master/Linear%20Stability%...
https://github.com/francispoulin/firedrakeQG/blob/master/LSA_SW_Jet_Firedrak...
The error comes from the assemble part. I have done this previously with QG and SW in 2D and managed to get it working. This is only a 1D problem, so a little bit different, but should be easier. Unfortunately, when I do try it I get the following error, see below. I have tried different combinations but am not having any success.
One issue is that there are three variables: u, v, eta (height). I define the solution space to be V*V*V where V is a CG space. This is the major difference compared to what I tried before.
I am pretty sure that I am getting the syntax wrong here. Can anyone offer any advice?
Cheers, Francis
(firedrake) fpoulin@fpoulin-Gazelle:~/Research/2015-Ben_Storer/Firedrake$ python LSA_SW_Jet_Firedrake.py ('Ro = ', 0.25) ('Bu = ', 0.0004) ('profile = ', 'bickley') Traceback (most recent call last): File "LSA_SW_Jet_Firedrake.py", line 103, in <module> petsc_a = assemble(a, mat_type='aij', bcs=bc).M.handle File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 102, in assemble allocate_only=allocate_only) File "<decorator-gen-279>", line 2, in _assemble File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/driver.py", line 44, in compile_form fd = ufl_utils.compute_form_data(form) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 387, in compute_form_data check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 152, in check_form_arity check_integrand_arity(itg.integrand(), arguments) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 145, in check_integrand_arity args = map_expr_dag(rules, expr, compress=False) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 37, in map_expr_dag result, = map_expr_dags(function, [expression], compress=compress) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 86, in map_expr_dags r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands]) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 42, in sum raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(a, b)) ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None), Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 1, None)) vs (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None),).
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hello again, Sorry for the bother but I just wanted to follow up and ask again if anyone might have any ideas on how to ask for an imaginary target in SLEPc or converge to eigenvalues with the largest imaginary part in general? I still have not been able to get past the error message from my previous email. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Francis Poulin [fpoulin@uwaterloo.ca] Sent: Tuesday, March 14, 2017 11:43 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Problems with SLEPc Hello again, I am trying to play with the SLEPc options to solve the eigenvalue problem that I asked about yesterday. The direct method works is when I have the following: opts = PETSc.Options() opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") opts.setValue("eps_type", "lapack") opts.setValue("eps_largest_imaginary", None) opts.setValue("eps_tol", 1e-10) I am trying to set up an indirect method and knowing that the eigenvalue for k=0.45 has an imaginary part of 0.11, I would like to give that as a seed. Ideally, I would like to find it without knowing the answer but I am happy to start with a known result for testing purposes. When I try the following, opts = PETSc.Options() opts.setValue("eps_gen_non_hermitian", None) opts.setValue("st_pc_factor_shift_type", "NONZERO") #opts.setValue("eps_type", "lapack") opts.setValue("eps_type", "krylovschur") #opts.setValue("eps_largest_imaginary", None) opts.setValue("eps_target_imaginary", None) opts.setValue("eps_target", 0.11) opts.setValue("st_type", "cayley") opts.setValue("eps_tol", 1e-10) I get the the error below. When I ask for largest_imaginary I get zero. Does someone have some advice as to how I might be able to ask for an imaginary target? I know this is similar to a question I asked a few months ago. There the solution was to ask for a magnitude. That worked ok when the spectrum was purely real or imaginary but if it is complex, as it is in this case, it seems less obvious how to do that. Any advice would be greatly appreciated. Cheers, Francis Traceback (most recent call last): File "LSA_SW_Jet_Firedrake.py", line 123, in <module> es.setFromOptions() File "SLEPc/EPS.pyx", line 376, in slepc4py.SLEPc.EPS.setFromOptions (src/slepc4py.SLEPc.c:23442) petsc4py.PETSc.Error: error code 63 [0] EPSSetFromOptions() line 220 in /tmp/pip-0gbT1h-build/src/eps/interface/epsopts.c [0] EPSSetWhichEigenpairs() line 550 in /tmp/pip-0gbT1h-build/src/eps/interface/epsopts.c [0] Argument out of range [0] Invalid 'which' value ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Francis Poulin [fpoulin@uwaterloo.ca] Sent: Monday, March 13, 2017 11:46 AM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] Problems with SLEPc Hello everyone, I am happy to say that I have the weak form working and, using lapack, I get virtually the same results as using a spectral method. I believe the FEM converges faster but at this point I'm just happy to get something working. Updated python scripts and notebooks can be found at the following: https://github.com/francispoulin/firedrakeQG/blob/master/LSA_SW_Jet_Firedrak... https://github.com/francispoulin/firedrakeQG/blob/master/Linear%20Stability%... The problem that I'm having is that I want to use an indirect solver it converges to zero, which might be an eigenvalue but the one I want is the eigenvalue with the largest imaginary part. I know that SLEPc can be a bit finicky but does anyone have any advice as to how I might get the direct solver converting to the eigenvalues with a non-zero imaginary party? I remember a few months ago I discovered that cayley seemed to converge better than shift_invert for some problems but I haven't had any success with this problem. opts.setValue("st_type", "cayley") @ Julian: Thanks for the advice. Sadly I'm not sure where to start looking to find the entry for the DirichletBC. I will have to think about that some more. Cheers, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Julian Andrej [juan@tf.uni-kiel.de] Sent: Sunday, March 12, 2017 10:53 AM To: Firedrake Project Subject: Re: [firedrake] Problems with SLEPc On Sun, Mar 12, 2017 at 3:35 PM, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Thanks Andrew an Julian for your helpful feedback.
@Andrew: Yes, you were right that I had improperly defined the form. I have since fixed that and don't get an error. I am now working on getting the correct result.
@Julian: This is a very interesting idea that you mention, which I admit I don't quite understand. What do you mean by shifting the boundary conditions out of my spectrum?
When I use a matrix method, say spectral collocation on a chebyschev grid, I neglect the rows and columns corresponding to the end points since those are zero. Do you mean something should/can be done with the FEM as well?
Yes, you would have to do something similar. Basically its the same approach since you end up with just a 1 on the diagonal of your matrix where the boundary condition is applied (rows/cols are zero otherwise). This causes SLEPc to find an eigenvalue at 1 (which is obvious and correct numerically, but not physically). So the tasks should be 1. Find the entry where you apply the Dirichlet BC (so the corresponding row/col pair for the node) 2. If you are looking for the smallest eigenvalues add a huge number to the diagonal entry (thats what i meant by shifting out of the spectrum) If you know you only have eigenvalues with large real part, the Dirichlet entries are not a problem imho.
To make things a bit more complicated, I am only imposing zero Dirichlet BCs on the across channel velocity.
If you had a reference you could point me to I will eagerly dig it up and try and better understand your concerns. I certainly want to do this right and avoid spurious eigenmodes if at all possible.
Cheers, Francis
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
________________________________________ From: firedrake-bounces@imperial.ac.uk [firedrake-bounces@imperial.ac.uk] on behalf of Julian Andrej [juan@tf.uni-kiel.de] Sent: Sunday, March 12, 2017 10:08 AM To: Firedrake Project Subject: Re: [firedrake] Problems with SLEPc
You are computing Eigenvalues from a problem which involves Dirichlet boundary conditions but don't shift them out of your spectrum. So your if your Eigenvalues are around 0 real part, this leads to spurious Eigenmodes which arise from Dirichlet values on the boundary.
On Sun, Mar 12, 2017 at 2:22 PM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
"Adding expressions with non-matching form arguments {0} vs {1}."
It looks like you declared one of the forms incorrectly, e.g. used a trial function instead of a test function, or tried to add a form with a test function to a form with no test function, etc.
On 12 March 2017 at 12:49, Francis Poulin <fpoulin@uwaterloo.ca> wrote:
Hello,
I am trying to build on the QG SLEPc demo that I put together a while ago to solve for the linear instability of a 1D Shallow Water Jet. I've done this in the past using spectral methods and have had much success, but I often get spurious modes, which requires very high resolution to remove. I have been told by a reliable source that the FEM should be less prone to these spurious modes, which would be great.
I have a notebook and a python script at the following links:
https://github.com/francispoulin/firedrakeQG/blob/master/Linear%20Stability%...
https://github.com/francispoulin/firedrakeQG/blob/master/LSA_SW_Jet_Firedrak...
The error comes from the assemble part. I have done this previously with QG and SW in 2D and managed to get it working. This is only a 1D problem, so a little bit different, but should be easier. Unfortunately, when I do try it I get the following error, see below. I have tried different combinations but am not having any success.
One issue is that there are three variables: u, v, eta (height). I define the solution space to be V*V*V where V is a CG space. This is the major difference compared to what I tried before.
I am pretty sure that I am getting the syntax wrong here. Can anyone offer any advice?
Cheers, Francis
(firedrake) fpoulin@fpoulin-Gazelle:~/Research/2015-Ben_Storer/Firedrake$ python LSA_SW_Jet_Firedrake.py ('Ro = ', 0.25) ('Bu = ', 0.0004) ('profile = ', 'bickley') Traceback (most recent call last): File "LSA_SW_Jet_Firedrake.py", line 103, in <module> petsc_a = assemble(a, mat_type='aij', bcs=bc).M.handle File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 102, in assemble allocate_only=allocate_only) File "<decorator-gen-279>", line 2, in _assemble File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/home/fpoulin/software/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/home/fpoulin/software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/driver.py", line 44, in compile_form fd = ufl_utils.compute_form_data(form) File "/home/fpoulin/software/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 387, in compute_form_data check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 152, in check_form_arity check_integrand_arity(itg.integrand(), arguments) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 145, in check_integrand_arity args = map_expr_dag(rules, expr, compress=False) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 37, in map_expr_dag result, = map_expr_dags(function, [expression], compress=compress) File "/home/fpoulin/software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 86, in map_expr_dags r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands]) File "/home/fpoulin/software/firedrake/src/ufl/ufl/algorithms/check_arities.py", line 42, in sum raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(a, b)) ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None), Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 1, None)) vs (Argument(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x7f2a458c5290>, FiniteElement('Lagrange', interval, 2), name=None), Mesh(VectorElement(FiniteElement('Lagrange', interval, 1), dim=1), 1)), 0, None),).
------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo
email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (3)
-
Andrew McRae
-
Francis Poulin
-
Julian Andrej