Re: [firedrake] local discontinuity
Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear David, Many thanks for this very detailed explanation. I'll try this way and let you know about the outcome. Best wishes, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56:45 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear David, I tried what you suggested and I have the following error when solving the Laplace equation: ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. I defined: V_d = FunctionSpace(mesh, "CG", 1) V_p = FunctionSpace(mesh, "CG", 1) W = V_d*V_p trial_d, trial_p = TrialFunction(W) v_d, v_p = TestFunction(W) result_mixed = Function(W) u_d, u_p = split(result_mixed) I define the normal pointing outward the disk interface as: n_vec = FacetNormal(mesh) n_int = I_d("+") * n_vec("+") + I_d("-") * n_vec("-") where I_d (and I_p) are the discontinuous (DG0) indicators defined as in Tomasz code. Then I define my weak forms as follows: A_u_d = sigma_d*dot(grad(trial_d),grad(v_d))*avg(I_d)*dx # Laplace in V_d A_u_p = sigma_p*dot(grad(trial_p),grad(v_p))*avg(I_p)*dx # Laplace in V_p L_u_d = sigma_d*v_d*dot(grad(tmp_d),n_int)*4*avg(I_d)*avg(I_p)*dx # interface L_u_p = -sigma_p*v_p*E*ds(domaine_bc_r) - sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx # right BC and interface A_u = A_u_d + A_u_p L_u = L_u_d + L_u_p cell_pb = LinearVariationalProblem(A_u, L_u, result_mixed) cell_solver = LinearVariationalSolver(cell_pb) cell_solver.solve() The full error is pasted below. Do you know where it comes from ? Thank you, Floriane UFL:ERROR Integral of type cell cannot contain a ReferenceNormal. Traceback (most recent call last): File "lin_cell.py", line 196, in <module> cell_solver = LinearVariationalSolver(cell_pb) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 286, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 156, in __init__ pre_function_callback=pre_f_callback) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving_utils.py", line 333, in __init__ form_compiler_parameters=fcp) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 143, in create_assembly_callable collect_loops=True) File "<decorator-gen-280>", line 2, in _assemble File "/Users/mmfg/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/Users/mmfg/firedrake/src/tsfc/tsfc/driver.py", line 46, in compile_form fd = ufl_utils.compute_form_data(form) File "/Users/mmfg/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 382, in compute_form_data _check_facet_geometry(self.integral_data) File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 168, in _check_facet_geometry error("Integral of type %s cannot contain a %s." % (it, cls.__name__)) File "/Users/mmfg/firedrake/src/ufl/ufl/log.py", line 172, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear Floriane, The error you're getting is due to the fact that you have an expression containing facet normals on a cell-wise integral. This one: `sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx` is the trouble-maker. Best regards, - Thomas ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Floriane Gidel [RPG] <mmfg@leeds.ac.uk> Sent: 08 October 2018 14:18:09 To: firedrake Subject: Re: [firedrake] local discontinuity Dear David, I tried what you suggested and I have the following error when solving the Laplace equation: ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. I defined: V_d = FunctionSpace(mesh, "CG", 1) V_p = FunctionSpace(mesh, "CG", 1) W = V_d*V_p trial_d, trial_p = TrialFunction(W) v_d, v_p = TestFunction(W) result_mixed = Function(W) u_d, u_p = split(result_mixed) I define the normal pointing outward the disk interface as: n_vec = FacetNormal(mesh) n_int = I_d("+") * n_vec("+") + I_d("-") * n_vec("-") where I_d (and I_p) are the discontinuous (DG0) indicators defined as in Tomasz code. Then I define my weak forms as follows: A_u_d = sigma_d*dot(grad(trial_d),grad(v_d))*avg(I_d)*dx # Laplace in V_d A_u_p = sigma_p*dot(grad(trial_p),grad(v_p))*avg(I_p)*dx # Laplace in V_p L_u_d = sigma_d*v_d*dot(grad(tmp_d),n_int)*4*avg(I_d)*avg(I_p)*dx # interface L_u_p = -sigma_p*v_p*E*ds(domaine_bc_r) - sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx # right BC and interface A_u = A_u_d + A_u_p L_u = L_u_d + L_u_p cell_pb = LinearVariationalProblem(A_u, L_u, result_mixed) cell_solver = LinearVariationalSolver(cell_pb) cell_solver.solve() The full error is pasted below. Do you know where it comes from ? Thank you, Floriane UFL:ERROR Integral of type cell cannot contain a ReferenceNormal. Traceback (most recent call last): File "lin_cell.py", line 196, in <module> cell_solver = LinearVariationalSolver(cell_pb) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 286, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 156, in __init__ pre_function_callback=pre_f_callback) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving_utils.py", line 333, in __init__ form_compiler_parameters=fcp) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 143, in create_assembly_callable collect_loops=True) File "<decorator-gen-280>", line 2, in _assemble File "/Users/mmfg/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/Users/mmfg/firedrake/src/tsfc/tsfc/driver.py", line 46, in compile_form fd = ufl_utils.compute_form_data(form) File "/Users/mmfg/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 382, in compute_form_data _check_facet_geometry(self.integral_data) File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 168, in _check_facet_geometry error("Integral of type %s cannot contain a %s." % (it, cls.__name__)) File "/Users/mmfg/firedrake/src/ufl/ufl/log.py", line 172, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear Thomas, That should indeed be dS and not dx, thank you! Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Gibson, Thomas <t.gibson15@imperial.ac.uk> Envoyé : lundi 8 octobre 2018 14:34 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, The error you're getting is due to the fact that you have an expression containing facet normals on a cell-wise integral. This one: `sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx` is the trouble-maker. Best regards, - Thomas ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Floriane Gidel [RPG] <mmfg@leeds.ac.uk> Sent: 08 October 2018 14:18:09 To: firedrake Subject: Re: [firedrake] local discontinuity Dear David, I tried what you suggested and I have the following error when solving the Laplace equation: ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. I defined: V_d = FunctionSpace(mesh, "CG", 1) V_p = FunctionSpace(mesh, "CG", 1) W = V_d*V_p trial_d, trial_p = TrialFunction(W) v_d, v_p = TestFunction(W) result_mixed = Function(W) u_d, u_p = split(result_mixed) I define the normal pointing outward the disk interface as: n_vec = FacetNormal(mesh) n_int = I_d("+") * n_vec("+") + I_d("-") * n_vec("-") where I_d (and I_p) are the discontinuous (DG0) indicators defined as in Tomasz code. Then I define my weak forms as follows: A_u_d = sigma_d*dot(grad(trial_d),grad(v_d))*avg(I_d)*dx # Laplace in V_d A_u_p = sigma_p*dot(grad(trial_p),grad(v_p))*avg(I_p)*dx # Laplace in V_p L_u_d = sigma_d*v_d*dot(grad(tmp_d),n_int)*4*avg(I_d)*avg(I_p)*dx # interface L_u_p = -sigma_p*v_p*E*ds(domaine_bc_r) - sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx # right BC and interface A_u = A_u_d + A_u_p L_u = L_u_d + L_u_p cell_pb = LinearVariationalProblem(A_u, L_u, result_mixed) cell_solver = LinearVariationalSolver(cell_pb) cell_solver.solve() The full error is pasted below. Do you know where it comes from ? Thank you, Floriane UFL:ERROR Integral of type cell cannot contain a ReferenceNormal. Traceback (most recent call last): File "lin_cell.py", line 196, in <module> cell_solver = LinearVariationalSolver(cell_pb) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 286, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 156, in __init__ pre_function_callback=pre_f_callback) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving_utils.py", line 333, in __init__ form_compiler_parameters=fcp) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 143, in create_assembly_callable collect_loops=True) File "<decorator-gen-280>", line 2, in _assemble File "/Users/mmfg/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/Users/mmfg/firedrake/src/tsfc/tsfc/driver.py", line 46, in compile_form fd = ufl_utils.compute_form_data(form) File "/Users/mmfg/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 382, in compute_form_data _check_facet_geometry(self.integral_data) File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 168, in _check_facet_geometry error("Integral of type %s cannot contain a %s." % (it, cls.__name__)) File "/Users/mmfg/firedrake/src/ufl/ufl/log.py", line 172, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Note also that using avg in a dx integral doesn’t make much sense (thought it doesn’t actually change the answer) From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Monday, 8 October 2018 at 14:43 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear Thomas, That should indeed be dS and not dx, thank you! Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Gibson, Thomas <t.gibson15@imperial.ac.uk> Envoyé : lundi 8 octobre 2018 14:34 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, The error you're getting is due to the fact that you have an expression containing facet normals on a cell-wise integral. This one: `sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx` is the trouble-maker. Best regards, - Thomas ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Floriane Gidel [RPG] <mmfg@leeds.ac.uk> Sent: 08 October 2018 14:18:09 To: firedrake Subject: Re: [firedrake] local discontinuity Dear David, I tried what you suggested and I have the following error when solving the Laplace equation: ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. I defined: V_d = FunctionSpace(mesh, "CG", 1) V_p = FunctionSpace(mesh, "CG", 1) W = V_d*V_p trial_d, trial_p = TrialFunction(W) v_d, v_p = TestFunction(W) result_mixed = Function(W) u_d, u_p = split(result_mixed) I define the normal pointing outward the disk interface as: n_vec = FacetNormal(mesh) n_int = I_d("+") * n_vec("+") + I_d("-") * n_vec("-") where I_d (and I_p) are the discontinuous (DG0) indicators defined as in Tomasz code. Then I define my weak forms as follows: A_u_d = sigma_d*dot(grad(trial_d),grad(v_d))*avg(I_d)*dx # Laplace in V_d A_u_p = sigma_p*dot(grad(trial_p),grad(v_p))*avg(I_p)*dx # Laplace in V_p L_u_d = sigma_d*v_d*dot(grad(tmp_d),n_int)*4*avg(I_d)*avg(I_p)*dx # interface L_u_p = -sigma_p*v_p*E*ds(domaine_bc_r) - sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx # right BC and interface A_u = A_u_d + A_u_p L_u = L_u_d + L_u_p cell_pb = LinearVariationalProblem(A_u, L_u, result_mixed) cell_solver = LinearVariationalSolver(cell_pb) cell_solver.solve() The full error is pasted below. Do you know where it comes from ? Thank you, Floriane UFL:ERROR Integral of type cell cannot contain a ReferenceNormal. Traceback (most recent call last): File "lin_cell.py", line 196, in <module> cell_solver = LinearVariationalSolver(cell_pb) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 286, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 156, in __init__ pre_function_callback=pre_f_callback) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving_utils.py", line 333, in __init__ form_compiler_parameters=fcp) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 143, in create_assembly_callable collect_loops=True) File "<decorator-gen-280>", line 2, in _assemble File "/Users/mmfg/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/Users/mmfg/firedrake/src/tsfc/tsfc/driver.py", line 46, in compile_form fd = ufl_utils.compute_form_data(form) File "/Users/mmfg/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 382, in compute_form_data _check_facet_geometry(self.integral_data) File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 168, in _check_facet_geometry error("Integral of type %s cannot contain a %s." % (it, cls.__name__)) File "/Users/mmfg/firedrake/src/ufl/ufl/log.py", line 172, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear all, Concerning the integral on the interface between the disk and the plane, is it not equivalent to write *4*avg(I_d)*avg(I_c)*dS and *dS(interface_id) where interface_id is the id of the interface ? I don't get the same results when using one or the other. Thank you, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : lundi 8 octobre 2018 14:48 À : firedrake Objet : Re: [firedrake] local discontinuity Note also that using avg in a dx integral doesn’t make much sense (thought it doesn’t actually change the answer) From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Monday, 8 October 2018 at 14:43 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear Thomas, That should indeed be dS and not dx, thank you! Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Gibson, Thomas <t.gibson15@imperial.ac.uk> Envoyé : lundi 8 octobre 2018 14:34 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, The error you're getting is due to the fact that you have an expression containing facet normals on a cell-wise integral. This one: `sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx` is the trouble-maker. Best regards, - Thomas ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Floriane Gidel [RPG] <mmfg@leeds.ac.uk> Sent: 08 October 2018 14:18:09 To: firedrake Subject: Re: [firedrake] local discontinuity Dear David, I tried what you suggested and I have the following error when solving the Laplace equation: ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. I defined: V_d = FunctionSpace(mesh, "CG", 1) V_p = FunctionSpace(mesh, "CG", 1) W = V_d*V_p trial_d, trial_p = TrialFunction(W) v_d, v_p = TestFunction(W) result_mixed = Function(W) u_d, u_p = split(result_mixed) I define the normal pointing outward the disk interface as: n_vec = FacetNormal(mesh) n_int = I_d("+") * n_vec("+") + I_d("-") * n_vec("-") where I_d (and I_p) are the discontinuous (DG0) indicators defined as in Tomasz code. Then I define my weak forms as follows: A_u_d = sigma_d*dot(grad(trial_d),grad(v_d))*avg(I_d)*dx # Laplace in V_d A_u_p = sigma_p*dot(grad(trial_p),grad(v_p))*avg(I_p)*dx # Laplace in V_p L_u_d = sigma_d*v_d*dot(grad(tmp_d),n_int)*4*avg(I_d)*avg(I_p)*dx # interface L_u_p = -sigma_p*v_p*E*ds(domaine_bc_r) - sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx # right BC and interface A_u = A_u_d + A_u_p L_u = L_u_d + L_u_p cell_pb = LinearVariationalProblem(A_u, L_u, result_mixed) cell_solver = LinearVariationalSolver(cell_pb) cell_solver.solve() The full error is pasted below. Do you know where it comes from ? Thank you, Floriane UFL:ERROR Integral of type cell cannot contain a ReferenceNormal. Traceback (most recent call last): File "lin_cell.py", line 196, in <module> cell_solver = LinearVariationalSolver(cell_pb) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 286, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 156, in __init__ pre_function_callback=pre_f_callback) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving_utils.py", line 333, in __init__ form_compiler_parameters=fcp) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 143, in create_assembly_callable collect_loops=True) File "<decorator-gen-280>", line 2, in _assemble File "/Users/mmfg/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/Users/mmfg/firedrake/src/tsfc/tsfc/driver.py", line 46, in compile_form fd = ufl_utils.compute_form_data(form) File "/Users/mmfg/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 382, in compute_form_data _check_facet_geometry(self.integral_data) File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 168, in _check_facet_geometry error("Integral of type %s cannot contain a %s." % (it, cls.__name__)) File "/Users/mmfg/firedrake/src/ufl/ufl/log.py", line 172, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Yes. Those are equivalent. You should be able to work out which of your results is correct by just integrating the unit function along the line and observing which solution approximates pi*d From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Wednesday, 10 October 2018 at 15:39 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear all, Concerning the integral on the interface between the disk and the plane, is it not equivalent to write *4*avg(I_d)*avg(I_c)*dS and *dS(interface_id) where interface_id is the id of the interface ? I don't get the same results when using one or the other. Thank you, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : lundi 8 octobre 2018 14:48 À : firedrake Objet : Re: [firedrake] local discontinuity Note also that using avg in a dx integral doesn’t make much sense (thought it doesn’t actually change the answer) From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Monday, 8 October 2018 at 14:43 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear Thomas, That should indeed be dS and not dx, thank you! Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Gibson, Thomas <t.gibson15@imperial.ac.uk> Envoyé : lundi 8 octobre 2018 14:34 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, The error you're getting is due to the fact that you have an expression containing facet normals on a cell-wise integral. This one: `sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx` is the trouble-maker. Best regards, - Thomas ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Floriane Gidel [RPG] <mmfg@leeds.ac.uk> Sent: 08 October 2018 14:18:09 To: firedrake Subject: Re: [firedrake] local discontinuity Dear David, I tried what you suggested and I have the following error when solving the Laplace equation: ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. I defined: V_d = FunctionSpace(mesh, "CG", 1) V_p = FunctionSpace(mesh, "CG", 1) W = V_d*V_p trial_d, trial_p = TrialFunction(W) v_d, v_p = TestFunction(W) result_mixed = Function(W) u_d, u_p = split(result_mixed) I define the normal pointing outward the disk interface as: n_vec = FacetNormal(mesh) n_int = I_d("+") * n_vec("+") + I_d("-") * n_vec("-") where I_d (and I_p) are the discontinuous (DG0) indicators defined as in Tomasz code. Then I define my weak forms as follows: A_u_d = sigma_d*dot(grad(trial_d),grad(v_d))*avg(I_d)*dx # Laplace in V_d A_u_p = sigma_p*dot(grad(trial_p),grad(v_p))*avg(I_p)*dx # Laplace in V_p L_u_d = sigma_d*v_d*dot(grad(tmp_d),n_int)*4*avg(I_d)*avg(I_p)*dx # interface L_u_p = -sigma_p*v_p*E*ds(domaine_bc_r) - sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx # right BC and interface A_u = A_u_d + A_u_p L_u = L_u_d + L_u_p cell_pb = LinearVariationalProblem(A_u, L_u, result_mixed) cell_solver = LinearVariationalSolver(cell_pb) cell_solver.solve() The full error is pasted below. Do you know where it comes from ? Thank you, Floriane UFL:ERROR Integral of type cell cannot contain a ReferenceNormal. Traceback (most recent call last): File "lin_cell.py", line 196, in <module> cell_solver = LinearVariationalSolver(cell_pb) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 286, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 156, in __init__ pre_function_callback=pre_f_callback) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving_utils.py", line 333, in __init__ form_compiler_parameters=fcp) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 143, in create_assembly_callable collect_loops=True) File "<decorator-gen-280>", line 2, in _assemble File "/Users/mmfg/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/Users/mmfg/firedrake/src/tsfc/tsfc/driver.py", line 46, in compile_form fd = ufl_utils.compute_form_data(form) File "/Users/mmfg/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 382, in compute_form_data _check_facet_geometry(self.integral_data) File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 168, in _check_facet_geometry error("Integral of type %s cannot contain a %s." % (it, cls.__name__)) File "/Users/mmfg/firedrake/src/ufl/ufl/log.py", line 172, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear all, Thank you, the code seems to work fine, which is also good news for Firedrake since my new team here (at INRIA Bordeaux, France) is now interested to use Firedrake. I'd like to verify my results with an analytical solution of the form u(r, theta) = r*exp(i theta), but I'm having some issues that you might be able to understand. As my mesh is cartesian I define r = sqrt(x[0]*x[0] + x[1]*x[1]), and theta = atan(x[1]/x[0]) (ensuring that x[0]!=0), but have an issue when computing u.interpolate(r*exp(1j*theta)) Traceback (most recent call last): File "sol_analytique.py", line 147, in <module> u.interpolate(r*exp(1j*theta)) TypeError: unsupported operand type(s) for *: 'complex' and 'Function' I also tried u.interpolate(r*exp(complex(0,theta))) But I get a UFL error: Traceback (most recent call last): File "sol_analytique.py", line 147, in <module> u.interpolate(r*exp(complex(0,theta))) File "/Users/mmfg/firedrake/src/ufl/ufl/operators.py", line 584, in exp return _mathfunction(f, Exp) File "/Users/mmfg/firedrake/src/ufl/ufl/operators.py", line 570, in _mathfunction f = as_ufl(f) File "/Users/mmfg/firedrake/src/ufl/ufl/constantvalue.py", line 417, in as_ufl " to any UFL type." % str(expression)) ufl.log.UFLValueError: Invalid type conversion: -1j can not be converted to any UFL type. It might be a very basic question but how can I deal with complex functions ? Thank you, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : mercredi 10 octobre 2018 15:41 À : firedrake Objet : Re: [firedrake] local discontinuity Yes. Those are equivalent. You should be able to work out which of your results is correct by just integrating the unit function along the line and observing which solution approximates pi*d From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Wednesday, 10 October 2018 at 15:39 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear all, Concerning the integral on the interface between the disk and the plane, is it not equivalent to write *4*avg(I_d)*avg(I_c)*dS and *dS(interface_id) where interface_id is the id of the interface ? I don't get the same results when using one or the other. Thank you, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : lundi 8 octobre 2018 14:48 À : firedrake Objet : Re: [firedrake] local discontinuity Note also that using avg in a dx integral doesn’t make much sense (thought it doesn’t actually change the answer) From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Monday, 8 October 2018 at 14:43 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear Thomas, That should indeed be dS and not dx, thank you! Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Gibson, Thomas <t.gibson15@imperial.ac.uk> Envoyé : lundi 8 octobre 2018 14:34 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, The error you're getting is due to the fact that you have an expression containing facet normals on a cell-wise integral. This one: `sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx` is the trouble-maker. Best regards, - Thomas ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Floriane Gidel [RPG] <mmfg@leeds.ac.uk> Sent: 08 October 2018 14:18:09 To: firedrake Subject: Re: [firedrake] local discontinuity Dear David, I tried what you suggested and I have the following error when solving the Laplace equation: ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. I defined: V_d = FunctionSpace(mesh, "CG", 1) V_p = FunctionSpace(mesh, "CG", 1) W = V_d*V_p trial_d, trial_p = TrialFunction(W) v_d, v_p = TestFunction(W) result_mixed = Function(W) u_d, u_p = split(result_mixed) I define the normal pointing outward the disk interface as: n_vec = FacetNormal(mesh) n_int = I_d("+") * n_vec("+") + I_d("-") * n_vec("-") where I_d (and I_p) are the discontinuous (DG0) indicators defined as in Tomasz code. Then I define my weak forms as follows: A_u_d = sigma_d*dot(grad(trial_d),grad(v_d))*avg(I_d)*dx # Laplace in V_d A_u_p = sigma_p*dot(grad(trial_p),grad(v_p))*avg(I_p)*dx # Laplace in V_p L_u_d = sigma_d*v_d*dot(grad(tmp_d),n_int)*4*avg(I_d)*avg(I_p)*dx # interface L_u_p = -sigma_p*v_p*E*ds(domaine_bc_r) - sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx # right BC and interface A_u = A_u_d + A_u_p L_u = L_u_d + L_u_p cell_pb = LinearVariationalProblem(A_u, L_u, result_mixed) cell_solver = LinearVariationalSolver(cell_pb) cell_solver.solve() The full error is pasted below. Do you know where it comes from ? Thank you, Floriane UFL:ERROR Integral of type cell cannot contain a ReferenceNormal. Traceback (most recent call last): File "lin_cell.py", line 196, in <module> cell_solver = LinearVariationalSolver(cell_pb) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 286, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 156, in __init__ pre_function_callback=pre_f_callback) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving_utils.py", line 333, in __init__ form_compiler_parameters=fcp) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 143, in create_assembly_callable collect_loops=True) File "<decorator-gen-280>", line 2, in _assemble File "/Users/mmfg/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/Users/mmfg/firedrake/src/tsfc/tsfc/driver.py", line 46, in compile_form fd = ufl_utils.compute_form_data(form) File "/Users/mmfg/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 382, in compute_form_data _check_facet_geometry(self.integral_data) File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 168, in _check_facet_geometry error("Integral of type %s cannot contain a %s." % (it, cls.__name__)) File "/Users/mmfg/firedrake/src/ufl/ufl/log.py", line 172, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear Floriane, Master Firedrake only works with real values. We are working on complex capabilities, and a beta version is available online (https://github.com/firedrakeproject/firedrake/projects). That said, if your problem is usually real and you’re just trying to use a complex function as a test case, you should probably use a real test case instead as you actually want to test real code paths. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Tuesday, 16 October 2018 at 08:12 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear all, Thank you, the code seems to work fine, which is also good news for Firedrake since my new team here (at INRIA Bordeaux, France) is now interested to use Firedrake. I'd like to verify my results with an analytical solution of the form u(r, theta) = r*exp(i theta), but I'm having some issues that you might be able to understand. As my mesh is cartesian I define r = sqrt(x[0]*x[0] + x[1]*x[1]), and theta = atan(x[1]/x[0]) (ensuring that x[0]!=0), but have an issue when computing u.interpolate(r*exp(1j*theta)) Traceback (most recent call last): File "sol_analytique.py", line 147, in <module> u.interpolate(r*exp(1j*theta)) TypeError: unsupported operand type(s) for *: 'complex' and 'Function' I also tried u.interpolate(r*exp(complex(0,theta))) But I get a UFL error: Traceback (most recent call last): File "sol_analytique.py", line 147, in <module> u.interpolate(r*exp(complex(0,theta))) File "/Users/mmfg/firedrake/src/ufl/ufl/operators.py", line 584, in exp return _mathfunction(f, Exp) File "/Users/mmfg/firedrake/src/ufl/ufl/operators.py", line 570, in _mathfunction f = as_ufl(f) File "/Users/mmfg/firedrake/src/ufl/ufl/constantvalue.py", line 417, in as_ufl " to any UFL type." % str(expression)) ufl.log.UFLValueError: Invalid type conversion: -1j can not be converted to any UFL type. It might be a very basic question but how can I deal with complex functions ? Thank you, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : mercredi 10 octobre 2018 15:41 À : firedrake Objet : Re: [firedrake] local discontinuity Yes. Those are equivalent. You should be able to work out which of your results is correct by just integrating the unit function along the line and observing which solution approximates pi*d From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Wednesday, 10 October 2018 at 15:39 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear all, Concerning the integral on the interface between the disk and the plane, is it not equivalent to write *4*avg(I_d)*avg(I_c)*dS and *dS(interface_id) where interface_id is the id of the interface ? I don't get the same results when using one or the other. Thank you, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : lundi 8 octobre 2018 14:48 À : firedrake Objet : Re: [firedrake] local discontinuity Note also that using avg in a dx integral doesn’t make much sense (thought it doesn’t actually change the answer) From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Monday, 8 October 2018 at 14:43 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear Thomas, That should indeed be dS and not dx, thank you! Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Gibson, Thomas <t.gibson15@imperial.ac.uk> Envoyé : lundi 8 octobre 2018 14:34 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, The error you're getting is due to the fact that you have an expression containing facet normals on a cell-wise integral. This one: `sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx` is the trouble-maker. Best regards, - Thomas ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Floriane Gidel [RPG] <mmfg@leeds.ac.uk> Sent: 08 October 2018 14:18:09 To: firedrake Subject: Re: [firedrake] local discontinuity Dear David, I tried what you suggested and I have the following error when solving the Laplace equation: ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. I defined: V_d = FunctionSpace(mesh, "CG", 1) V_p = FunctionSpace(mesh, "CG", 1) W = V_d*V_p trial_d, trial_p = TrialFunction(W) v_d, v_p = TestFunction(W) result_mixed = Function(W) u_d, u_p = split(result_mixed) I define the normal pointing outward the disk interface as: n_vec = FacetNormal(mesh) n_int = I_d("+") * n_vec("+") + I_d("-") * n_vec("-") where I_d (and I_p) are the discontinuous (DG0) indicators defined as in Tomasz code. Then I define my weak forms as follows: A_u_d = sigma_d*dot(grad(trial_d),grad(v_d))*avg(I_d)*dx # Laplace in V_d A_u_p = sigma_p*dot(grad(trial_p),grad(v_p))*avg(I_p)*dx # Laplace in V_p L_u_d = sigma_d*v_d*dot(grad(tmp_d),n_int)*4*avg(I_d)*avg(I_p)*dx # interface L_u_p = -sigma_p*v_p*E*ds(domaine_bc_r) - sigma_p*v_p*dot(grad(tmp_p),n_int)*4*avg(I_d)*avg(I_p)*dx # right BC and interface A_u = A_u_d + A_u_p L_u = L_u_d + L_u_p cell_pb = LinearVariationalProblem(A_u, L_u, result_mixed) cell_solver = LinearVariationalSolver(cell_pb) cell_solver.solve() The full error is pasted below. Do you know where it comes from ? Thank you, Floriane UFL:ERROR Integral of type cell cannot contain a ReferenceNormal. Traceback (most recent call last): File "lin_cell.py", line 196, in <module> cell_solver = LinearVariationalSolver(cell_pb) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 286, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 156, in __init__ pre_function_callback=pre_f_callback) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving_utils.py", line 333, in __init__ form_compiler_parameters=fcp) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 143, in create_assembly_callable collect_loops=True) File "<decorator-gen-280>", line 2, in _assemble File "/Users/mmfg/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/assemble.py", line 192, in _assemble kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 193, in compile_form number_map).kernels File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__ obj = make_obj() File "/Users/mmfg/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj obj.__init__(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 112, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/Users/mmfg/firedrake/src/tsfc/tsfc/driver.py", line 46, in compile_form fd = ufl_utils.compute_form_data(form) File "/Users/mmfg/firedrake/src/tsfc/tsfc/ufl_utils.py", line 56, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 382, in compute_form_data _check_facet_geometry(self.integral_data) File "/Users/mmfg/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 168, in _check_facet_geometry error("Integral of type %s cannot contain a %s." % (it, cls.__name__)) File "/Users/mmfg/firedrake/src/ufl/ufl/log.py", line 172, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Integral of type cell cannot contain a ReferenceNormal. ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear all, I'd like to extend the problem, below to k disks (k>=1) in one plane. As a reminder, the objective is to solve the Laplace equation in each disk and in the plane, with a discontinuous condition on the jump at the disks interfaces. The case with k=1 works when I solve a WF on a mixed function space with one component being V_d (disk) and the other one being V_p (plane). The Heaviside indicators I_d and I_p in the disk and in the plane enable integration on the disk/plane interface with 4*avg(I_d)*avg(I_p)*dS. I tried various ways to extend the case k=1 to k>1. For simplification; let's assume k=2 (ie 2 disks); - 1st way: I solve the same problem as in the case k=1, except that now my disk indicator function is nonzero in both disks (same with the BCs). - 2nd way: I define my mixed system as Vd*Vd*Vp (2disks, one plane), with an indicator function for each disk, and one for the plane. - 3rd way: I define my mixed system as Vd*Vp where Vd is a VectorFunctionSpace of dim k, and Vp is of dim 1. The disk indicator functions are defined on a VectorFunctionSpace as well; I have tried the 3 methods and none of them converged. Should one of these methods work or should I solve my problem differently ? Many thanks, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear Floriane, Any of those will work to define the problem. I suggest the first is the most scalable solution. Failing to converge can happen for many reasons. You might have a bug which makes your problem ill-posed. Or you might have made the system less well-conditioned by introducing aa second disk. Are you using a direct solver? (I would strongly suggest you do, both because it’s probably faster in 2D and because it’s essentially guaranteed to succeed if your problem is well-posed). Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Monday, 29 October 2018 at 17:11 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear all, I'd like to extend the problem, below to k disks (k>=1) in one plane. As a reminder, the objective is to solve the Laplace equation in each disk and in the plane, with a discontinuous condition on the jump at the disks interfaces. The case with k=1 works when I solve a WF on a mixed function space with one component being V_d (disk) and the other one being V_p (plane). The Heaviside indicators I_d and I_p in the disk and in the plane enable integration on the disk/plane interface with 4*avg(I_d)*avg(I_p)*dS. I tried various ways to extend the case k=1 to k>1. For simplification; let's assume k=2 (ie 2 disks); - 1st way: I solve the same problem as in the case k=1, except that now my disk indicator function is nonzero in both disks (same with the BCs). - 2nd way: I define my mixed system as Vd*Vd*Vp (2disks, one plane), with an indicator function for each disk, and one for the plane. - 3rd way: I define my mixed system as Vd*Vp where Vd is a VectorFunctionSpace of dim k, and Vp is of dim 1. The disk indicator functions are defined on a VectorFunctionSpace as well; I have tried the 3 methods and none of them converged. Should one of these methods work or should I solve my problem differently ? Many thanks, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear David, Thank you for your fast answer. I'll check my first method then and if I cannot find what's wrong I'll let you know. Many thanks, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : lundi 29 octobre 2018 17:26 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, Any of those will work to define the problem. I suggest the first is the most scalable solution. Failing to converge can happen for many reasons. You might have a bug which makes your problem ill-posed. Or you might have made the system less well-conditioned by introducing aa second disk. Are you using a direct solver? (I would strongly suggest you do, both because it’s probably faster in 2D and because it’s essentially guaranteed to succeed if your problem is well-posed). Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Monday, 29 October 2018 at 17:11 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear all, I'd like to extend the problem, below to k disks (k>=1) in one plane. As a reminder, the objective is to solve the Laplace equation in each disk and in the plane, with a discontinuous condition on the jump at the disks interfaces. The case with k=1 works when I solve a WF on a mixed function space with one component being V_d (disk) and the other one being V_p (plane). The Heaviside indicators I_d and I_p in the disk and in the plane enable integration on the disk/plane interface with 4*avg(I_d)*avg(I_p)*dS. I tried various ways to extend the case k=1 to k>1. For simplification; let's assume k=2 (ie 2 disks); - 1st way: I solve the same problem as in the case k=1, except that now my disk indicator function is nonzero in both disks (same with the BCs). - 2nd way: I define my mixed system as Vd*Vd*Vp (2disks, one plane), with an indicator function for each disk, and one for the plane. - 3rd way: I define my mixed system as Vd*Vp where Vd is a VectorFunctionSpace of dim k, and Vp is of dim 1. The disk indicator functions are defined on a VectorFunctionSpace as well; I have tried the 3 methods and none of them converged. Should one of these methods work or should I solve my problem differently ? Many thanks, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear all, I went back to the problem with one disk because I realized it only converges when initialized with the exact solution (in which case I made a test of convergence that confirmed convergence to the exact solution). When I initialise it with u=0 both in the plane and disk, then the solver diverges (DIVERGED_MAX_IT). My solver is defined as follows: V = V_p*V_d u = Function(V) u_d, u_p = split(u) v_d, v_p = TestFunction(V) a = dot(grad(u_p),grad(v_p))*I_p*dx + dot(grad(u_d),grad(v_d))*I_d*dx + (u_d('+')-u_p('+'))*(v_d('+')-v_p('+'))*4*avg(I_d)*avg(I_p)*dS solve(a==0, u, bcs=[bc_left, bc_right, bc_disk_only, bc_plane_only]) where bc_disk_only and bc_plane_only force nodes of V.sub(0) to be zeros in the plane and nodes of V.sub(1) to be zeros in the disk, respectively, and are defined as in Tomasz fluid/structure code, using continuous indicators I_cg_d and I_cg_p (while I_d and I_p are discontinuous). The solver above only converges when u_d and u_p are initialized with the exact analytical solution, or when I replace I_d and I_p in the first two integrals by their continuous forms I_cg_d and I_cg_p (but then the obtained solution is wrong). This is quite a simple equation so I don't know what I'm doing wrong. Does anyone understand what's wrong ? Many thanks, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : lundi 29 octobre 2018 17:26 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, Any of those will work to define the problem. I suggest the first is the most scalable solution. Failing to converge can happen for many reasons. You might have a bug which makes your problem ill-posed. Or you might have made the system less well-conditioned by introducing aa second disk. Are you using a direct solver? (I would strongly suggest you do, both because it’s probably faster in 2D and because it’s essentially guaranteed to succeed if your problem is well-posed). Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Monday, 29 October 2018 at 17:11 To: firedrake <firedrake@imperial.ac.uk> Subject: Re: [firedrake] local discontinuity Dear all, I'd like to extend the problem, below to k disks (k>=1) in one plane. As a reminder, the objective is to solve the Laplace equation in each disk and in the plane, with a discontinuous condition on the jump at the disks interfaces. The case with k=1 works when I solve a WF on a mixed function space with one component being V_d (disk) and the other one being V_p (plane). The Heaviside indicators I_d and I_p in the disk and in the plane enable integration on the disk/plane interface with 4*avg(I_d)*avg(I_p)*dS. I tried various ways to extend the case k=1 to k>1. For simplification; let's assume k=2 (ie 2 disks); - 1st way: I solve the same problem as in the case k=1, except that now my disk indicator function is nonzero in both disks (same with the BCs). - 2nd way: I define my mixed system as Vd*Vd*Vp (2disks, one plane), with an indicator function for each disk, and one for the plane. - 3rd way: I define my mixed system as Vd*Vp where Vd is a VectorFunctionSpace of dim k, and Vp is of dim 1. The disk indicator functions are defined on a VectorFunctionSpace as well; I have tried the 3 methods and none of them converged. Should one of these methods work or should I solve my problem differently ? Many thanks, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Ham, David A <david.ham@imperial.ac.uk> Envoyé : vendredi 5 octobre 2018 15:56 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane, I think you can do this in a manner analogous to the way that Thomasz did fluid-structure interaction. We define a mixed function space with two continuous components: (V_d, V_p). We’re going to use the first component for the solution in the disk and the second component for the solution outside the disk. We further define two DG0 functions I_d and I_p such that I_d is 1 inside the disk and 0 outside, and I_p = 1 – I_d. You can now write the Laplace equation for the two parts using essentially normal Firedrake code except that you multiply the test function by the appropriate indicator function. This ensures that you only actually assemble integral contributions on the correct part of the domain. Now we just need the surface integral over the facets on the edge of the disk. Observe that this is the only place in the domain where both indicator functions are positive. You can write your jump integrals using the *dS measure. You can restrict the intervals to the relevant edges by multiplying the integrand by 4*avg(I_p)*avg(I_d). If W = V_d * V_p and u_d, u_p = TrialFunctions(W) then you can write the jump as: avg(u_d) – avg(u_p) The avg is mathematically unnecessary and does nothing, however in UFL all terminals that appear in dS integrals have to be restricted (you could just as well write u_d(‘+’) – u_p(‘+’)). Finally, when you come to the solve, you need to ensure that all the V_d nodes outside the disk are eliminated, and all the V_p nodes inside the disk too. You do this using the same DirichletBC subclassing trick which is in https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py... Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. Linear mixed fluid-structure interaction system ...<https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html> www.firedrakeproject.org Linear mixed fluid-structure interaction system¶. This tutorial demonstrates the use of subdomain functionality and show how to describe a system consisting of multiple materials in Firedrake. I realise that is slightly involved, but it should work and the syntax will be fairly clean once you have done it. Feel free to ask for clarifications. Regards, David From: <firedrake-bounces@imperial.ac.uk> on behalf of "Floriane Gidel [RPG]" <mmfg@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Friday, 5 October 2018 at 15:33 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] local discontinuity Dear Firedrakers, I would like to solve equations in a domain made of two subdomains: - a disk, inside which test functions are continuous (Omega_d) - a plane around the disk, in which test functions are continuous (Omega_p) However, I would like the test functions to be discontinuous across the interface Gamma between the disk and the surrounding plane. This is to solve Laplace equation in each domain, augmented by a dynamic boundary condition on the jump between the two domains, that is Delta u = 0 in Omega_d Delta u = 0 in Omega_p partial_t [u] = ... on Gamma, where [u]= u(+)-u(-) is the jump Is there a way to define a functionspace on the full mesh, with continuous basis functions in each subdomain but discontinuous across the interface ? Or is there any better way to solve this kind of system ? Thank you for your help, Floriane
Dear Floriane,
On 5 Nov 2018, at 09:18, Floriane Gidel [RPG] <mmfg@leeds.ac.uk> wrote:
Dear all,
I went back to the problem with one disk because I realized it only converges when initialized with the exact solution (in which case I made a test of convergence that confirmed convergence to the exact solution). When I initialise it with u=0 both in the plane and disk, then the solver diverges (DIVERGED_MAX_IT). My solver is defined as follows:
V = V_p*V_d u = Function(V) u_d, u_p = split(u) v_d, v_p = TestFunction(V)
a = dot(grad(u_p),grad(v_p))*I_p*dx + dot(grad(u_d),grad(v_d))*I_d*dx + (u_d('+')-u_p('+'))*(v_d('+')-v_p('+'))*4*avg(I_d)*avg(I_p)*dS solve(a==0, u, bcs=[bc_left, bc_right, bc_disk_only, bc_plane_only])
where bc_disk_only and bc_plane_only force nodes of V.sub(0) to be zeros in the plane and nodes of V.sub(1) to be zeros in the disk, respectively, and are defined as in Tomasz fluid/structure code, using continuous indicators I_cg_d and I_cg_p (while I_d and I_p are discontinuous). The solver above only converges when u_d and u_p are initialized with the exact analytical solution, or when I replace I_d and I_p in the first two integrals by their continuous forms I_cg_d and I_cg_p (but then the obtained solution is wrong).
This is quite a simple equation so I don't know what I'm doing wrong. Does anyone understand what's wrong ?
Let's start by confirming that you can solve this equation with an exact (direct) solver. If you say: parameters = {"snes_type": "ksponly", "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps", "ksp_converged_reason": True, "snes_converged_reason": True} and run with: solve(..., solver_parameters=parameters) can you obtain a solution? Cheers, Lawrence
Dear Lawrence, No, I get the following error: Traceback (most recent call last): File "main.py", line 103, in <module> solve(a==0,u,bcs=[bc_left, bc_right, bc_disk_only,bc_plane_only], solver_parameter=parameters) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving.py", line 124, in solve _solve_varproblem(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving.py", line 172, in _solve_varproblem solver.solve() File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 220, in solve self.snes.solve(None, v) File "PETSc/SNES.pyx", line 537, in petsc4py.PETSc.SNES.solve (src/petsc4py.PETSc.c:172544) petsc4py.PETSc.Error: error code 56 [0] SNESSolve() line 3967 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/snes/interface/snes.c [0] SNESSolve_KSPONLY() line 38 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/snes/impls/ksponly/ksponly.c [0] KSPSolve() line 598 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/ksp/interface/itfunc.c [0] KSPSetUp() line 378 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 924 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 84 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 250 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 17 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Lawrence Mitchell <wence@gmx.li> Envoyé : lundi 5 novembre 2018 09:22 À : firedrake Objet : Re: [firedrake] local discontinuity Dear Floriane,
On 5 Nov 2018, at 09:18, Floriane Gidel [RPG] <mmfg@leeds.ac.uk> wrote:
Dear all,
I went back to the problem with one disk because I realized it only converges when initialized with the exact solution (in which case I made a test of convergence that confirmed convergence to the exact solution). When I initialise it with u=0 both in the plane and disk, then the solver diverges (DIVERGED_MAX_IT). My solver is defined as follows:
V = V_p*V_d u = Function(V) u_d, u_p = split(u) v_d, v_p = TestFunction(V)
a = dot(grad(u_p),grad(v_p))*I_p*dx + dot(grad(u_d),grad(v_d))*I_d*dx + (u_d('+')-u_p('+'))*(v_d('+')-v_p('+'))*4*avg(I_d)*avg(I_p)*dS solve(a==0, u, bcs=[bc_left, bc_right, bc_disk_only, bc_plane_only])
where bc_disk_only and bc_plane_only force nodes of V.sub(0) to be zeros in the plane and nodes of V.sub(1) to be zeros in the disk, respectively, and are defined as in Tomasz fluid/structure code, using continuous indicators I_cg_d and I_cg_p (while I_d and I_p are discontinuous). The solver above only converges when u_d and u_p are initialized with the exact analytical solution, or when I replace I_d and I_p in the first two integrals by their continuous forms I_cg_d and I_cg_p (but then the obtained solution is wrong).
This is quite a simple equation so I don't know what I'm doing wrong. Does anyone understand what's wrong ?
Let's start by confirming that you can solve this equation with an exact (direct) solver. If you say: parameters = {"snes_type": "ksponly", "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps", "ksp_converged_reason": True, "snes_converged_reason": True} and run with: solve(..., solver_parameters=parameters) can you obtain a solution? Cheers, Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Ah, sorry!
On 5 Nov 2018, at 09:30, Floriane Gidel [RPG] <mmfg@leeds.ac.uk> wrote:
Dear Lawrence,
No, I get the following error:
Traceback (most recent call last): File "main.py", line 103, in <module> solve(a==0,u,bcs=[bc_left, bc_right, bc_disk_only,bc_plane_only], solver_parameter=parameters) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving.py", line 124, in solve _solve_varproblem(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving.py", line 172, in _solve_varproblem solver.solve() File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 220, in solve self.snes.solve(None, v) File "PETSc/SNES.pyx", line 537, in petsc4py.PETSc.SNES.solve (src/petsc4py.PETSc.c:172544) petsc4py.PETSc.Error: error code 56 [0] SNESSolve() line 3967 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/snes/interface/snes.c [0] SNESSolve_KSPONLY() line 38 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/snes/impls/ksponly/ksponly.c [0] KSPSolve() line 598 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/ksp/interface/itfunc.c [0] KSPSetUp() line 378 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 924 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 84 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 250 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 17 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest
I forgot you had a mixed space. Please additionally add the option "mat_type": "aij" to the solver parameters, and try again. Thanks, Lawrence
Yes, that works ! Linear firedrake_0_ solve converged due to CONVERGED_ITS iterations 1 Nonlinear firedrake_0_ solve converged due to CONVERGED_ITS iterations 1 David also suggested to use a direct solver but I did not use the option "mat_type": "aij". I'll read the doc about the options you suggested so I (hopefully) don't need to ask next time. Many many thanks, also for answering so fast. I'll do the test of convergence again to make sure it converges to the exact solution, but I expect it still does. Thank you very much, Floriane ________________________________ De : firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> de la part de Lawrence Mitchell <wence@gmx.li> Envoyé : lundi 5 novembre 2018 09:50 À : firedrake@imperial.ac.uk Objet : Re: [firedrake] local discontinuity Ah, sorry!
On 5 Nov 2018, at 09:30, Floriane Gidel [RPG] <mmfg@leeds.ac.uk> wrote:
Dear Lawrence,
No, I get the following error:
Traceback (most recent call last): File "main.py", line 103, in <module> solve(a==0,u,bcs=[bc_left, bc_right, bc_disk_only,bc_plane_only], solver_parameter=parameters) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving.py", line 124, in solve _solve_varproblem(*args, **kwargs) File "/Users/mmfg/firedrake/src/firedrake/firedrake/solving.py", line 172, in _solve_varproblem solver.solve() File "/Users/mmfg/firedrake/src/firedrake/firedrake/variational_solver.py", line 220, in solve self.snes.solve(None, v) File "PETSc/SNES.pyx", line 537, in petsc4py.PETSc.SNES.solve (src/petsc4py.PETSc.c:172544) petsc4py.PETSc.Error: error code 56 [0] SNESSolve() line 3967 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/snes/interface/snes.c [0] SNESSolve_KSPONLY() line 38 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/snes/impls/ksponly/ksponly.c [0] KSPSolve() line 598 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/ksp/interface/itfunc.c [0] KSPSetUp() line 378 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 924 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/pc/interface/precon.c [0] PCSetUp_LU() line 84 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/ksp/pc/impls/factor/lu/lu.c [0] MatGetOrdering() line 250 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/mat/order/sorder.c [0] MatGetOrdering_ND() line 17 in /private/var/folders/7h/wbj8xp7n3g5cfbr32ctcmwzcy3jf53/T/pip-yGXqPh-build/src/mat/order/spnd.c [0] No support for this operation for this object type [0] Cannot get rows for matrix type nest
I forgot you had a mixed space. Please additionally add the option "mat_type": "aij" to the solver parameters, and try again. Thanks, Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 5 Nov 2018, at 10:14, Floriane Gidel [RPG] <mmfg@leeds.ac.uk> wrote:
Yes, that works !
Linear firedrake_0_ solve converged due to CONVERGED_ITS iterations 1 Nonlinear firedrake_0_ solve converged due to CONVERGED_ITS iterations 1
David also suggested to use a direct solver but I did not use the option "mat_type": "aij". I'll read the doc about the options you suggested so I (hopefully) don't need to ask next time. Many many thanks, also for answering so fast. I'll do the test of convergence again to make sure it converges to the exact solution, but I expect it still does.
OK, so this suggests that the iterative solver you get by default for this problem is rather weak. It's "just" a laplacian with a discontinuous coefficient and one internal jump condition. So it may well be that you can use: parameters = {"mat_type": "aij", "pc_type": "hypre", "ksp_type": "cg", "snes_type": "ksponly"} You should also run with "ksp_monitor": True while testing to ensure that you converge in a small number of iterations. If you're only going to do 2D problems, depending on the size, up to around 100000 dofs in total is probably fastest with a direct solver. Cheers, Lawrence
participants (5)
-
Floriane Gidel [RPG]
-
Gibson, Thomas
-
Ham, David A
-
Lawrence Mitchell
-
Lawrence Mitchell