Skip to content
Snippets Groups Projects
Commit 2f567ddb authored by Dominic Kempf's avatar Dominic Kempf
Browse files

Symbolic jacobian is now symmetric too, (but still different)

parent da822b6a
No related branches found
No related tags found
No related merge requests found
......@@ -165,45 +165,27 @@ def transform_accumulation_term(term, measure, subdomain_id):
expr_tv = temporary_variable(expr_tv_name)
# The data that is used to collect the arguments for the accumulate function
accumargs = []
residual_shape = {}
accumargs = [None] * (2 * len(test_ma))
residual_shape = [None] * len(test_ma)
arg_restr = [None] * len(test_ma)
# Determine whether the inames in the jacobian case clash. This is the case, when
# we have a galerkin method.
galerkin = False
if len(test_ma) == 2:
galerkin = (test_ma[0].argexpr.element() == test_ma[1].argexpr.element())
# Generate the code for the modified arguments:
for arg in test_ma:
from dune.perftool.pdelab.argument import pymbolic_argument
from dune.perftool.pdelab.basis import name_lfs
accumargs.append(name_lfs(arg.argexpr.element()))
if galerkin and arg.argexpr.number() == 1:
iname = lfs_iname(arg.argexpr.element(), context="arg")
else:
iname = lfs_iname(arg.argexpr.element())
accumargs.append(iname)
# Add this iname to the set of used inames for dependencies later on
acc_inames = acc_inames.union(frozenset({iname}))
for ma in test_ma:
count = ma.argexpr.number()
lfs = name_lfs(ma.argexpr.element())
lfsi = lfs_iname(ma.argexpr.element(), argnumber=count)
# Determine the shape
residual_shape[arg.argexpr.number()] = name_lfs_bound(name_lfs(arg.argexpr.element()))
accumargs[2 * count] = lfs
accumargs[2 * count + 1] = lfsi
# Determine the restriction for later
arg_restr[arg.argexpr.number()] = arg.restriction
arg_restr[count] = ma.restriction
residual_shape[count] = name_lfs_bound(lfs)
from dune.perftool.pdelab.argument import name_accumulation_variable
accumvar = name_accumulation_variable(arg_restr)
# The residual/the jacobian should be represented through a loopy global argument
# TODO this seems still a bit hacky, esp. w.r.t. systems
shape = tuple(v for k, v in sorted(residual_shape.items(), key=lambda (k, v): k))
globalarg(accumvar, shape=shape)
globalarg(accumvar, shape=tuple(residual_shape))
from dune.perftool.pdelab.quadrature import name_factor
factor = name_factor()
......
......@@ -36,9 +36,9 @@ def pymbolic_testfunction(ma):
assert bool(ma.index) == ma.grad
if ma.grad:
return Subscript(Variable(name_testfunction_gradient(ma)), (Variable(lfs_iname(ma.argexpr.element())), Variable(name_index(ma.index))))
return Subscript(Variable(name_testfunction_gradient(ma)), (Variable(lfs_iname(ma.argexpr.element(), argnumber=ma.argexpr.number())), Variable(name_index(ma.index))))
else:
return Subscript(Variable(name_testfunction(ma)), Variable(lfs_iname(ma.argexpr.element())))
return Subscript(Variable(name_testfunction(ma)), Variable(lfs_iname(ma.argexpr.element(), argnumber=ma.argexpr.number())))
@symbol
......
......@@ -141,7 +141,7 @@ def _lfs_iname(element, context):
return name
def lfs_iname(element, context=''):
def lfs_iname(element, context='', argnumber=None):
""" Get the iname to iterate over the local function space given by element
Arguments:
......@@ -155,6 +155,8 @@ def lfs_iname(element, context=''):
a given purpose, see the 'Loops and dependencies' of the loopy docs:
https://documen.tician.de/loopy/tutorial.html#loops-and-dependencies
"""
if argnumber is not None:
context = 'arg{}_{}'.format(argnumber, context)
return _lfs_iname(element, context)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment