From 155d82f650ccff88a3e8e7fbbad289fe0ef8e6d3 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 12 Sep 2019 12:37:17 -0500 Subject: [PATCH 1/2] Add robin operator --- pytential/symbolic/pde/scalar.py | 113 +++++++++++++++++++++++++++++++ 1 file changed, 113 insertions(+) diff --git a/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py index 882b1b3c..9d5b275c 100644 --- a/pytential/symbolic/pde/scalar.py +++ b/pytential/symbolic/pde/scalar.py @@ -26,6 +26,7 @@ __doc__ = """ .. autoclass:: L2WeightedPDEOperator .. autoclass:: DirichletOperator .. autoclass:: NeumannOperator +.. autoclass:: RobinOperator """ @@ -308,4 +309,116 @@ class NeumannOperator(L2WeightedPDEOperator): # }}} +# {{{ robin + +class RobinOperator(L2WeightedPDEOperator): + """IE operator and field representation for solving Robin boundary + value problems with scalar kernels a * u + b * u_n = g (e.g. + :class:`sumpy.kernel.LaplaceKernel`) + + .. automethod:: is_unique_only_up_to_constant + .. automethod:: representation + .. automethod:: operator + """ + + def __init__(self, kernel, + loc_sign, a, b, + use_l2_weighting=False, kernel_arguments=None): + """ + :arg a: coefficient for the function value u. + :arg b: coefficient for the normal derivative u_n. + :arg loc_sign: +1 for exterior, -1 for interior + """ + + if not a**2 + b**2 > np.finfo(np.float64).eps: + raise ValueError("a and b cannot be both zero") + + if a * b < -np.finfo(np.float64).eps: + raise ValueError("illposed operator since ab < 0") + + assert loc_sign in [-1, 1] + + from sumpy.kernel import LaplaceKernel + if not isinstance(kernel, LaplaceKernel): + raise NotImplementedError("Unsupported kernel %s" % str(kernel)) + + if kernel_arguments is None: + kernel_arguments = {} + self.kernel_arguments = kernel_arguments + + self.a = a + self.b = b + self.alpha = b / (a**2 + b**2) + self.beta = -1 * a / (a**2 + b**2) + self.loc_sign = loc_sign + + L2WeightedPDEOperator.__init__(self, kernel, use_l2_weighting) + + def is_unique_only_up_to_constant(self): + from sumpy.kernel import LaplaceKernel + laplace_interior_neumann = ( + isinstance(self.kernel, LaplaceKernel) + and self.loc_sign < 0 + and self.a == 0) + laplace_exterior_dirichlet = ( + isinstance(self.kernel, LaplaceKernel) + and self.loc_sign > 0 + and self.b == 0) + return laplace_interior_neumann or laplace_exterior_dirichlet + + def representation(self, u, map_potentials=None, qbx_forced_limit=None, + **kwargs): + sqrt_w = self.get_sqrt_weight() + inv_sqrt_w_u = cse(u/sqrt_w) + + if map_potentials is None: + def map_potentials(x): # pylint:disable=function-redefined + return x + + kwargs["qbx_forced_limit"] = qbx_forced_limit + kwargs["kernel_arguments"] = self.kernel_arguments + + return ( + self.alpha * map_potentials( + sym.S(self.kernel, inv_sqrt_w_u, **kwargs)) + + self.beta * map_potentials( + sym.D(self.kernel, inv_sqrt_w_u, **kwargs)) + ) + + def operator(self, u): + sqrt_w = self.get_sqrt_weight() + inv_sqrt_w_u = cse(u/sqrt_w) + + knl = self.kernel + + knl_kwargs = {} + knl_kwargs["kernel_arguments"] = self.kernel_arguments + + if self.is_unique_only_up_to_constant(): + # The interior Neumann operator in this representation + # has a nullspace. The mean of the density must be matched + # to the desired solution separately. As is, this operator + # returns a mean that is not well-specified. + + amb_dim = self.kernel.dim + ones_contribution = ( + sym.Ones() * sym.mean(amb_dim, amb_dim-1, inv_sqrt_w_u)) + else: + ones_contribution = 0 + + return (-self.loc_sign * 0.5 * u + + sqrt_w * ( + self.a * self.alpha * sym.S( + self.kernel, inv_sqrt_w_u, **knl_kwargs) + + self.a * self.beta * sym.D( + self.kernel, inv_sqrt_w_u, **knl_kwargs) + + self.b * self.alpha * sym.Sp( + knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs) + + self.b * self.beta * sym.Dp( + knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs) + + ones_contribution + )) + +# }}} + # vim: foldmethod=marker -- GitLab From 0ba5354d8e5f8e54469ce0a0c61522ab7eedf771 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 12 Sep 2019 13:25:46 -0500 Subject: [PATCH 2/2] Fix wellposedness tests --- pytential/symbolic/pde/scalar.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py index 9d5b275c..0ddf9baf 100644 --- a/pytential/symbolic/pde/scalar.py +++ b/pytential/symbolic/pde/scalar.py @@ -333,8 +333,9 @@ class RobinOperator(L2WeightedPDEOperator): if not a**2 + b**2 > np.finfo(np.float64).eps: raise ValueError("a and b cannot be both zero") - if a * b < -np.finfo(np.float64).eps: - raise ValueError("illposed operator since ab < 0") + if a * b * loc_sign > np.finfo(np.float64).eps: + raise ValueError("illposed operator with a = %f, b = %f, loc_sign = %d" + % (a, b, loc_sign)) assert loc_sign in [-1, 1] @@ -409,9 +410,9 @@ class RobinOperator(L2WeightedPDEOperator): return (-self.loc_sign * 0.5 * u + sqrt_w * ( self.a * self.alpha * sym.S( - self.kernel, inv_sqrt_w_u, **knl_kwargs) + knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs) + self.a * self.beta * sym.D( - self.kernel, inv_sqrt_w_u, **knl_kwargs) + knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs) + self.b * self.alpha * sym.Sp( knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs) + self.b * self.beta * sym.Dp( -- GitLab