From 4638de23e20903e8c5540b0d8e12abdc7008a56d Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 1 Jun 2017 18:32:07 -0400 Subject: [PATCH 01/11] adding volume potential --- examples/cahn-hilliard.py | 297 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 283 insertions(+), 14 deletions(-) diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index f7cd6559..1f589405 100644 --- a/examples/cahn-hilliard.py +++ b/examples/cahn-hilliard.py @@ -18,6 +18,9 @@ qbx_order = bdry_quad_order bdry_ovsmp_quad_order = 4*bdry_quad_order fmm_order = 8 +vol_quad_order = 5 +vol_qbx_order = 2 + # }}} @@ -28,37 +31,303 @@ def main(): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - from meshmode.mesh.generation import ellipse, make_curve_mesh - from functools import partial +# {{{ volume mesh generation + + from enum import Enum + class Geometry(Enum): + RegularRectangle = 1 + Circle = 2 + + shape = Geometry.Circle + + if shape == Geometry.RegularRectangle: + from meshmode.mesh.generation import generate_regular_rect_mesh + ext = 1. + h = 0.05 + mesh = generate_regular_rect_mesh( + a=(-ext/2., -ext/2.), + b=( ext/2., ext/2.), + n=(int(ext/h), int(ext/h)) + ) + elif shape == Geometry.Circle: + from meshmode.mesh.io import import generate_gmsh + h = 0.05 + mesh = generate_emsh( + FileSource("circle.step"), + 2, + order=mesh_order, + force_ambient_dim=2, + orther_options=["-string", + "Mesh.CharacteristicLengthMax = %g;" % h] + ) + else: + 1/0 + + logger.info("%d elements" % mesh.nelements) + +# }}} + +# {{{ discretization and connections - mesh = make_curve_mesh( - partial(ellipse, 2), - np.linspace(0, 1, nelements+1), - mesh_order) + vol_discr = Discretization(ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(vol_quad_order)) - pre_density_discr = Discretization( - cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order)) + from meshmode.mesh import BTAG_ALL + from meshmode.discretization.connection import make_face_restriction + pre_density_connection = make_face_restriction( + vol_discr, + InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order), + BTAG_ALL) + pre_density_discr = pre_density_connection.to_discr from pytential.qbx import ( QBXLayerPotentialSource, QBXTargetAssociationFailedException) + qbx, _ = QBXLayerPotentialSource( - pre_density_discr, fine_order=bdry_ovsmp_quad_order, qbx_order=qbx_order, + pre_density_discr, + fine_order=bdry_ovsmp_quad_order, qbx_order=qbx_order, fmm_order=fmm_order, expansion_disks_in_tree_have_extent=True, ).with_refinement() + density_discr = qbx.density_discr + # composition of connetions + # vol_discr --> pre_density_discr --> density_discr + # via ChainedDiscretizationConnection + + +# from meshmode.mesh.generation import ellipse, make_curve_mesh +# from functools import partial +# +# mesh = make_curve_mesh( +# partial(ellipse, 2), +# np.linspace(0, 1, nelements+1), +# mesh_order) +# +# pre_density_discr = Discretization( +# cl_ctx, mesh, +# InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order)) +# +# from pytential.qbx import ( +# QBXLayerPotentialSource, QBXTargetAssociationFailedException) +# qbx, _ = QBXLayerPotentialSource( +# pre_density_discr, fine_order=bdry_ovsmp_quad_order, qbx_order=qbx_order, +# fmm_order=fmm_order, +# expansion_disks_in_tree_have_extent=True, +# ).with_refinement() +# density_discr = qbx.density_discr + +# }}} + +# {{{ a kernel class for G0 +# FIXME: will the expressions work when lambda is complex? +# (may need conversion from 1j to var("I")) + class ShidongKernel(ExpressionKernel): + init_arg_names = ("dim", "lambda1", "lambda2") + + def __init__(self, dim=None, lambda1=0., lambda2=1.): + """ + :arg lambda1,lambda2: The roots of the quadratic equation w.r.t + laplacian. + """ + # Assert against repeated roots. + if abs(lambda1**2 - lambda2**2) < 1e-9: + raise RuntimeError("illposed since input roots are too close") + +# Based on http://mathworld.wolfram.com/ModifiedBesselFunctionoftheSecondKind.html + if dim == 2: + r = pymbolic_real_norm_2(make_sym_vector("d", dim)) + expr = var("hankel_1")(0, var("I") * lambda1 * r) - \ + var("hankel_1")(0, var("I") * lambda2 * r) + scaling = 1. / ( 4. * var("I") * (lambda1**2 - lambda2**2) ) + else: + raise RuntimeError("unsupported dimensionality") + + ExpressionKernel.__init__( + self, + dim, + expression=expr, + scaling = scaling, + is_complex_valued=True) + + self.lambda1 = lambda1 + self.lambda2 = lambda2 + + def __getinitargs__(self): + return(self._dim, self.lambda1, self.lambda2) + + def update_persistent_hash(self, key_hash, key_builder): + key_hash.update(type(self).__name__.encode("utf8")) + key_builder.rec(key_hash, + (self._dim, self.lambda1, self.lambda2) + ) + + def __repr__(self): + if self._dim is not None: + return "ShdgKnl%dD(%f, %f)" % ( + self._dim, self.lambda1, self.lambda2) + else: + return "ShdgKnl(%f, %f)" % (self.lambda1, self.lambda2) + + def prepare_loopy_kernel(self, loopy_knl): + from sumpy.codegen import (bessel_preamble_generator, bessel_mangler) + loopy_knl = lp.register_function_manglers(loopy_knl, + [bessel_mangler]) + loopy_knl = lp.register_preamble_generators(loopy_knl, + [bessel_preamble_generator]) + return loopy_knl + + def get_args(self): + k_dtype = np.complex128 + return [ + KernelArgument( + loopy_arg=lp.ValueArg("shidong_kernel", k_dtype), + )] + + mapper_method = "map_shidong_kernel" + +# }}} + +# {{{ extended kernel getters + def get_extkernel_for_G0(lambda1, lambda2): + from sumpy.symbolic import pymbolic_real_norm_2 + from pymbolic.primitives import make_sym_vector + from pymbolic import var + + d = make_sym_vector("d", 3) + r2 = pymbolic_real_norm_2(d[:-1]) + expr = var("hankel_1")(0, var("I") * lambda1 * r2 + + var("I") * d[-1]**2) \ + - var("hankel_1")(0, var("I") * lambda2 * r2 + + var("I") * d[-1]**2) + scaling = 1. / ( 4. * var("I") * (lambda1**2 - lambda2**2) ) + + from sumpy.kernel import ExpressionKernel + return ExpressionKernel( + dim=3, + expression=expr, + scaling=scaling, + is_complex_valued=True) + + def get_extkernel_for_G1(lamb): + from sumpy.symbolic import pymbolic_real_norm_2 + from pymbolic.primitives import make_sym_vector + from pymbolic import var + + d = make_sym_vector("d", 3) + r2 = pymbolic_real_norm_2(d[:-1]) + expr = var("hankel_1")(0, var("I") * lamb * r2 + + var("I") * d[-1]**2) + scaling = - var("I") / 4. + + from sumpy.kernel import ExpressionKernel + return ExpressionKernel( + dim=3, + expression=expr, + scaling=scaling, + is_complex_valued=True) + + +# }}} + +# {{{ equation info + + s = 1.5 + epsilon = 0.01 + delta_t = 0.05 + + b = s / (epsilon**2) + c = 1. / (epsilon * delta_t) + + sqdet = np.sqrt( b**2 - 4. * c ) + assert np.abs(sqdet) > 1e-6 + + lambda1 = ( b + sqdet ) / 2. + lambda2 = ( b - sqdet ) / 2. + from pytential.symbolic.pde.cahn_hilliard import CahnHilliardOperator chop = CahnHilliardOperator( # FIXME: Constants? - lambda1=1.5, - lambda2=1.25, - c=1) + lambda1=lambda1, + lambda2=lambda2, + c=c) unk = chop.make_unknown("sigma") bound_op = bind(qbx, chop.operator(unk)) + yukawa_2d_in_3d_kernel_1 = get_extkernel_for_G1(lambda1) + shidong_2d_in_3d_kernel = get_extkernel_for_G0(lambda1, lambda2) + + from sumpy.qbx import LayerPotential + from sumpy.qbx import LineTaylorLocalExpansion + layer_pot_v0f1 = LayerPotential(ctx, [ + LineTaylorLocalExpansion(shidong_2d_in_3d_kernel, + order=vol_qbx_order)]) + layer_pot_v1f1 = LayerPotential(ctx, [ + LineTaylorLocalExpansion(yukawa_2d_in_3d_kernel_1, + order=vol_qbx_order)]) +# }}} + +# {{{ volume integral + + vol_x = vol_discr.nodes(),with_queue(queue) + + # target points + targets = cl.array.zeros(queue, (3,) + vol_x.shape[1:], vol_x.dtype) + targets[:2] = vol_x + + # expansion centers + center_dist = 0.125*np.min( + cl.clmath.sqrt( + bind(vol_discr, + p.area_element(mesh.ambient_dim, mesh.dim)) + (queue)).get()) + centers = make_obj_array([ci.copy().reshape(vol_discr.nnodes) for ci in targets]) + centers[2][:] = center_dist + print(center_dist) + + # source points + # FIXME: use over sampled source points? + sources = cl.array.zeros(queue, (3,) + vol_x.shape[1:], vol_x.dtype) + sources[:2] = vol_x + + # a manufactured f1 + x_sin_factor = 30 + y_sin_factor = 10 + def f1_func(x, y) + return 0.1 * cl.clmath.sin(x_sin_factor*x) * cl.clmath.sin(y_sin_factor*y) + + # strengths (with quadrature weights) + f1 = f1_func(vol_x[0], vol_x[1]) + vol_weights = bind(vol_discr, + p.area_element(mesh.ambient_dim, mesh.dim) * p.QWeight() + )(queue) + + print("volume: %d source nodes, %d target nodes" % ( + vol_discr.nnodes, vol_discr.nnodes)) + + evt, (vol_pot_v0f1,) = layer_pot_v0f1( + queue, + targets=targets.reshape(3, vol_discr.nnodes), + centers=centers, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=( + (vol_weights * f1).reshape(vol_discr.nnodes),) + ) + + evt, (vol_pot_v1f1,) = layer_pot_v1f1( + queue, + targets=targets.reshape(3, vol_discr.nnodes), + centers=centers, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=( + (vol_weights * f1).reshape(vol_discr.nnodes),) + ) + +# }}} + + # {{{ fix rhs and solve nodes = density_discr.nodes().with_queue(queue) @@ -101,7 +370,7 @@ def main(): indicator = bind( (indicator_qbx, PointsTarget(targets)), sym.D(LaplaceKernel(2), sym.var("sigma")))( - queue, sigma=ones_density).get() + queue, sigma=ones_density).get() try: fld_in_vol = bind( -- GitLab From 1128669a22ef31f8a54c03a5641af3adbaddc28c Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 1 Jun 2017 23:08:35 -0400 Subject: [PATCH 02/11] minor fixes --- examples/cahn-hilliard.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index 16f6091a..1495bddb 100644 --- a/examples/cahn-hilliard.py +++ b/examples/cahn-hilliard.py @@ -50,9 +50,9 @@ def main(): n=(int(ext/h), int(ext/h)) ) elif shape == Geometry.Circle: - from meshmode.mesh.io import import generate_gmsh + from meshmode.mesh.io import generate_gmsh, FileSource h = 0.05 - mesh = generate_emsh( + mesh = generate_gmsh( FileSource("circle.step"), 2, order=mesh_order, @@ -291,7 +291,7 @@ def main(): # a manufactured f1 x_sin_factor = 30 y_sin_factor = 10 - def f1_func(x, y) + def f1_func(x, y): return 0.1 * cl.clmath.sin(x_sin_factor*x) * cl.clmath.sin(y_sin_factor*y) # strengths (with quadrature weights) -- GitLab From d02d76fcdc3d2862cb41d7ec9082347605b7f539 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 2 Jun 2017 10:17:02 -0400 Subject: [PATCH 03/11] improved root finding algorithm --- pytential/symbolic/pde/cahn_hilliard.py | 33 +++++++++++++++++++++++-- 1 file changed, 31 insertions(+), 2 deletions(-) diff --git a/pytential/symbolic/pde/cahn_hilliard.py b/pytential/symbolic/pde/cahn_hilliard.py index 459f1993..3bbf1bf0 100644 --- a/pytential/symbolic/pde/cahn_hilliard.py +++ b/pytential/symbolic/pde/cahn_hilliard.py @@ -38,12 +38,41 @@ class CahnHilliardOperator(L2WeightedPDEOperator): self.b = b self.c = c - lam1 = np.sqrt(-(np.sqrt(b**2-4*c)-b)/2) - lam2 = np.sqrt((np.sqrt(b**2-4*c)+b)/2) + # Issue: + # - let crat = np.abs(4.*c) / ( b**2 + 1e-12 ) + # - when crat is close to zero, sqrt(b**2-4*c) is close to abs(b), then for b>=0, + # sqrt(b**2-4*c) - b is inaccurate. + # - similarly, when crat is close to one, sqrt(b**2-4*c) is close to zero, then for b>0, + # sqrt(b**2-4*c) + b is inaccurate. + # Solution: + # - find a criteria for crat to choose from formulae, or + # - use the computed root with smaller residual + def quadratic_formula_1(a, b, c): + return ( -b + np.sqrt(b**2 - 4*a*c) ) / (2*a) + def quadratic_formula_2(a, b, c): + return ( -b - np.sqrt(b**2 - 4*a*c) ) / (2*a) + def citardauq_formula_1(a, b, c): + return 2*c / (-b - np.sqrt(b**2-4*a*c)) + def citardauq_formula_2(a, b, c): + return 2*c / (-b + np.sqrt(b**2-4*a*c)) def f(x): return x**2 - b*x + c + root11 = quadratic_formula_1(1, -b, c) + root12 = citardauq_formula_1(1, -b, c) + if np.abs(f(root11)) < np.abs(f(root12)): + lam1 = np.sqrt(root11) + else: + lam1 = np.sqrt(root12) + + root21 = quadratic_formula_2(1, -b, c) + root22 = citardauq_formula_2(1, -b, c) + if np.abs(f(root21)) < np.abs(f(root22)): + lam1 = np.sqrt(root21) + else: + lam2 = np.sqrt(root22) + assert np.abs(f(lam1**2)) < 1e-12 assert np.abs(f(lam2**2)) < 1e-12 -- GitLab From 495deb6d75f292ae290561a31f0cef1af731082f Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 2 Jun 2017 10:17:23 -0400 Subject: [PATCH 04/11] more fixes --- examples/cahn-hilliard.py | 237 +++++++++++++++++++------------------- 1 file changed, 119 insertions(+), 118 deletions(-) diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index 1495bddb..7c8e9ff8 100644 --- a/examples/cahn-hilliard.py +++ b/examples/cahn-hilliard.py @@ -9,6 +9,7 @@ from meshmode.discretization.poly_element import \ from pytential import bind, sym, norm # noqa from pytential.target import PointsTarget + # {{{ set some constants for use below nelements = 20 @@ -23,10 +24,121 @@ vol_qbx_order = 2 # }}} +# {{{ a kernel class for G0 +# FIXME: will the expressions work when lambda is complex? +# (may need conversion from 1j to var("I")) +from sumpy.kernel import ExpressionKernel +class ShidongKernel(ExpressionKernel): + init_arg_names = ("dim", "lambda1", "lambda2") + + def __init__(self, dim=None, lambda1=0., lambda2=1.): + """ + :arg lambda1,lambda2: The roots of the quadratic equation w.r.t + laplacian. + """ + # Assert against repeated roots. + if abs(lambda1**2 - lambda2**2) < 1e-9: + raise RuntimeError("illposed since input roots are too close") + +# Based on http://mathworld.wolfram.com/ModifiedBesselFunctionoftheSecondKind.html + if dim == 2: + r = pymbolic_real_norm_2(make_sym_vector("d", dim)) + expr = var("hankel_1")(0, var("I") * lambda1 * r) - \ + var("hankel_1")(0, var("I") * lambda2 * r) + scaling = 1. / ( 4. * var("I") * (lambda1**2 - lambda2**2) ) + else: + raise RuntimeError("unsupported dimensionality") + + ExpressionKernel.__init__( + self, + dim, + expression=expr, + scaling = scaling, + is_complex_valued=True) + + self.lambda1 = lambda1 + self.lambda2 = lambda2 + + def __getinitargs__(self): + return(self._dim, self.lambda1, self.lambda2) + + def update_persistent_hash(self, key_hash, key_builder): + key_hash.update(type(self).__name__.encode("utf8")) + key_builder.rec(key_hash, + (self._dim, self.lambda1, self.lambda2) + ) + + def __repr__(self): + if self._dim is not None: + return "ShdgKnl%dD(%f, %f)" % ( + self._dim, self.lambda1, self.lambda2) + else: + return "ShdgKnl(%f, %f)" % (self.lambda1, self.lambda2) + + def prepare_loopy_kernel(self, loopy_knl): + from sumpy.codegen import (bessel_preamble_generator, bessel_mangler) + loopy_knl = lp.register_function_manglers(loopy_knl, + [bessel_mangler]) + loopy_knl = lp.register_preamble_generators(loopy_knl, + [bessel_preamble_generator]) + return loopy_knl + + def get_args(self): + k_dtype = np.complex128 + return [ + KernelArgument( + loopy_arg=lp.ValueArg("shidong_kernel", k_dtype), + )] + + mapper_method = "map_shidong_kernel" + +# }}} + +# {{{ extended kernel getters +def get_extkernel_for_G0(lambda1, lambda2): + from sumpy.symbolic import pymbolic_real_norm_2 + from pymbolic.primitives import make_sym_vector + from pymbolic import var + + d = make_sym_vector("d", 3) + r2 = pymbolic_real_norm_2(d[:-1]) + expr = var("hankel_1")(0, var("I") * lambda1 * r2 + + var("I") * d[-1]**2) \ + - var("hankel_1")(0, var("I") * lambda2 * r2 + + var("I") * d[-1]**2) + scaling = 1. / ( 4. * var("I") * (lambda1**2 - lambda2**2) ) + + from sumpy.kernel import ExpressionKernel + return ExpressionKernel( + dim=3, + expression=expr, + scaling=scaling, + is_complex_valued=True) + +def get_extkernel_for_G1(lamb): + from sumpy.symbolic import pymbolic_real_norm_2 + from pymbolic.primitives import make_sym_vector + from pymbolic import var + + d = make_sym_vector("d", 3) + r2 = pymbolic_real_norm_2(d[:-1]) + expr = var("hankel_1")(0, var("I") * lamb * r2 + + var("I") * d[-1]**2) + scaling = - var("I") / 4. + + from sumpy.kernel import ExpressionKernel + return ExpressionKernel( + dim=3, + expression=expr, + scaling=scaling, + is_complex_valued=True) + +# }}} def main(): import logging logging.basicConfig(level=logging.INFO) + logger = logging.getLogger(__name__) cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) @@ -57,11 +169,11 @@ def main(): 2, order=mesh_order, force_ambient_dim=2, - orther_options=["-string", + other_options=["-string", "Mesh.CharacteristicLengthMax = %g;" % h] ) else: - 1/0 + RuntimeError("unsupported geometry") logger.info("%d elements" % mesh.nelements) @@ -69,7 +181,7 @@ def main(): # {{{ discretization and connections - vol_discr = Discretization(ctx, mesh, + vol_discr = Discretization(cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(vol_quad_order)) from meshmode.mesh import BTAG_ALL @@ -118,117 +230,6 @@ def main(): # ).with_refinement() # density_discr = qbx.density_discr -# }}} - -# {{{ a kernel class for G0 -# FIXME: will the expressions work when lambda is complex? -# (may need conversion from 1j to var("I")) - class ShidongKernel(ExpressionKernel): - init_arg_names = ("dim", "lambda1", "lambda2") - - def __init__(self, dim=None, lambda1=0., lambda2=1.): - """ - :arg lambda1,lambda2: The roots of the quadratic equation w.r.t - laplacian. - """ - # Assert against repeated roots. - if abs(lambda1**2 - lambda2**2) < 1e-9: - raise RuntimeError("illposed since input roots are too close") - -# Based on http://mathworld.wolfram.com/ModifiedBesselFunctionoftheSecondKind.html - if dim == 2: - r = pymbolic_real_norm_2(make_sym_vector("d", dim)) - expr = var("hankel_1")(0, var("I") * lambda1 * r) - \ - var("hankel_1")(0, var("I") * lambda2 * r) - scaling = 1. / ( 4. * var("I") * (lambda1**2 - lambda2**2) ) - else: - raise RuntimeError("unsupported dimensionality") - - ExpressionKernel.__init__( - self, - dim, - expression=expr, - scaling = scaling, - is_complex_valued=True) - - self.lambda1 = lambda1 - self.lambda2 = lambda2 - - def __getinitargs__(self): - return(self._dim, self.lambda1, self.lambda2) - - def update_persistent_hash(self, key_hash, key_builder): - key_hash.update(type(self).__name__.encode("utf8")) - key_builder.rec(key_hash, - (self._dim, self.lambda1, self.lambda2) - ) - - def __repr__(self): - if self._dim is not None: - return "ShdgKnl%dD(%f, %f)" % ( - self._dim, self.lambda1, self.lambda2) - else: - return "ShdgKnl(%f, %f)" % (self.lambda1, self.lambda2) - - def prepare_loopy_kernel(self, loopy_knl): - from sumpy.codegen import (bessel_preamble_generator, bessel_mangler) - loopy_knl = lp.register_function_manglers(loopy_knl, - [bessel_mangler]) - loopy_knl = lp.register_preamble_generators(loopy_knl, - [bessel_preamble_generator]) - return loopy_knl - - def get_args(self): - k_dtype = np.complex128 - return [ - KernelArgument( - loopy_arg=lp.ValueArg("shidong_kernel", k_dtype), - )] - - mapper_method = "map_shidong_kernel" - -# }}} - -# {{{ extended kernel getters - def get_extkernel_for_G0(lambda1, lambda2): - from sumpy.symbolic import pymbolic_real_norm_2 - from pymbolic.primitives import make_sym_vector - from pymbolic import var - - d = make_sym_vector("d", 3) - r2 = pymbolic_real_norm_2(d[:-1]) - expr = var("hankel_1")(0, var("I") * lambda1 * r2 - + var("I") * d[-1]**2) \ - - var("hankel_1")(0, var("I") * lambda2 * r2 - + var("I") * d[-1]**2) - scaling = 1. / ( 4. * var("I") * (lambda1**2 - lambda2**2) ) - - from sumpy.kernel import ExpressionKernel - return ExpressionKernel( - dim=3, - expression=expr, - scaling=scaling, - is_complex_valued=True) - - def get_extkernel_for_G1(lamb): - from sumpy.symbolic import pymbolic_real_norm_2 - from pymbolic.primitives import make_sym_vector - from pymbolic import var - - d = make_sym_vector("d", 3) - r2 = pymbolic_real_norm_2(d[:-1]) - expr = var("hankel_1")(0, var("I") * lamb * r2 - + var("I") * d[-1]**2) - scaling = - var("I") / 4. - - from sumpy.kernel import ExpressionKernel - return ExpressionKernel( - dim=3, - expression=expr, - scaling=scaling, - is_complex_valued=True) - - # }}} # {{{ equation info @@ -239,6 +240,7 @@ def main(): b = s / (epsilon**2) c = 1. / (epsilon * delta_t) + print ( ["CH operator b = ", b, ", c = ", c] ) sqdet = np.sqrt( b**2 - 4. * c ) assert np.abs(sqdet) > 1e-6 @@ -247,7 +249,7 @@ def main(): lambda2 = ( b - sqdet ) / 2. from pytential.symbolic.pde.cahn_hilliard import CahnHilliardOperator - chop = CahnHilliardOperator(b=5, c=1) + chop = CahnHilliardOperator(b=b, c=c) unk = chop.make_unknown("sigma") bound_op = bind(qbx, chop.operator(unk)) @@ -257,10 +259,10 @@ def main(): from sumpy.qbx import LayerPotential from sumpy.qbx import LineTaylorLocalExpansion - layer_pot_v0f1 = LayerPotential(ctx, [ + layer_pot_v0f1 = LayerPotential(cl_ctx, [ LineTaylorLocalExpansion(shidong_2d_in_3d_kernel, order=vol_qbx_order)]) - layer_pot_v1f1 = LayerPotential(ctx, [ + layer_pot_v1f1 = LayerPotential(cl_ctx, [ LineTaylorLocalExpansion(yukawa_2d_in_3d_kernel_1, order=vol_qbx_order)]) # }}} @@ -323,7 +325,6 @@ def main(): # }}} - # {{{ fix rhs and solve nodes = density_discr.nodes().with_queue(queue) -- GitLab From ccdbc506c8e28e540da4ab0d4204562e42ec0528 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 2 Jun 2017 15:08:21 -0400 Subject: [PATCH 05/11] style fixes --- examples/cahn-hilliard.py | 85 ++++++++++++++++--------- pytential/symbolic/pde/cahn_hilliard.py | 15 +++-- 2 files changed, 64 insertions(+), 36 deletions(-) diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index 7c8e9ff8..6c531e1c 100644 --- a/examples/cahn-hilliard.py +++ b/examples/cahn-hilliard.py @@ -10,6 +10,12 @@ from meshmode.discretization.poly_element import \ from pytential import bind, sym, norm # noqa from pytential.target import PointsTarget +from pytools.obj_array import make_obj_array +import pytential.symbolic.primitives as p + +from sumpy.kernel import ExpressionKernel +import loopy as lp + # {{{ set some constants for use below nelements = 20 @@ -95,6 +101,17 @@ class ShidongKernel(ExpressionKernel): # }}} # {{{ extended kernel getters +class HankelBasedKernel(ExpressionKernel): + def prepare_loopy_kernel(self, loopy_knl): + from sumpy.codegen import (bessel_preamble_generator, bessel_mangler) + loopy_knl = lp.register_function_manglers(loopy_knl, + [bessel_mangler]) + loopy_knl = lp.register_preamble_generators(loopy_knl, + [bessel_preamble_generator]) + + return loopy_knl + + def get_extkernel_for_G0(lambda1, lambda2): from sumpy.symbolic import pymbolic_real_norm_2 from pymbolic.primitives import make_sym_vector @@ -108,8 +125,7 @@ def get_extkernel_for_G0(lambda1, lambda2): + var("I") * d[-1]**2) scaling = 1. / ( 4. * var("I") * (lambda1**2 - lambda2**2) ) - from sumpy.kernel import ExpressionKernel - return ExpressionKernel( + return HankelBasedKernel( dim=3, expression=expr, scaling=scaling, @@ -126,8 +142,7 @@ def get_extkernel_for_G1(lamb): + var("I") * d[-1]**2) scaling = - var("I") / 4. - from sumpy.kernel import ExpressionKernel - return ExpressionKernel( + return HankelBasedKernel( dim=3, expression=expr, scaling=scaling, @@ -240,25 +255,22 @@ def main(): b = s / (epsilon**2) c = 1. / (epsilon * delta_t) - print ( ["CH operator b = ", b, ", c = ", c] ) - - sqdet = np.sqrt( b**2 - 4. * c ) - assert np.abs(sqdet) > 1e-6 - - lambda1 = ( b + sqdet ) / 2. - lambda2 = ( b - sqdet ) / 2. + print("-- setup Cahn-Hilliard operator") from pytential.symbolic.pde.cahn_hilliard import CahnHilliardOperator - chop = CahnHilliardOperator(b=b, c=c) + chop = CahnHilliardOperator(b=b, c=c) unk = chop.make_unknown("sigma") bound_op = bind(qbx, chop.operator(unk)) - yukawa_2d_in_3d_kernel_1 = get_extkernel_for_G1(lambda1) - shidong_2d_in_3d_kernel = get_extkernel_for_G0(lambda1, lambda2) + print ("-- construct kernels") + yukawa_2d_in_3d_kernel_1 = get_extkernel_for_G1(chop.lambdas[0]) + shidong_2d_in_3d_kernel = get_extkernel_for_G0(chop.lambdas[0], + chop.lambdas[1]) + print("-- construct layer potentials") from sumpy.qbx import LayerPotential - from sumpy.qbx import LineTaylorLocalExpansion + from sumpy.expansion.local import LineTaylorLocalExpansion layer_pot_v0f1 = LayerPotential(cl_ctx, [ LineTaylorLocalExpansion(shidong_2d_in_3d_kernel, order=vol_qbx_order)]) @@ -269,7 +281,8 @@ def main(): # {{{ volume integral - vol_x = vol_discr.nodes(),with_queue(queue) + print("-- perform volume integrals") + vol_x = vol_discr.nodes().with_queue(queue) # target points targets = cl.array.zeros(queue, (3,) + vol_x.shape[1:], vol_x.dtype) @@ -354,26 +367,38 @@ def main(): def check_pde(): from sumpy.point_calculus import CalculusPatch - cp = CalculusPatch(np.zeros(2), order=4, h=0.1) - targets = cl.array.to_device(queue, cp.points) + vec_h = [1e-1, 1e-2, 1e-3, 1e-4] + vec_ru = [] + vec_rv = [] + for dx in vec_h: + cp = CalculusPatch(np.zeros(2), order=4, h=dx) + targets = cl.array.to_device(queue, cp.points) - u, v = bind( - (qbx, PointsTarget(targets)), - chop.representation(unk))(queue, sigma=sigma) + u, v = bind( + (qbx, PointsTarget(targets)), + chop.representation(unk))(queue, sigma=sigma) + + u = u.get().real + v = v.get().real + + lap_u = -(v - chop.b*u) - u = u.get().real - v = v.get().real + vec_ru.append(la.norm( + cp.laplace(lap_u) - chop.b * cp.laplace(u) + chop.c*u)) + vec_rv.append(la.norm( + v + cp.laplace(u) - chop.b*u)) - lap_u = -(v - chop.b*u) + # print(la.norm(u), la.norm(v)) - print(la.norm(u), la.norm(v)) + #print(la.norm( + # cp.laplace(lap_u) - chop.b * cp.laplace(u) + chop.c*u)) - print(la.norm( - cp.laplace(lap_u) - chop.b * cp.laplace(u) + chop.c*u)) + #print(la.norm( + # v + cp.laplace(u) - chop.b*u)) - print(la.norm( - v + cp.laplace(u) - chop.b*u)) - 1/0 + from tabulate import tabulate + print(tabulate([vec_h, vec_ru, vec_rv], + headers=['h', 'resid_u', 'resid_v'])) check_pde() diff --git a/pytential/symbolic/pde/cahn_hilliard.py b/pytential/symbolic/pde/cahn_hilliard.py index 3bbf1bf0..2ac6167c 100644 --- a/pytential/symbolic/pde/cahn_hilliard.py +++ b/pytential/symbolic/pde/cahn_hilliard.py @@ -40,19 +40,22 @@ class CahnHilliardOperator(L2WeightedPDEOperator): # Issue: # - let crat = np.abs(4.*c) / ( b**2 + 1e-12 ) - # - when crat is close to zero, sqrt(b**2-4*c) is close to abs(b), then for b>=0, - # sqrt(b**2-4*c) - b is inaccurate. - # - similarly, when crat is close to one, sqrt(b**2-4*c) is close to zero, then for b>0, - # sqrt(b**2-4*c) + b is inaccurate. + # - when crat is close to zero, sqrt(b**2-4*c) is close to abs(b), + # then for b>=0, sqrt(b**2-4*c) - b is inaccurate. + # - similarly, when crat is close to one, sqrt(b**2-4*c) is close to zero, + # then for b>0, sqrt(b**2-4*c) + b is inaccurate. # Solution: # - find a criteria for crat to choose from formulae, or # - use the computed root with smaller residual def quadratic_formula_1(a, b, c): - return ( -b + np.sqrt(b**2 - 4*a*c) ) / (2*a) + return (-b + np.sqrt(b**2 - 4*a*c)) / (2*a) + def quadratic_formula_2(a, b, c): - return ( -b - np.sqrt(b**2 - 4*a*c) ) / (2*a) + return (-b - np.sqrt(b**2 - 4*a*c)) / (2*a) + def citardauq_formula_1(a, b, c): return 2*c / (-b - np.sqrt(b**2-4*a*c)) + def citardauq_formula_2(a, b, c): return 2*c / (-b + np.sqrt(b**2-4*a*c)) -- GitLab From dc48b69231cfc8e29008aec70b11151d5f6fa66b Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 2 Jun 2017 15:59:49 -0400 Subject: [PATCH 06/11] pass check_pde --- examples/cahn-hilliard.py | 20 ++++++++++++++++++-- examples/check_pde.txt.cache | 5 +++++ 2 files changed, 23 insertions(+), 2 deletions(-) create mode 100644 examples/check_pde.txt.cache diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index 6c531e1c..6dfbcf71 100644 --- a/examples/cahn-hilliard.py +++ b/examples/cahn-hilliard.py @@ -304,6 +304,7 @@ def main(): sources[:2] = vol_x # a manufactured f1 + # FIXME: use f1 from the previous solution x_sin_factor = 30 y_sin_factor = 10 def f1_func(x, y): @@ -397,8 +398,13 @@ def main(): # v + cp.laplace(u) - chop.b*u)) from tabulate import tabulate - print(tabulate([vec_h, vec_ru, vec_rv], - headers=['h', 'resid_u', 'resid_v'])) + # overwrite if file exists + # using .cache simply because it is in .gitignore + with open('check_pde.txt.cache', 'w') as f: + print(tabulate([["h"] + vec_h, + ["residual_u"] + vec_ru, + ["residual_v"] +, vec_rv] + ), file=f) check_pde() @@ -422,6 +428,16 @@ def main(): sym.D(LaplaceKernel(2), sym.var("sigma")))( queue, sigma=ones_density).get() + # clean up the mess + def clean_file(filename): + import os + try: + os.remove(filename) + except OSError: + pass + clean_file("failed-targets.vts") + clean_file("potential.vts") + try: u, v = bind( (qbx_stick_out, PointsTarget(targets)), diff --git a/examples/check_pde.txt.cache b/examples/check_pde.txt.cache new file mode 100644 index 00000000..9ae95c22 --- /dev/null +++ b/examples/check_pde.txt.cache @@ -0,0 +1,5 @@ +---------- ----------- ----------- ----------- ----------- +h 0.1 0.01 0.001 0.0001 +residual_u 0.264407 0.000264403 4.94253e-07 3.97963e-05 +residual_v 1.76263e-05 1.76271e-08 3.29687e-11 2.55261e-09 +---------- ----------- ----------- ----------- ----------- -- GitLab From effc070b300f32aa5e1a7f93b27afb1c8f9b0c39 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 2 Jun 2017 16:00:49 -0400 Subject: [PATCH 07/11] deleted redundant file --- examples/check_pde.txt.cache | 5 ----- 1 file changed, 5 deletions(-) delete mode 100644 examples/check_pde.txt.cache diff --git a/examples/check_pde.txt.cache b/examples/check_pde.txt.cache deleted file mode 100644 index 9ae95c22..00000000 --- a/examples/check_pde.txt.cache +++ /dev/null @@ -1,5 +0,0 @@ ----------- ----------- ----------- ----------- ----------- -h 0.1 0.01 0.001 0.0001 -residual_u 0.264407 0.000264403 4.94253e-07 3.97963e-05 -residual_v 1.76263e-05 1.76271e-08 3.29687e-11 2.55261e-09 ----------- ----------- ----------- ----------- ----------- -- GitLab From 34e411d14beb9f85c941797e6c1861bf36d44001 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 2 Jun 2017 16:03:58 -0400 Subject: [PATCH 08/11] output ascii to .dat file --- .gitignore | 1 + examples/cahn-hilliard.py | 5 ++--- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.gitignore b/.gitignore index 3d5b8042..43a3d720 100644 --- a/.gitignore +++ b/.gitignore @@ -20,6 +20,7 @@ examples/*.pdf *.vtu *.vts +*.dat .cache diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index 6dfbcf71..f628f743 100644 --- a/examples/cahn-hilliard.py +++ b/examples/cahn-hilliard.py @@ -399,11 +399,10 @@ def main(): from tabulate import tabulate # overwrite if file exists - # using .cache simply because it is in .gitignore - with open('check_pde.txt.cache', 'w') as f: + with open('check_pde.dat', 'w') as f: print(tabulate([["h"] + vec_h, ["residual_u"] + vec_ru, - ["residual_v"] +, vec_rv] + ["residual_v"] + vec_rv] ), file=f) check_pde() -- GitLab From 1094a8e59b91c4eebe9d83f41d79310958916148 Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 3 Jun 2017 12:25:02 -0400 Subject: [PATCH 09/11] still struggling with dtype --- examples/cahn-hilliard.py | 213 +++++++++++++++++++++++++++++--------- 1 file changed, 164 insertions(+), 49 deletions(-) diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index f628f743..7373acf4 100644 --- a/examples/cahn-hilliard.py +++ b/examples/cahn-hilliard.py @@ -18,7 +18,7 @@ import loopy as lp # {{{ set some constants for use below -nelements = 20 +# {{{ all kinds of orders bdry_quad_order = 4 mesh_order = bdry_quad_order qbx_order = bdry_quad_order @@ -27,6 +27,57 @@ fmm_order = 8 vol_quad_order = 5 vol_qbx_order = 2 +# }}} +# {{{ mesh generation +nelements = 20 + +from enum import Enum +class Geometry(Enum): + RegularRectangle = 1 + Circle = 2 + +shape = Geometry.Circle +# }}} +# {{{ physical parameters +s = 1.5 +epsilon = 0.01 +delta_t = 0.05 +final_t = delta_t * 1 +theta_y = 60. / 180. * np.pi + +b = s / (epsilon**2) +c = 1. / (epsilon * delta_t) +# }}} +# {{{ initial phi + +# This phi function is also used to do PDE check +import pymbolic as pmbl +x = pmbl.var("x") +y = pmbl.var("y") + +# FIXME: modify pymbolic to use tanh function +# phi = tanh(x / sqrt (2 * epsilon)) +phi = x**2 + x*y + 2*y**2 +phi3 = phi**3 +laphi = pmbl.differentiate(pmbl.differentiate(phi, 'x'), 'x') + \ + pmbl.differentiate(pmbl.differentiate(phi, 'y'), 'y') +laphi3 = pmbl.differentiate(pmbl.differentiate(phi3, 'x'), 'x') + \ + pmbl.differentiate(pmbl.differentiate(phi3, 'y'), 'y') +f1_expr = c * phi - (1+s) / epsilon**2 * laphi + 1 / epsilon**2 * laphi3 +f2_expr = (phi3 - (1+s) * phi) / epsilon**2 + +def f1_func(x,y): + return pmbl.evaluate(f1_expr, {"x": x, "y": y}) + +def f2_func(x,y): + return pmbl.evaluate(f2_expr, {"x": x, "y": y}) + +def initial_phi(x,y): + return pmbl.evaluate(phi, {"x": x, "y": y}) + +#def initial_phi(x, y): +# return np.tanh(x / np.sqrt(2. * initial_epsilon)) +# }}} # }}} @@ -160,13 +211,6 @@ def main(): # {{{ volume mesh generation - from enum import Enum - class Geometry(Enum): - RegularRectangle = 1 - Circle = 2 - - shape = Geometry.Circle - if shape == Geometry.RegularRectangle: from meshmode.mesh.generation import generate_regular_rect_mesh ext = 1. @@ -247,28 +291,21 @@ def main(): # }}} -# {{{ equation info +# {{{ setup operator and potentials - s = 1.5 - epsilon = 0.01 - delta_t = 0.05 - - b = s / (epsilon**2) - c = 1. / (epsilon * delta_t) - - print("-- setup Cahn-Hilliard operator") + #print("-- setup Cahn-Hilliard operator") from pytential.symbolic.pde.cahn_hilliard import CahnHilliardOperator chop = CahnHilliardOperator(b=b, c=c) unk = chop.make_unknown("sigma") bound_op = bind(qbx, chop.operator(unk)) - print ("-- construct kernels") + #print ("-- construct kernels") yukawa_2d_in_3d_kernel_1 = get_extkernel_for_G1(chop.lambdas[0]) shidong_2d_in_3d_kernel = get_extkernel_for_G0(chop.lambdas[0], chop.lambdas[1]) - print("-- construct layer potentials") + #print("-- construct layer potentials") from sumpy.qbx import LayerPotential from sumpy.expansion.local import LineTaylorLocalExpansion layer_pot_v0f1 = LayerPotential(cl_ctx, [ @@ -279,9 +316,8 @@ def main(): order=vol_qbx_order)]) # }}} -# {{{ volume integral +# {{{ setup for volume integral - print("-- perform volume integrals") vol_x = vol_discr.nodes().with_queue(queue) # target points @@ -299,19 +335,10 @@ def main(): print(center_dist) # source points - # FIXME: use over sampled source points? + # TODO: use over sampled source points? sources = cl.array.zeros(queue, (3,) + vol_x.shape[1:], vol_x.dtype) sources[:2] = vol_x - # a manufactured f1 - # FIXME: use f1 from the previous solution - x_sin_factor = 30 - y_sin_factor = 10 - def f1_func(x, y): - return 0.1 * cl.clmath.sin(x_sin_factor*x) * cl.clmath.sin(y_sin_factor*y) - - # strengths (with quadrature weights) - f1 = f1_func(vol_x[0], vol_x[1]) vol_weights = bind(vol_discr, p.area_element(mesh.ambient_dim, mesh.dim) * p.QWeight() )(queue) @@ -319,26 +346,53 @@ def main(): print("volume: %d source nodes, %d target nodes" % ( vol_discr.nnodes, vol_discr.nnodes)) - evt, (vol_pot_v0f1,) = layer_pot_v0f1( - queue, - targets=targets.reshape(3, vol_discr.nnodes), - centers=centers, - sources=sources.reshape(3, vol_discr.nnodes), - strengths=( - (vol_weights * f1).reshape(vol_discr.nnodes),) - ) - evt, (vol_pot_v1f1,) = layer_pot_v1f1( - queue, - targets=targets.reshape(3, vol_discr.nnodes), - centers=centers, - sources=sources.reshape(3, vol_discr.nnodes), - strengths=( - (vol_weights * f1).reshape(vol_discr.nnodes),) - ) +# }}} +# {{{ prepare for time stepping + timestep_number = 0 + time = 0 + def get_vts_filename(tmstp_num): + return "solution-" + '{0:03}'.format(tmstp_num) + ".vts" + output_vts_filename = get_vts_filename(timestep_number) # }}} +# {{{ [[TIME STEPPING]] + while time < final_t: + timestep_number += 1 + time += delta_t + output_vts_filename = get_vts_filename(timestep_number) + + # a manufactured f1 function + #x_sin_factor = 30 + #y_sin_factor = 10 + #def f1_func(x, y): + # return 0.1 * cl.clmath.sin(x_sin_factor*x) \ + # * cl.clmath.sin(y_sin_factor*y) + + # get f1 to compute strengths + f1 = f1_func(vol_x[0], vol_x[1]) + f2 = f2_func(vol_x[0], vol_x[1]) + + evt, (vol_pot_v0f1,) = layer_pot_v0f1( + queue, + targets=targets.reshape(3, vol_discr.nnodes), + centers=centers, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=( + (vol_weights * f1).reshape(vol_discr.nnodes),) + ) + + evt, (vol_pot_v1f1,) = layer_pot_v1f1( + queue, + targets=targets.reshape(3, vol_discr.nnodes), + centers=centers, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=( + (vol_weights * f1).reshape(vol_discr.nnodes),) + ) + + # {{{ fix rhs and solve nodes = density_discr.nodes().with_queue(queue) @@ -371,6 +425,8 @@ def main(): vec_h = [1e-1, 1e-2, 1e-3, 1e-4] vec_ru = [] vec_rv = [] + vec_rp = [] + vec_rm = [] for dx in vec_h: cp = CalculusPatch(np.zeros(2), order=4, h=dx) targets = cl.array.to_device(queue, cp.points) @@ -384,11 +440,65 @@ def main(): lap_u = -(v - chop.b*u) + # Check for homogeneous PDEs for u and v vec_ru.append(la.norm( cp.laplace(lap_u) - chop.b * cp.laplace(u) + chop.c*u)) vec_rv.append(la.norm( v + cp.laplace(u) - chop.b*u)) + # Check for inhomogeneous PDEs for phi and mu + + targets_in_3d = cl.array.zeros( + queue, + (3,) + cp.points.shape[1:], + vol_x.dtype) + targets_in_3d[:2] = cp.points + + #center_dist = 0.125*np.min( + # cl.clmath.sqrt( + # bind(vol_discr, + # p.area_element(mesh.ambient_dim, mesh.dim)) + # (queue)).get()) + + centers_in_3d = make_obj_array( + [ci.copy() + for ci in targets_in_3d]) + centers_in_3d[2][:] = center_dist + + evt, (v0f1,) = layer_pot_v0f1( + queue, + targets=targets_in_3d, + centers=centers_in_3d, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=( + (vol_weights * f1).reshape(vol_discr.nnodes),) + ) + + evt, (v1f1,) = layer_pot_v1f1( + queue, + targets=targets_in_3d, + centers=centers_in_3d, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=( + (vol_weights * f1).reshape(vol_discr.nnodes),) + ) + + f14ck = f1_func(cp.points[0], cp.points[1]) + f24ck = f2_func(cp.points[0], cp.points[1]) + utild = v0f1 + vtild = f24ck - v1f1 + lambda1**2 * v0f1 + + ph = utild + u + mu = epsilon * (vtild + v) + lap_ph = f24ck + chop.b * ph - mu / epsilon + + vec_rp.append(la.norm( + cp.laplace(lap_ph) - chop.b * cp.laplace(ph) + chop.c*ph + - f14ck)) + vec_rm.append(la.norm( + mu / epsilon + cp.laplace(ph) - chop.b*ph + - f24ck)) + # print(la.norm(u), la.norm(v)) #print(la.norm( @@ -402,7 +512,10 @@ def main(): with open('check_pde.dat', 'w') as f: print(tabulate([["h"] + vec_h, ["residual_u"] + vec_ru, - ["residual_v"] + vec_rv] + ["residual_v"] + vec_rv, + ["residual_phi"] + vec_rp, + ["residual_mu"] + vec_rm, + ] ), file=f) check_pde() @@ -412,7 +525,7 @@ def main(): # {{{ postprocess/visualize from sumpy.visualization import FieldPlotter - fplot = FieldPlotter(np.zeros(2), extent=5, npoints=500) + fplot = FieldPlotter(np.zeros(2), extent=1.5, npoints=500) targets = cl.array.to_device(queue, fplot.points) @@ -464,6 +577,8 @@ def main(): # }}} +# }}} + if __name__ == "__main__": main() -- GitLab From adaa8e65888ff264f1de55df48440eb5ad33f37b Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 4 Jun 2017 16:23:35 -0400 Subject: [PATCH 10/11] debugging volume integral --- examples/cahn-hilliard.py | 487 ++++++++++++++++++++------------------ 1 file changed, 256 insertions(+), 231 deletions(-) diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index 7373acf4..4cbd1680 100644 --- a/examples/cahn-hilliard.py +++ b/examples/cahn-hilliard.py @@ -22,19 +22,22 @@ import loopy as lp bdry_quad_order = 4 mesh_order = bdry_quad_order qbx_order = bdry_quad_order -bdry_ovsmp_quad_order = 4*bdry_quad_order +bdry_ovsmp_quad_order = 4 * bdry_quad_order fmm_order = 8 vol_quad_order = 5 -vol_qbx_order = 2 +vol_qbx_order = 2 # }}} # {{{ mesh generation nelements = 20 from enum import Enum + + class Geometry(Enum): RegularRectangle = 1 - Circle = 2 + Circle = 2 + shape = Geometry.Circle # }}} @@ -57,24 +60,28 @@ y = pmbl.var("y") # FIXME: modify pymbolic to use tanh function # phi = tanh(x / sqrt (2 * epsilon)) -phi = x**2 + x*y + 2*y**2 +phi = x**2 + x * y + 2 * y**2 phi3 = phi**3 laphi = pmbl.differentiate(pmbl.differentiate(phi, 'x'), 'x') + \ pmbl.differentiate(pmbl.differentiate(phi, 'y'), 'y') laphi3 = pmbl.differentiate(pmbl.differentiate(phi3, 'x'), 'x') + \ pmbl.differentiate(pmbl.differentiate(phi3, 'y'), 'y') -f1_expr = c * phi - (1+s) / epsilon**2 * laphi + 1 / epsilon**2 * laphi3 -f2_expr = (phi3 - (1+s) * phi) / epsilon**2 +f1_expr = c * phi - (1 + s) / epsilon**2 * laphi + 1 / epsilon**2 * laphi3 +f2_expr = (phi3 - (1 + s) * phi) / epsilon**2 + -def f1_func(x,y): +def f1_func(x, y): return pmbl.evaluate(f1_expr, {"x": x, "y": y}) -def f2_func(x,y): + +def f2_func(x, y): return pmbl.evaluate(f2_expr, {"x": x, "y": y}) -def initial_phi(x,y): + +def initial_phi(x, y): return pmbl.evaluate(phi, {"x": x, "y": y}) + #def initial_phi(x, y): # return np.tanh(x / np.sqrt(2. * initial_epsilon)) # }}} @@ -85,6 +92,8 @@ def initial_phi(x,y): # FIXME: will the expressions work when lambda is complex? # (may need conversion from 1j to var("I")) from sumpy.kernel import ExpressionKernel + + class ShidongKernel(ExpressionKernel): init_arg_names = ("dim", "lambda1", "lambda2") @@ -99,66 +108,64 @@ class ShidongKernel(ExpressionKernel): # Based on http://mathworld.wolfram.com/ModifiedBesselFunctionoftheSecondKind.html if dim == 2: - r = pymbolic_real_norm_2(make_sym_vector("d", dim)) + r = pymbolic_real_norm_2(make_sym_vector("d", dim)) expr = var("hankel_1")(0, var("I") * lambda1 * r) - \ var("hankel_1")(0, var("I") * lambda2 * r) - scaling = 1. / ( 4. * var("I") * (lambda1**2 - lambda2**2) ) + scaling = 1. / (4. * var("I") * (lambda1**2 - lambda2**2)) else: raise RuntimeError("unsupported dimensionality") ExpressionKernel.__init__( - self, - dim, - expression=expr, - scaling = scaling, - is_complex_valued=True) + self, + dim, + expression=expr, + scaling=scaling, + is_complex_valued=True) self.lambda1 = lambda1 self.lambda2 = lambda2 def __getinitargs__(self): - return(self._dim, self.lambda1, self.lambda2) + return (self._dim, self.lambda1, self.lambda2) def update_persistent_hash(self, key_hash, key_builder): key_hash.update(type(self).__name__.encode("utf8")) - key_builder.rec(key_hash, - (self._dim, self.lambda1, self.lambda2) - ) + key_builder.rec(key_hash, (self._dim, self.lambda1, self.lambda2)) def __repr__(self): if self._dim is not None: - return "ShdgKnl%dD(%f, %f)" % ( - self._dim, self.lambda1, self.lambda2) + return "ShdgKnl%dD(%f, %f)" % (self._dim, self.lambda1, + self.lambda2) else: return "ShdgKnl(%f, %f)" % (self.lambda1, self.lambda2) def prepare_loopy_kernel(self, loopy_knl): from sumpy.codegen import (bessel_preamble_generator, bessel_mangler) - loopy_knl = lp.register_function_manglers(loopy_knl, - [bessel_mangler]) - loopy_knl = lp.register_preamble_generators(loopy_knl, - [bessel_preamble_generator]) + loopy_knl = lp.register_function_manglers(loopy_knl, [bessel_mangler]) + loopy_knl = lp.register_preamble_generators( + loopy_knl, [bessel_preamble_generator]) return loopy_knl def get_args(self): k_dtype = np.complex128 return [ - KernelArgument( - loopy_arg=lp.ValueArg("shidong_kernel", k_dtype), - )] + KernelArgument( + loopy_arg=lp.ValueArg("shidong_kernel", k_dtype), ) + ] mapper_method = "map_shidong_kernel" + # }}} + # {{{ extended kernel getters class HankelBasedKernel(ExpressionKernel): def prepare_loopy_kernel(self, loopy_knl): from sumpy.codegen import (bessel_preamble_generator, bessel_mangler) - loopy_knl = lp.register_function_manglers(loopy_knl, - [bessel_mangler]) - loopy_knl = lp.register_preamble_generators(loopy_knl, - [bessel_preamble_generator]) + loopy_knl = lp.register_function_manglers(loopy_knl, [bessel_mangler]) + loopy_knl = lp.register_preamble_generators( + loopy_knl, [bessel_preamble_generator]) return loopy_knl @@ -174,13 +181,11 @@ def get_extkernel_for_G0(lambda1, lambda2): + var("I") * d[-1]**2) \ - var("hankel_1")(0, var("I") * lambda2 * r2 + var("I") * d[-1]**2) - scaling = 1. / ( 4. * var("I") * (lambda1**2 - lambda2**2) ) + scaling = 1. / (4. * var("I") * (lambda1**2 - lambda2**2)) return HankelBasedKernel( - dim=3, - expression=expr, - scaling=scaling, - is_complex_valued=True) + dim=3, expression=expr, scaling=scaling, is_complex_valued=True) + def get_extkernel_for_G1(lamb): from sumpy.symbolic import pymbolic_real_norm_2 @@ -189,18 +194,16 @@ def get_extkernel_for_G1(lamb): d = make_sym_vector("d", 3) r2 = pymbolic_real_norm_2(d[:-1]) - expr = var("hankel_1")(0, var("I") * lamb * r2 - + var("I") * d[-1]**2) - scaling = - var("I") / 4. + expr = var("hankel_1")(0, var("I") * lamb * r2 + var("I") * d[-1]**2) + scaling = -var("I") / 4. return HankelBasedKernel( - dim=3, - expression=expr, - scaling=scaling, - is_complex_valued=True) + dim=3, expression=expr, scaling=scaling, is_complex_valued=True) + # }}} + def main(): import logging logging.basicConfig(level=logging.INFO) @@ -209,57 +212,56 @@ def main(): cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) -# {{{ volume mesh generation + # {{{ volume mesh generation if shape == Geometry.RegularRectangle: - from meshmode.mesh.generation import generate_regular_rect_mesh - ext = 1. - h = 0.05 - mesh = generate_regular_rect_mesh( - a=(-ext/2., -ext/2.), - b=( ext/2., ext/2.), - n=(int(ext/h), int(ext/h)) - ) + from meshmode.mesh.generation import generate_regular_rect_mesh + ext = 1. + h = 0.05 + mesh = generate_regular_rect_mesh( + a=(-ext / 2., -ext / 2.), + b=(ext / 2., ext / 2.), + n=(int(ext / h), int(ext / h))) elif shape == Geometry.Circle: from meshmode.mesh.io import generate_gmsh, FileSource - h = 0.05 + h = 0.05 mesh = generate_gmsh( - FileSource("circle.step"), - 2, - order=mesh_order, - force_ambient_dim=2, - other_options=["-string", - "Mesh.CharacteristicLengthMax = %g;" % h] - ) + FileSource("circle.step"), + 2, + order=mesh_order, + force_ambient_dim=2, + other_options=[ + "-string", "Mesh.CharacteristicLengthMax = %g;" % h + ]) else: RuntimeError("unsupported geometry") logger.info("%d elements" % mesh.nelements) -# }}} + # }}} -# {{{ discretization and connections + # {{{ discretization and connections - vol_discr = Discretization(cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(vol_quad_order)) + vol_discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(vol_quad_order)) from meshmode.mesh import BTAG_ALL from meshmode.discretization.connection import make_face_restriction pre_density_connection = make_face_restriction( - vol_discr, - InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order), - BTAG_ALL) + vol_discr, + InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order), BTAG_ALL) pre_density_discr = pre_density_connection.to_discr - from pytential.qbx import ( - QBXLayerPotentialSource, QBXTargetAssociationFailedException) + from pytential.qbx import (QBXLayerPotentialSource, + QBXTargetAssociationFailedException) qbx, _ = QBXLayerPotentialSource( - pre_density_discr, - fine_order=bdry_ovsmp_quad_order, qbx_order=qbx_order, - fmm_order=fmm_order, - expansion_disks_in_tree_have_extent=True, - ).with_refinement() + pre_density_discr, + fine_order=bdry_ovsmp_quad_order, + qbx_order=qbx_order, + fmm_order=fmm_order, + expansion_disks_in_tree_have_extent=True, ).with_refinement() density_discr = qbx.density_discr @@ -267,97 +269,98 @@ def main(): # vol_discr --> pre_density_discr --> density_discr # via ChainedDiscretizationConnection + # from meshmode.mesh.generation import ellipse, make_curve_mesh + # from functools import partial + # + # mesh = make_curve_mesh( + # partial(ellipse, 2), + # np.linspace(0, 1, nelements+1), + # mesh_order) + # + # pre_density_discr = Discretization( + # cl_ctx, mesh, + # InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order)) + # + # from pytential.qbx import ( + # QBXLayerPotentialSource, QBXTargetAssociationFailedException) + # qbx, _ = QBXLayerPotentialSource( + # pre_density_discr, fine_order=bdry_ovsmp_quad_order, qbx_order=qbx_order, + # fmm_order=fmm_order, + # expansion_disks_in_tree_have_extent=True, + # ).with_refinement() + # density_discr = qbx.density_discr -# from meshmode.mesh.generation import ellipse, make_curve_mesh -# from functools import partial -# -# mesh = make_curve_mesh( -# partial(ellipse, 2), -# np.linspace(0, 1, nelements+1), -# mesh_order) -# -# pre_density_discr = Discretization( -# cl_ctx, mesh, -# InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order)) -# -# from pytential.qbx import ( -# QBXLayerPotentialSource, QBXTargetAssociationFailedException) -# qbx, _ = QBXLayerPotentialSource( -# pre_density_discr, fine_order=bdry_ovsmp_quad_order, qbx_order=qbx_order, -# fmm_order=fmm_order, -# expansion_disks_in_tree_have_extent=True, -# ).with_refinement() -# density_discr = qbx.density_discr - -# }}} + # }}} -# {{{ setup operator and potentials + # {{{ setup operator and potentials #print("-- setup Cahn-Hilliard operator") from pytential.symbolic.pde.cahn_hilliard import CahnHilliardOperator - chop = CahnHilliardOperator(b=b, c=c) + chop = CahnHilliardOperator(b=b, c=c) unk = chop.make_unknown("sigma") bound_op = bind(qbx, chop.operator(unk)) #print ("-- construct kernels") yukawa_2d_in_3d_kernel_1 = get_extkernel_for_G1(chop.lambdas[0]) - shidong_2d_in_3d_kernel = get_extkernel_for_G0(chop.lambdas[0], - chop.lambdas[1]) + shidong_2d_in_3d_kernel = get_extkernel_for_G0(chop.lambdas[0], + chop.lambdas[1]) #print("-- construct layer potentials") from sumpy.qbx import LayerPotential from sumpy.expansion.local import LineTaylorLocalExpansion layer_pot_v0f1 = LayerPotential(cl_ctx, [ - LineTaylorLocalExpansion(shidong_2d_in_3d_kernel, - order=vol_qbx_order)]) + LineTaylorLocalExpansion(shidong_2d_in_3d_kernel, order=vol_qbx_order) + ]) layer_pot_v1f1 = LayerPotential(cl_ctx, [ - LineTaylorLocalExpansion(yukawa_2d_in_3d_kernel_1, - order=vol_qbx_order)]) -# }}} + LineTaylorLocalExpansion( + yukawa_2d_in_3d_kernel_1, order=vol_qbx_order) + ]) + # }}} -# {{{ setup for volume integral + # {{{ setup for volume integral vol_x = vol_discr.nodes().with_queue(queue) # target points - targets = cl.array.zeros(queue, (3,) + vol_x.shape[1:], vol_x.dtype) + targets = cl.array.zeros(queue, (3, ) + vol_x.shape[1:], vol_x.dtype) targets[:2] = vol_x # expansion centers - center_dist = 0.125*np.min( - cl.clmath.sqrt( - bind(vol_discr, - p.area_element(mesh.ambient_dim, mesh.dim)) - (queue)).get()) - centers = make_obj_array([ci.copy().reshape(vol_discr.nnodes) for ci in targets]) + center_dist = 0.125 * np.min( + cl.clmath.sqrt( + bind(vol_discr, p.area_element(mesh.ambient_dim, mesh.dim))(queue)) + .get()) + centers = make_obj_array( + [ci.copy().reshape(vol_discr.nnodes) for ci in targets]) centers[2][:] = center_dist print(center_dist) # source points # TODO: use over sampled source points? - sources = cl.array.zeros(queue, (3,) + vol_x.shape[1:], vol_x.dtype) + sources = cl.array.zeros(queue, (3, ) + vol_x.shape[1:], vol_x.dtype) sources[:2] = vol_x - vol_weights = bind(vol_discr, - p.area_element(mesh.ambient_dim, mesh.dim) * p.QWeight() - )(queue) + vol_weights = bind( + vol_discr, + p.area_element(mesh.ambient_dim, mesh.dim) * p.QWeight())(queue) - print("volume: %d source nodes, %d target nodes" % ( - vol_discr.nnodes, vol_discr.nnodes)) + print("volume: %d source nodes, %d target nodes" % (vol_discr.nnodes, + vol_discr.nnodes)) + # }}} -# }}} - -# {{{ prepare for time stepping + # {{{ prepare for time stepping timestep_number = 0 - time = 0 + time = 0 + def get_vts_filename(tmstp_num): return "solution-" + '{0:03}'.format(tmstp_num) + ".vts" + output_vts_filename = get_vts_filename(timestep_number) -# }}} + # }}} -# {{{ [[TIME STEPPING]] + # {{{ [[TIME STEPPING]] while time < final_t: timestep_number += 1 time += delta_t @@ -374,24 +377,19 @@ def main(): f1 = f1_func(vol_x[0], vol_x[1]) f2 = f2_func(vol_x[0], vol_x[1]) - evt, (vol_pot_v0f1,) = layer_pot_v0f1( - queue, - targets=targets.reshape(3, vol_discr.nnodes), - centers=centers, - sources=sources.reshape(3, vol_discr.nnodes), - strengths=( - (vol_weights * f1).reshape(vol_discr.nnodes),) - ) - - evt, (vol_pot_v1f1,) = layer_pot_v1f1( - queue, - targets=targets.reshape(3, vol_discr.nnodes), - centers=centers, - sources=sources.reshape(3, vol_discr.nnodes), - strengths=( - (vol_weights * f1).reshape(vol_discr.nnodes),) - ) + evt, (vol_pot_v0f1, ) = layer_pot_v0f1( + queue, + targets=targets.reshape(3, vol_discr.nnodes), + centers=centers, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=((vol_weights * f1).reshape(vol_discr.nnodes), )) + evt, (vol_pot_v1f1, ) = layer_pot_v1f1( + queue, + targets=targets.reshape(3, vol_discr.nnodes), + centers=centers, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=((vol_weights * f1).reshape(vol_discr.nnodes), )) # {{{ fix rhs and solve @@ -399,20 +397,22 @@ def main(): def g(xvec): x, y = xvec - return cl.clmath.cos(5*cl.clmath.atan2(y, x)) + return cl.clmath.cos(5 * cl.clmath.atan2(y, x)) bc = sym.make_obj_array([ # FIXME: Realistic BC - 5+g(nodes), - 3-g(nodes), - ]) + 5 + g(nodes), + 3 - g(nodes), + ]) from pytential.solve import gmres gmres_result = gmres( - bound_op.scipy_op(queue, "sigma", dtype=np.complex128), - bc, tol=1e-8, progress=True, - stall_iterations=0, - hard_failure=True) + bound_op.scipy_op(queue, "sigma", dtype=np.complex128), + bc, + tol=1e-8, + progress=True, + stall_iterations=0, + hard_failure=True) sigma = gmres_result.solution @@ -422,36 +422,38 @@ def main(): def check_pde(): from sumpy.point_calculus import CalculusPatch - vec_h = [1e-1, 1e-2, 1e-3, 1e-4] + vec_h = [1e-1, 1e-2, 1e-3, 1e-4] vec_ru = [] vec_rv = [] vec_rp = [] vec_rm = [] + vec_rf1 = [] + vec_rf2 = [] + vec_rutld = [] + vec_rvtld = [] for dx in vec_h: cp = CalculusPatch(np.zeros(2), order=4, h=dx) targets = cl.array.to_device(queue, cp.points) - u, v = bind( - (qbx, PointsTarget(targets)), - chop.representation(unk))(queue, sigma=sigma) + u, v = bind((qbx, PointsTarget(targets)), + chop.representation(unk))( + queue, sigma=sigma) u = u.get().real v = v.get().real - lap_u = -(v - chop.b*u) + lap_u = -(v - chop.b * u) # Check for homogeneous PDEs for u and v - vec_ru.append(la.norm( - cp.laplace(lap_u) - chop.b * cp.laplace(u) + chop.c*u)) - vec_rv.append(la.norm( - v + cp.laplace(u) - chop.b*u)) + vec_ru.append( + la.norm( + cp.laplace(lap_u) - chop.b * cp.laplace(u) + chop.c * u)) + vec_rv.append(la.norm(v + cp.laplace(u) - chop.b * u)) # Check for inhomogeneous PDEs for phi and mu - targets_in_3d = cl.array.zeros( - queue, - (3,) + cp.points.shape[1:], - vol_x.dtype) + targets_in_3d = cl.array.zeros(queue, (3, ) + cp.points.shape[1:], + vol_x.dtype) targets_in_3d[:2] = cp.points #center_dist = 0.125*np.min( @@ -460,63 +462,93 @@ def main(): # p.area_element(mesh.ambient_dim, mesh.dim)) # (queue)).get()) - centers_in_3d = make_obj_array( - [ci.copy() - for ci in targets_in_3d]) + centers_in_3d = make_obj_array([ci.copy() for ci in targets_in_3d]) centers_in_3d[2][:] = center_dist - evt, (v0f1,) = layer_pot_v0f1( - queue, - targets=targets_in_3d, - centers=centers_in_3d, - sources=sources.reshape(3, vol_discr.nnodes), - strengths=( - (vol_weights * f1).reshape(vol_discr.nnodes),) - ) - - evt, (v1f1,) = layer_pot_v1f1( - queue, - targets=targets_in_3d, - centers=centers_in_3d, - sources=sources.reshape(3, vol_discr.nnodes), - strengths=( - (vol_weights * f1).reshape(vol_discr.nnodes),) - ) - - f14ck = f1_func(cp.points[0], cp.points[1]) - f24ck = f2_func(cp.points[0], cp.points[1]) + evt, (v0f1, ) = layer_pot_v0f1( + queue, + targets=targets_in_3d, + centers=centers_in_3d, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=((vol_weights * f1).reshape(vol_discr.nnodes), )) + v0f1 = v0f1.get().real + + evt, (v1f1, ) = layer_pot_v1f1( + queue, + targets=targets_in_3d, + centers=centers_in_3d, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=((vol_weights * f1).reshape(vol_discr.nnodes), )) + v1f1 = v1f1.get().real + + f14ck = cl.array.to_device(queue, + f1_func(cp.points[0], cp.points[1])) + f24ck = cl.array.to_device(queue, + f2_func(cp.points[0], cp.points[1])) + phi_init = cl.array.to_device(queue, + initial_phi(cp.points[0], + cp.points[1])) + + f14ck = f14ck.get().real + f24ck = f24ck.get().real + phi_init = phi_init.get().real + vec_rf1.append( + la.norm(chop.c * phi_init - + (1 + s) / epsilon**2 * cp.laplace(phi_init) + + 1. / epsilon**2 * cp.laplace(phi_init**3) - f14ck)) + vec_rf2.append( + la.norm((phi_init**3 - (1 + s) * phi_init) / epsilon**2 - + f24ck)) + utild = v0f1 - vtild = f24ck - v1f1 + lambda1**2 * v0f1 + vtild = f24ck - v1f1 + chop.lambdas[0]**2 * v0f1 + lap_utild = -vtild + chop.b * utild + f24ck + # FIXME: Not passing this check. (bugs in volume integral?) + vec_rutld.append( + la.norm( + cp.laplace(lap_utild) - chop.b * lap_utild + + chop.c * utild - f14ck)) + print(f14ck) + print(f1[0:10]) + print(cp.laplace(lap_utild)) + print((lap_utild)) + print((utild)) + print(chop.b) + print(chop.c) + vec_rvtld.append( + la.norm(vtild + cp.laplace(utild) - chop.b * utild - f24ck)) ph = utild + u mu = epsilon * (vtild + v) lap_ph = f24ck + chop.b * ph - mu / epsilon - vec_rp.append(la.norm( - cp.laplace(lap_ph) - chop.b * cp.laplace(ph) + chop.c*ph - - f14ck)) - vec_rm.append(la.norm( - mu / epsilon + cp.laplace(ph) - chop.b*ph - - f24ck)) - - # print(la.norm(u), la.norm(v)) - - #print(la.norm( - # cp.laplace(lap_u) - chop.b * cp.laplace(u) + chop.c*u)) - - #print(la.norm( - # v + cp.laplace(u) - chop.b*u)) + vec_rp.append( + la.norm( + cp.laplace(lap_ph) - chop.b * cp.laplace(ph) + chop.c * ph + - f14ck)) + vec_rm.append( + la.norm(mu / epsilon + cp.laplace(ph) - chop.b * ph - f24ck)) from tabulate import tabulate # overwrite if file exists with open('check_pde.dat', 'w') as f: - print(tabulate([["h"] + vec_h, - ["residual_u"] + vec_ru, - ["residual_v"] + vec_rv, - ["residual_phi"] + vec_rp, - ["residual_mu"] + vec_rm, - ] - ), file=f) + print( + "Residuals of PDE and numerical vs. symbolic differentiation:", + file=f) + print( + tabulate([ + ["h"] + vec_h, + ["residual_u"] + vec_ru, + ["residual_v"] + vec_rv, + ["residual_f1"] + vec_rf1, + ["residual_f2"] + vec_rf2, + ["residual_u_tld"] + vec_rutld, + ["residual_v_tld"] + vec_rvtld, + ["residual_phi"] + vec_rp, + ["residual_mu"] + vec_rm, + ]), + file=f) + 1 / 0 check_pde() @@ -535,10 +567,9 @@ def main(): from sumpy.kernel import LaplaceKernel ones_density = density_discr.zeros(queue) ones_density.fill(1) - indicator = bind( - (indicator_qbx, PointsTarget(targets)), - sym.D(LaplaceKernel(2), sym.var("sigma")))( - queue, sigma=ones_density).get() + indicator = bind((indicator_qbx, PointsTarget(targets)), + sym.D(LaplaceKernel(2), sym.var("sigma")))( + queue, sigma=ones_density).get() # clean up the mess def clean_file(filename): @@ -547,33 +578,27 @@ def main(): os.remove(filename) except OSError: pass + clean_file("failed-targets.vts") clean_file("potential.vts") try: - u, v = bind( - (qbx_stick_out, PointsTarget(targets)), - chop.representation(unk))(queue, sigma=sigma) + u, v = bind((qbx_stick_out, PointsTarget(targets)), + chop.representation(unk))( + queue, sigma=sigma) except QBXTargetAssociationFailedException as e: - fplot.write_vtk_file( - "failed-targets.vts", - [ - ("failed", e.failed_target_flags.get(queue)) - ] - ) + fplot.write_vtk_file("failed-targets.vts", + [("failed", e.failed_target_flags.get(queue))]) raise u = u.get().real v = v.get().real #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) - fplot.write_vtk_file( - "potential.vts", - [ - ("u", u), - ("v", v), - ("indicator", indicator), - ] - ) + fplot.write_vtk_file("potential.vts", [ + ("u", u), + ("v", v), + ("indicator", indicator), + ]) # }}} -- GitLab From 0a249aa5825d705e79a2a0ec8a6da5f2f2485622 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 9 Jun 2017 15:42:58 -0400 Subject: [PATCH 11/11] current state for merge --- examples/cahn-hilliard.py | 258 ++++++++++++++++++++++++-------------- 1 file changed, 167 insertions(+), 91 deletions(-) diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index 4cbd1680..2308ce8e 100644 --- a/examples/cahn-hilliard.py +++ b/examples/cahn-hilliard.py @@ -16,6 +16,15 @@ import pytential.symbolic.primitives as p from sumpy.kernel import ExpressionKernel import loopy as lp + +# clean up the mess +def clean_file(filename): + import os + try: + os.remove(filename) + except OSError: + pass + # {{{ set some constants for use below # {{{ all kinds of orders @@ -58,7 +67,9 @@ import pymbolic as pmbl x = pmbl.var("x") y = pmbl.var("y") -# FIXME: modify pymbolic to use tanh function +# b, c are known constants + +# TODO: modify pymbolic to use tanh function # phi = tanh(x / sqrt (2 * epsilon)) phi = x**2 + x * y + 2 * y**2 phi3 = phi**3 @@ -69,27 +80,35 @@ laphi3 = pmbl.differentiate(pmbl.differentiate(phi3, 'x'), 'x') + \ f1_expr = c * phi - (1 + s) / epsilon**2 * laphi + 1 / epsilon**2 * laphi3 f2_expr = (phi3 - (1 + s) * phi) / epsilon**2 - def f1_func(x, y): return pmbl.evaluate(f1_expr, {"x": x, "y": y}) - def f2_func(x, y): return pmbl.evaluate(f2_expr, {"x": x, "y": y}) - def initial_phi(x, y): return pmbl.evaluate(phi, {"x": x, "y": y}) +# }}} +# {{{ manufactured rhs such that the solution is phi given above + +lalaphi = pmbl.differentiate(pmbl.differentiate(laphi, 'x'), 'x') + \ + pmbl.differentiate(pmbl.differentiate(laphi, 'y'), 'y') +mu = - epsilon * laphi - phi / epsilon + phi3 / epsilon + +rhs1_expr = lalaphi - b * laphi + c * phi +rhs2_expr = mu / epsilon + laphi - b * phi +def rhs1_func(x, y): + return pmbl.evaluate(rhs1_expr, {"x": x, "y": y}) -#def initial_phi(x, y): -# return np.tanh(x / np.sqrt(2. * initial_epsilon)) +def rhs2_func(x, y): + return pmbl.evaluate(rhs2_expr, {"x": x, "y": y}) # }}} # }}} # {{{ a kernel class for G0 -# FIXME: will the expressions work when lambda is complex? +# TODO: will the expressions work when lambda is complex? # (may need conversion from 1j to var("I")) from sumpy.kernel import ExpressionKernel @@ -158,7 +177,6 @@ class ShidongKernel(ExpressionKernel): # }}} - # {{{ extended kernel getters class HankelBasedKernel(ExpressionKernel): def prepare_loopy_kernel(self, loopy_knl): @@ -176,11 +194,14 @@ def get_extkernel_for_G0(lambda1, lambda2): from pymbolic import var d = make_sym_vector("d", 3) - r2 = pymbolic_real_norm_2(d[:-1]) - expr = var("hankel_1")(0, var("I") * lambda1 * r2 - + var("I") * d[-1]**2) \ - - var("hankel_1")(0, var("I") * lambda2 * r2 - + var("I") * d[-1]**2) + # r2 = pymbolic_real_norm_2(d[:-1]) + r3 = pymbolic_real_norm_2(d) + #expr = var("hankel_1")(0, var("I") * lambda1 * r2 + # + var("I") * d[-1]**2) \ + # - var("hankel_1")(0, var("I") * lambda2 * r2 + # + var("I") * d[-1]**2) + expr = var("hankel_1")(0, var("I") * lambda1 * r3) \ + - var("hankel_1")(0, var("I") * lambda2 * r3) scaling = 1. / (4. * var("I") * (lambda1**2 - lambda2**2)) return HankelBasedKernel( @@ -193,8 +214,10 @@ def get_extkernel_for_G1(lamb): from pymbolic import var d = make_sym_vector("d", 3) - r2 = pymbolic_real_norm_2(d[:-1]) - expr = var("hankel_1")(0, var("I") * lamb * r2 + var("I") * d[-1]**2) + #r2 = pymbolic_real_norm_2(d[:-1]) + r3 = pymbolic_real_norm_2(d) + #expr = var("hankel_1")(0, var("I") * lamb * r2 + var("I") * d[-1]**2) + expr = var("hankel_1")(0, var("I") * lamb * r3) scaling = -var("I") / 4. return HankelBasedKernel( @@ -292,30 +315,39 @@ def main(): # }}} + # {{{ visualizers + + from meshmode.discretization.visualization import make_visualizer + vol_vis = make_visualizer(queue, vol_discr, 20) + bdry_vis = make_visualizer(queue, density_discr, 20) + + # }}} + # {{{ setup operator and potentials - #print("-- setup Cahn-Hilliard operator") from pytential.symbolic.pde.cahn_hilliard import CahnHilliardOperator chop = CahnHilliardOperator(b=b, c=c) unk = chop.make_unknown("sigma") bound_op = bind(qbx, chop.operator(unk)) - #print ("-- construct kernels") yukawa_2d_in_3d_kernel_1 = get_extkernel_for_G1(chop.lambdas[0]) shidong_2d_in_3d_kernel = get_extkernel_for_G0(chop.lambdas[0], chop.lambdas[1]) - #print("-- construct layer potentials") from sumpy.qbx import LayerPotential from sumpy.expansion.local import LineTaylorLocalExpansion + layer_pot_v0f1 = LayerPotential(cl_ctx, [ LineTaylorLocalExpansion(shidong_2d_in_3d_kernel, order=vol_qbx_order) ]) layer_pot_v1f1 = LayerPotential(cl_ctx, [ LineTaylorLocalExpansion( yukawa_2d_in_3d_kernel_1, order=vol_qbx_order) - ]) + ]) + + #print(layer_pot_v0f1.get_kernel()) + # print(layer_pot_v1f1.get_kernel()) # }}} # {{{ setup for volume integral @@ -341,6 +373,9 @@ def main(): sources = cl.array.zeros(queue, (3, ) + vol_x.shape[1:], vol_x.dtype) sources[:2] = vol_x + #np.set_printoptions(threshold=np.inf) + #print(sources) + vol_weights = bind( vol_discr, p.area_element(mesh.ambient_dim, mesh.dim) * p.QWeight())(queue) @@ -366,16 +401,9 @@ def main(): time += delta_t output_vts_filename = get_vts_filename(timestep_number) - # a manufactured f1 function - #x_sin_factor = 30 - #y_sin_factor = 10 - #def f1_func(x, y): - # return 0.1 * cl.clmath.sin(x_sin_factor*x) \ - # * cl.clmath.sin(y_sin_factor*y) - # get f1 to compute strengths - f1 = f1_func(vol_x[0], vol_x[1]) - f2 = f2_func(vol_x[0], vol_x[1]) + f1 = rhs1_func(vol_x[0], vol_x[1]) + f2 = rhs2_func(vol_x[0], vol_x[1]) evt, (vol_pot_v0f1, ) = layer_pot_v0f1( queue, @@ -391,6 +419,34 @@ def main(): sources=sources.reshape(3, vol_discr.nnodes), strengths=((vol_weights * f1).reshape(vol_discr.nnodes), )) + plot_volume_potentials = True + if plot_volume_potentials: + + evt, (vol_pot_v0id, ) = layer_pot_v0f1( + queue, + targets=targets.reshape(3, vol_discr.nnodes), + centers=centers, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=((vol_weights * 1).reshape(vol_discr.nnodes), )) + + evt, (vol_pot_v1id, ) = layer_pot_v1f1( + queue, + targets=targets.reshape(3, vol_discr.nnodes), + centers=centers, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=((vol_weights * 1).reshape(vol_discr.nnodes), )) + + #print(vol_pot_v1f1.get().real[0:20]) + #print(vol_pot_v1f1.get().real[0:20]) + #print(vol_pot_v0id.get().real[0:20]) + #print(vol_pot_v1id.get().real[0:20]) + + clean_file("ch-volume.vtu") + vol_vis.write_vtk_file("ch-volume.vtu", [ + ("vol_pot_v0f1", vol_pot_v0f1), ("vol_pot_v1f1", vol_pot_v1f1), + ("vol_pot_v0id", vol_pot_v0id), ("vol_pot_v1id", vol_pot_v1id) + ]) + # {{{ fix rhs and solve nodes = density_discr.nodes().with_queue(queue) @@ -400,7 +456,7 @@ def main(): return cl.clmath.cos(5 * cl.clmath.atan2(y, x)) bc = sym.make_obj_array([ - # FIXME: Realistic BC + # TODO: Realistic BC 5 + g(nodes), 3 - g(nodes), ]) @@ -418,11 +474,56 @@ def main(): # }}} + # }}} + + # {{{ postprocess/visualize + + immersed_visualization = False + if immersed_visualization == True: + from sumpy.visualization import FieldPlotter + fplot = FieldPlotter(np.zeros(2), extent=1.5, npoints=500) + + targets = cl.array.to_device(queue, fplot.points) + + qbx_stick_out = qbx.copy(target_stick_out_factor=0.05) + indicator_qbx = qbx_stick_out.copy(qbx_order=2) + + from sumpy.kernel import LaplaceKernel + ones_density = density_discr.zeros(queue) + ones_density.fill(1) + indicator = bind((indicator_qbx, PointsTarget(targets)), + sym.D(LaplaceKernel(2), sym.var("sigma")))( + queue, sigma=ones_density).get() + + clean_file("failed-targets.vts") + clean_file("potential.vts") + + try: + u, v = bind((qbx_stick_out, PointsTarget(targets)), + chop.representation(unk))( + queue, sigma=sigma) + except QBXTargetAssociationFailedException as e: + fplot.write_vtk_file("failed-targets.vts", + [("failed", e.failed_target_flags.get(queue))]) + raise + u = u.get().real + v = v.get().real + + #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) + fplot.write_vtk_file("potential.vts", [ + ("u", u), + ("v", v), + ("indicator", indicator), + ]) + + # }}} + # {{{ check pde def check_pde(): from sumpy.point_calculus import CalculusPatch - vec_h = [1e-1, 1e-2, 1e-3, 1e-4] + vec_h = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7] + #vec_h = [1e-1] vec_ru = [] vec_rv = [] vec_rp = [] @@ -431,6 +532,8 @@ def main(): vec_rf2 = [] vec_rutld = [] vec_rvtld = [] + vec_rv1id = [] # V1[id] + vec_rv0id = [] # V0[id] for dx in vec_h: cp = CalculusPatch(np.zeros(2), order=4, h=dx) targets = cl.array.to_device(queue, cp.points) @@ -447,7 +550,8 @@ def main(): # Check for homogeneous PDEs for u and v vec_ru.append( la.norm( - cp.laplace(lap_u) - chop.b * cp.laplace(u) + chop.c * u)) + cp.laplace(lap_u) - chop.b * cp.laplace(u) + + chop.c * u)) vec_rv.append(la.norm(v + cp.laplace(u) - chop.b * u)) # Check for inhomogeneous PDEs for phi and mu @@ -472,6 +576,16 @@ def main(): sources=sources.reshape(3, vol_discr.nnodes), strengths=((vol_weights * f1).reshape(vol_discr.nnodes), )) v0f1 = v0f1.get().real + print(v0f1) + + evt, (v0id, ) = layer_pot_v0f1( + queue, + targets=targets_in_3d, + centers=centers_in_3d, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=((vol_weights * 1).reshape(vol_discr.nnodes), )) + v0id = v0id.get().real + print(v0id) evt, (v1f1, ) = layer_pot_v1f1( queue, @@ -480,6 +594,23 @@ def main(): sources=sources.reshape(3, vol_discr.nnodes), strengths=((vol_weights * f1).reshape(vol_discr.nnodes), )) v1f1 = v1f1.get().real + print(v1f1) + + evt, (v1id, ) = layer_pot_v1f1( + queue, + targets=targets_in_3d, + centers=centers_in_3d, + sources=sources.reshape(3, vol_discr.nnodes), + strengths=((vol_weights * 1).reshape(vol_discr.nnodes), )) + v1id = v1id.get().real + print(v1id) + + # FIXME: Why this is equal to 0 (instead of 1) + vec_rv1id.append( + la.norm(cp.laplace(v1id) - chop.lambdas[1]**2 * v1id + - 1)) + vec_rv0id.append( + la.norm(cp.laplace(v0id) - chop.lambdas[1]**2 * v0id - v1id)) f14ck = cl.array.to_device(queue, f1_func(cp.points[0], cp.points[1])) @@ -506,15 +637,8 @@ def main(): # FIXME: Not passing this check. (bugs in volume integral?) vec_rutld.append( la.norm( - cp.laplace(lap_utild) - chop.b * lap_utild + - chop.c * utild - f14ck)) - print(f14ck) - print(f1[0:10]) - print(cp.laplace(lap_utild)) - print((lap_utild)) - print((utild)) - print(chop.b) - print(chop.c) + cp.laplace(lap_utild) - chop.b * lap_utild + chop.c * utild + - f14ck)) vec_rvtld.append( la.norm(vtild + cp.laplace(utild) - chop.b * utild - f24ck)) @@ -542,6 +666,8 @@ def main(): ["residual_v"] + vec_rv, ["residual_f1"] + vec_rf1, ["residual_f2"] + vec_rf2, + ["residual_v1id"] + vec_rv1id, + ["residual_v0id"] + vec_rv0id, ["residual_u_tld"] + vec_rutld, ["residual_v_tld"] + vec_rvtld, ["residual_phi"] + vec_rp, @@ -550,60 +676,10 @@ def main(): file=f) 1 / 0 - check_pde() - - # }}} - - # {{{ postprocess/visualize - - from sumpy.visualization import FieldPlotter - fplot = FieldPlotter(np.zeros(2), extent=1.5, npoints=500) - - targets = cl.array.to_device(queue, fplot.points) - - qbx_stick_out = qbx.copy(target_stick_out_factor=0.05) - indicator_qbx = qbx_stick_out.copy(qbx_order=2) - - from sumpy.kernel import LaplaceKernel - ones_density = density_discr.zeros(queue) - ones_density.fill(1) - indicator = bind((indicator_qbx, PointsTarget(targets)), - sym.D(LaplaceKernel(2), sym.var("sigma")))( - queue, sigma=ones_density).get() - - # clean up the mess - def clean_file(filename): - import os - try: - os.remove(filename) - except OSError: - pass - - clean_file("failed-targets.vts") - clean_file("potential.vts") - - try: - u, v = bind((qbx_stick_out, PointsTarget(targets)), - chop.representation(unk))( - queue, sigma=sigma) - except QBXTargetAssociationFailedException as e: - fplot.write_vtk_file("failed-targets.vts", - [("failed", e.failed_target_flags.get(queue))]) - raise - u = u.get().real - v = v.get().real - - #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) - fplot.write_vtk_file("potential.vts", [ - ("u", u), - ("v", v), - ("indicator", indicator), - ]) + # check_pde() # }}} -# }}} - if __name__ == "__main__": main() -- GitLab