...
 
Commits (3)
......@@ -158,7 +158,7 @@ def cse(kernel):
writer_map = writer_map_with_aliasing(kernel)
pdelab_variables = set("r r_n r_s jac jac_s_s jac_s_n jac_n_n jac_n_s".split())
pdelab_variables = set("r r_n r_s jac jac_s_s jac_s_n jac_n_n jac_n_s x x_n x_s z z_n z_s".split())
name_generator = UniqueNameGenerator(kernel.all_variable_names() | pdelab_variables)
......
......@@ -1010,10 +1010,12 @@ def generate_residual_kernels(form, original_form):
operator_kernels = {}
if refinement:
from ufl import Form
bdry_integrals = form.integrals_by_type('exterior_facet')
extruded_bdry_form = (Form([i.reconstruct(integral_type='bottom_facet') for i in bdry_integrals]) +
Form([i.reconstruct(integral_type='top_facet') for i in bdry_integrals]))
form = form + extruded_bdry_form
from meshstructure import MeshExtrusion
if isinstance(refinement, MeshExtrusion):
bdry_integrals = form.integrals_by_type('exterior_facet')
extruded_bdry_form = (Form([i.reconstruct(integral_type='bottom_facet') for i in bdry_integrals]) +
Form([i.reconstruct(integral_type='top_facet') for i in bdry_integrals]))
form = form + extruded_bdry_form
used_measures = {"volume": ("cell", "interior_facet", "bottom_facet", "top_facet"),
"skeleton": ("interior_facet",),
......@@ -1034,6 +1036,59 @@ def generate_residual_kernels(form, original_form):
pattern_baseclass()
enum_alpha()
# Maybe add numerical differentiation
if get_form_option("numerical_jacobian"):
# Include headers for numerical methods
include_file("dune/pdelab/localoperator/defaultimp.hh", filetag="operatorfile")
# Numerical jacobian base class
_, loptype = class_type_from_cache("operator")
from dune.codegen.pdelab.signatures import ufl_measure_to_pdelab_measure
which = pdelab_kernel.capitalize()
base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype),
classtag="operator")
# Numerical jacobian initializer list
ini = name_initree_member()
ini_constructor = name_initree_constructor_param()
initializer_list("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype),
["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(
ini_constructor, ini, which.lower())],
classtag="operator",
)
# In the case of matrix free operator evaluation we need jacobian apply methods
if get_form_option("generate_jacobian_apply"):
from dune.codegen.pdelab.driver import is_linear
if is_linear(original_form):
# Numeical jacobian apply base class
base_class("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype),
classtag="operator")
# Numerical jacobian apply initializer list
initializer_list(
"Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype),
["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor,
ini,
which.lower())],
classtag="operator",
)
else:
# Numerical nonlinear jacobian apply base class
base_class("Dune::PDELab::NumericalNonlinearJacobianApply{}<{}>".format(which,
loptype),
classtag="operator")
# Numerical nonlinear jacobian apply initializer list
initializer_list(
"Dune::PDELab::NumericalNonlinearJacobianApply{}<{}>".format(which,
loptype),
["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor,
ini,
which.lower())],
classtag="operator",
)
operator_kernels[(pdelab_kernel, 'residual')] = kernel
return operator_kernels
# Generate the necessary residual methods
......@@ -1162,11 +1217,12 @@ def generate_jacobian_kernels(form, original_form):
if refinement:
with global_context(form_type="jacobian_apply"):
from ufl import Form
bdry_integrals = jac_apply_form.integrals_by_type('exterior_facet')
extruded_bdry_form = (Form([i.reconstruct(integral_type='bottom_facet') for i in bdry_integrals]) +
Form([i.reconstruct(integral_type='top_facet') for i in bdry_integrals]))
jac_apply_form = jac_apply_form + extruded_bdry_form
from meshstructure import MeshExtrusion
if isinstance(refinement, MeshExtrusion):
bdry_integrals = jac_apply_form.integrals_by_type('exterior_facet')
extruded_bdry_form = (Form([i.reconstruct(integral_type='bottom_facet') for i in bdry_integrals]) +
Form([i.reconstruct(integral_type='top_facet') for i in bdry_integrals]))
jac_apply_form = jac_apply_form + extruded_bdry_form
used_measures = {"volume": ("cell", "interior_facet", "bottom_facet", "top_facet"),
"skeleton": ("interior_facet",),
"boundary": ("exterior_facet",)}
......
from meshstructure import *
from meshstructure.unstructured import UnstructuredHyperCube
quadrilateral = UnstructuredHyperCube(2)
dx = dx(metadata={"refinement": HyperCubeRefinement(quadrilateral, 4, 4), "apply_cse": True})
x = SpatialCoordinate(cell)
g = sum([c**2 for c in x])
f = sum([-2 for _ in x])
V = FiniteElement("CG", cell, degree)
V = FiniteElement("CG", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
......
__name = blockstructured_poisson_neumann_{__exec_suffix}
__exec_suffix = numdiff, symdiff | expand num
__name = blockstructured_poisson_neumann
cells = 20 20
extension = 1. 1.
......@@ -13,7 +12,8 @@ extension = vtu
compare_l2errorsquared = 1e-8
[formcompiler.r]
numerical_jacobian = 1, 0 | expand num
matrix_free = 1
generate_jacobians = 0
blockstructured = 1
number_of_blocks = 4
geometry_mixins = blockstructured_equidistant
from meshstructure import *
from meshstructure.unstructured import UnstructuredHyperCube
cell = quadrilateral
x = SpatialCoordinate(cell)
quadrilateral = UnstructuredHyperCube(2)
dx = dx(metadata={"refinement": HyperCubeRefinement(quadrilateral, 4, 4), "apply_cse": True})
ds = ds(metadata={"refinement": HyperCubeRefinement(quadrilateral, 4, 4), "apply_cse": True})
c = (0.5-x[0])**2 + (0.5-x[1])**2
g = exp(-1.*c)
f = 4*(1.-c)*g
......
from meshstructure import *
from meshstructure.unstructured import UnstructuredHyperCube
x = SpatialCoordinate(cell)
quadrilateral = UnstructuredHyperCube(2)
dx = dx(metadata={"refinement": HyperCubeRefinement(quadrilateral, 8, 8), "apply_cse": True})
I = Identity(dim)
A = as_matrix([[x[i]*x[j] + I[i,j] for j in range(dim)] for i in range(dim)])
g = sum([c**2 for c in x])
......
from meshstructure import *
from meshstructure.unstructured import UnstructuredHyperCube
cell = quadrilateral
x = SpatialCoordinate(cell)
dirichlet = conditional(x[0] + x[1] < 1.0, 1, 0)
dx = dx(metadata={"levels":4, "height":1., "quadrature_order": 2})
dS = dS(metadata={"levels":4, "height":1., "quadrature_order": 2})
ds = ds(metadata={"levels":4, "height":1., "quadrature_order": 2})
line = UnstructuredHyperCube(1)
metadata = {"refinement": MeshExtrusion(line, 10), "levels": 10, "height": 1., "quadrature_order": 2, "apply_cse": True}
dx = dx(metadata=metadata)
dS = dS(metadata=metadata)
ds = ds(metadata=metadata)
dsd = ds(subdomain_data=dirichlet, subdomain_id=1)
dso = ds(subdomain_data=dirichlet, subdomain_id=0)
V = FiniteElement("DG", cell, 0)
V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
beta = as_vector((1., 1.))
beta = as_vector((1., 0.))
n = FacetNormal(cell)('+')
beta_n_pos = 0.5 * (inner(beta, n) + abs(inner(beta, n)))
beta_n_neg = 0.5 * (-1 * inner(beta, n) + abs(inner(beta, n)))
def numerical_flux(normal, inside, outside):
return conditional(inner(beta, n) > 0, inside, outside)*inner(beta, n)
def numerical_flux(inside, outside):
return beta_n_pos * inside - beta_n_neg * outside
mass = u*v*dx
r = 1.*u*inner(beta, grad(v))*dx - numerical_flux(n, u('+'), u('-'))*jump(v)*dS - inner(beta, n)*u*v*dso
r = 1.*u*inner(beta, grad(v))*dx - numerical_flux(u('+'), u('-'))*jump(v)*dS - inner(beta, n)*u*v*dso - inner(beta,n)*v*dsd