From f77a9fbee1b5685092cead4426297ff37f70e20a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Thu, 24 Nov 2016 17:22:42 +0100
Subject: [PATCH] 3D sumfactorized poisson
---
python/dune/perftool/loopy/buffer.py | 10 +++++++++-
python/dune/perftool/sumfact/basis.py | 20 +++++++++++++++-----
python/dune/perftool/sumfact/sumfact.py | 17 ++++++++++++-----
test/sumfact/mass/CMakeLists.txt | 5 +++++
test/sumfact/mass/mass_3d.mini | 19 +++++++++++++++++++
test/sumfact/mass/mass_3d.ufl | 10 ++++++++++
test/sumfact/poisson/CMakeLists.txt | 11 +++++++++++
test/sumfact/poisson/poisson_3d_order1.mini | 20 ++++++++++++++++++++
test/sumfact/poisson/poisson_3d_order1.ufl | 10 ++++++++++
test/sumfact/poisson/poisson_3d_order2.mini | 20 ++++++++++++++++++++
test/sumfact/poisson/poisson_3d_order2.ufl | 10 ++++++++++
11 files changed, 141 insertions(+), 11 deletions(-)
create mode 100644 test/sumfact/mass/mass_3d.mini
create mode 100644 test/sumfact/mass/mass_3d.ufl
create mode 100644 test/sumfact/poisson/poisson_3d_order1.mini
create mode 100644 test/sumfact/poisson/poisson_3d_order1.ufl
create mode 100644 test/sumfact/poisson/poisson_3d_order2.mini
create mode 100644 test/sumfact/poisson/poisson_3d_order2.ufl
diff --git a/python/dune/perftool/loopy/buffer.py b/python/dune/perftool/loopy/buffer.py
index e7207c2b..d7e6c97a 100644
--- a/python/dune/perftool/loopy/buffer.py
+++ b/python/dune/perftool/loopy/buffer.py
@@ -1,5 +1,6 @@
from dune.perftool.error import PerftoolLoopyError
from dune.perftool.generation import (generator_factory,
+ get_global_context_value,
temporary_variable,
)
@@ -33,10 +34,17 @@ class FlipFlopBuffer(object):
name = "{}_{}".format(self.identifier, self._counter)
self._counter = self._counter + 1
+ # Get geometric dimension
+ formdata = get_global_context_value('formdata')
+ dim = formdata.geometric_dimension
+
+ # TODO: doc!
+ storage_shape = (self.base_storage_size,) + (1,)*(dim-1)
+
# Construct the temporary and return it
temporary_variable(name,
base_storage=base,
- storage_shape=(self.base_storage_size, 1),
+ storage_shape=storage_shape,
managed=True,
**kwargs
)
diff --git a/python/dune/perftool/sumfact/basis.py b/python/dune/perftool/sumfact/basis.py
index 0af2d38c..03665d47 100644
--- a/python/dune/perftool/sumfact/basis.py
+++ b/python/dune/perftool/sumfact/basis.py
@@ -3,6 +3,7 @@
NB: Basis evaluation is only needed for the trial function argument in jacobians, as the
multiplication with the test function is part of the sum factorization kernel.
"""
+
from dune.perftool.generation import (backend,
domain,
get_counter,
@@ -69,13 +70,15 @@ def sumfact_evaluate_coefficient_gradient(element, name, restriction, component)
a_matrices[i] = dtheta_matrix
a_matrices = tuple(a_matrices)
- # Do stage 1
buffer_name = name_sumfact_base_buffer()
initialize_buffer(buffer_name,
base_storage_size=product(max(mat.rows, mat.cols) for mat in a_matrices),
num=2
)
+ # Add a sum factorization kernel that implements the
+ # evaluation of the gradients of basis functions at quadrature
+ # points (stage 1)
insn_dep = setup_theta(element, restriction, component, a_matrices, buffer_name)
var, _ = sum_factorization_kernel(a_matrices,
buffer_name,
@@ -110,20 +113,26 @@ def pymbolic_trialfunction_gradient(element, restriction, component):
@kernel_cached
def pymbolic_trialfunction(element, restriction, component):
+ # Get geometric dimension
+ formdata = get_global_context_value('formdata')
+ dim = formdata.geometric_dimension
+
+ # Setup sumfactorization
theta = name_theta()
rows = quadrature_points_per_direction()
cols = basis_functions_per_direction()
a_matrix = AMatrix(theta, rows, cols)
- # TODO: this is 2D only
- a_matrices = (a_matrix, a_matrix)
- buffer_name = name_sumfact_base_buffer()
+ a_matrices = (a_matrix,)*dim
- # Do stage 1
+ # Flip flop buffers for sumfactorization
+ buffer_name = name_sumfact_base_buffer()
initialize_buffer(buffer_name,
base_storage_size=product(max(mat.rows, mat.cols) for mat in a_matrices),
num=2
)
+ # Add a sum factorization kernel that implements the evaluation of
+ # the basis functions at quadrature points (stage 1)
insn_dep = setup_theta(element, restriction, component, a_matrices, buffer_name)
var, _ = sum_factorization_kernel(a_matrices,
buffer_name,
@@ -193,6 +202,7 @@ def evaluate_reference_gradient(element, name, restriction):
inames = lfs_inames(element, restriction)
assert(len(quad_inames) == len(inames))
+ # Matrices for sumfactorization
theta = name_theta()
dtheta = name_dtheta()
diff --git a/python/dune/perftool/sumfact/sumfact.py b/python/dune/perftool/sumfact/sumfact.py
index 36240228..04a18a7b 100644
--- a/python/dune/perftool/sumfact/sumfact.py
+++ b/python/dune/perftool/sumfact/sumfact.py
@@ -32,6 +32,7 @@ from dune.perftool.pdelab.restriction import restricted_name
from dune.perftool.pdelab.spaces import name_lfs
from dune.perftool.sumfact.amatrix import (AMatrix,
quadrature_points_per_direction,
+ name_dtheta_transposed,
basis_functions_per_direction,
name_theta,
name_theta_transposed,
@@ -124,7 +125,6 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
# If this is a gradient we need different matrices
if accterm.argument.index:
- from dune.perftool.sumfact.amatrix import name_dtheta_transposed
dtheta_transposed = name_dtheta_transposed()
rows = basis_functions_per_direction()
cols = quadrature_points_per_direction()
@@ -135,13 +135,14 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
a_matrices = tuple(a_matrices)
else:
- a_matrices = (theta_matrix, theta_matrix)
+ a_matrices = (theta_matrix,)*dim
+ # Initialize a base storage for this buffer and get a temporay pointing to it
temp = initialize_buffer(buf,
base_storage_size=product(max(mat.rows, mat.cols) for mat in a_matrices),
num=2
).get_temporary(shape=(quadrature_points_per_direction(),) * dim,
- dim_tags="f,f")
+ dim_tags=",".join(['f']*dim))
# Replace gradient iname with correct index for assignement
replace_dict = {}
@@ -162,7 +163,7 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
contribution_ids.append(contrib_dep)
# Add a sum factorization kernel that implements the multiplication
- # with the test function (step 3)
+ # with the test function (stage 3)
result, insn_dep = sum_factorization_kernel(a_matrices,
buf,
insn_dep=frozenset({contrib_dep}),
@@ -274,9 +275,15 @@ def sum_factorization_kernel(a_matrices, buf, insn_dep=frozenset({}), additional
)
})
+ # Get geometric dimension
+ formdata = get_global_context_value('formdata')
+ dim = formdata.geometric_dimension
+
+ out_shape = tuple(mat.rows for mat in a_matrices)
+ dim_tags = ",".join(['f']*dim)
out = get_buffer_temporary(buf,
shape=out_shape,
- dim_tags="f,f")
+ dim_tags=dim_tags)
silenced_warning('read_no_write({})'.format(out))
return out, insn_dep
diff --git a/test/sumfact/mass/CMakeLists.txt b/test/sumfact/mass/CMakeLists.txt
index a3357bef..a640f020 100644
--- a/test/sumfact/mass/CMakeLists.txt
+++ b/test/sumfact/mass/CMakeLists.txt
@@ -3,3 +3,8 @@ dune_add_formcompiler_system_test(UFLFILE mass.ufl
BASENAME sumfact_mass
INIFILE mass.mini
)
+
+dune_add_formcompiler_system_test(UFLFILE mass_3d.ufl
+ BASENAME sumfact_mass_3d
+ INIFILE mass_3d.mini
+ )
diff --git a/test/sumfact/mass/mass_3d.mini b/test/sumfact/mass/mass_3d.mini
new file mode 100644
index 00000000..17da3e4f
--- /dev/null
+++ b/test/sumfact/mass/mass_3d.mini
@@ -0,0 +1,19 @@
+__name = sumfact_mass_3d_{__exec_suffix}
+__exec_suffix = {diff_suffix}_{vec_suffix}
+
+diff_suffix = numdiff, symdiff | expand num
+vec_suffix = vec, nonvec | expand vec
+
+cells = 1 1 1
+extension = 1. 1. 1.
+
+printmatrix = true
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 1, 0 | expand num
+vectorize = 1, 0 | expand vec
+sumfact = 1
diff --git a/test/sumfact/mass/mass_3d.ufl b/test/sumfact/mass/mass_3d.ufl
new file mode 100644
index 00000000..bd56d296
--- /dev/null
+++ b/test/sumfact/mass/mass_3d.ufl
@@ -0,0 +1,10 @@
+cell = "hexahedron"
+
+V = FiniteElement("DG", cell, 1)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+r = u * v * dx
+
+forms = [r]
diff --git a/test/sumfact/poisson/CMakeLists.txt b/test/sumfact/poisson/CMakeLists.txt
index 74cce96e..70a9791a 100644
--- a/test/sumfact/poisson/CMakeLists.txt
+++ b/test/sumfact/poisson/CMakeLists.txt
@@ -11,6 +11,17 @@ dune_add_formcompiler_system_test(UFLFILE poisson_order2.ufl
INIFILE poisson_order2.mini
)
+# 3D Test cases
+dune_add_formcompiler_system_test(UFLFILE poisson_3d_order1.ufl
+ BASENAME sumfact_poisson_3d_order1
+ INIFILE poisson_3d_order1.mini
+ )
+dune_add_formcompiler_system_test(UFLFILE poisson_3d_order2.ufl
+ BASENAME sumfact_poisson_3d_order2
+ INIFILE poisson_3d_order2.mini
+ )
+
+
# # 2. Poisson Test Case: DG, f + pure dirichlet
# dune_add_formcompiler_system_test(UFLFILE poisson_dg.ufl
# BASENAME sumfact_poisson_dg
diff --git a/test/sumfact/poisson/poisson_3d_order1.mini b/test/sumfact/poisson/poisson_3d_order1.mini
new file mode 100644
index 00000000..6250546e
--- /dev/null
+++ b/test/sumfact/poisson/poisson_3d_order1.mini
@@ -0,0 +1,20 @@
+__name = sumfact_poisson_3d_order1_{__exec_suffix}
+__exec_suffix = {diff_suffix}_{vec_suffix}
+
+diff_suffix = numdiff, symdiff | expand num
+vec_suffix = vec, nonvec | expand vec
+
+cells = 8 8 8
+extension = 1. 1. 1.
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_ref
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 1, 0 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-4
+sumfact = 1
+vectorize = 1, 0 | expand vec
diff --git a/test/sumfact/poisson/poisson_3d_order1.ufl b/test/sumfact/poisson/poisson_3d_order1.ufl
new file mode 100644
index 00000000..4de9c239
--- /dev/null
+++ b/test/sumfact/poisson/poisson_3d_order1.ufl
@@ -0,0 +1,10 @@
+cell = "hexahedron"
+
+f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());", cell=cell)
+g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", cell=cell)
+
+V = FiniteElement("CG", cell, 1, dirichlet_expression=g)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+forms = [(inner(grad(u), grad(v)) - f*v)*dx]
diff --git a/test/sumfact/poisson/poisson_3d_order2.mini b/test/sumfact/poisson/poisson_3d_order2.mini
new file mode 100644
index 00000000..2db0990d
--- /dev/null
+++ b/test/sumfact/poisson/poisson_3d_order2.mini
@@ -0,0 +1,20 @@
+__name = sumfact_poisson_3d_order2_{__exec_suffix}
+__exec_suffix = {diff_suffix}_{vec_suffix}
+
+diff_suffix = numdiff, symdiff | expand num
+vec_suffix = vec, nonvec | expand vec
+
+cells = 8 8 8
+extension = 1. 1. 1.
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_ref
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 1, 0 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-4
+sumfact = 1
+vectorize = 1, 0 | expand vec
diff --git a/test/sumfact/poisson/poisson_3d_order2.ufl b/test/sumfact/poisson/poisson_3d_order2.ufl
new file mode 100644
index 00000000..66870b07
--- /dev/null
+++ b/test/sumfact/poisson/poisson_3d_order2.ufl
@@ -0,0 +1,10 @@
+cell = "hexahedron"
+
+f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());", cell=cell)
+g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", cell=cell)
+
+V = FiniteElement("CG", cell, 2, dirichlet_expression=g)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+forms = [(inner(grad(u), grad(v)) - f*v)*dx]
--
GitLab