I figured out from http://hplgit.github.io/INF5620/doc/pub/fenics_tutorial1.1/tu2.html in the section avoiding assembly, that I could ease the load on the form compiler by precompiling the forms for the Neumann boundary conditions.
I have zero flux Neumann BCs at the limits of a 1D interval and an image charge potential due to the charged species in the domain. The image charges are allocated at each boundary according to the centroid of the concentration profile in the domain. I have an initial concentration at a finite thickness close to the left boundary.
Since I am solving a Poisson-Nernst-Planck problem in a mixed space,
\frac{\partial c}{\partial t}
=
\nabla \cdot \left(D \nabla c \right)
+ \mu \left(\nabla \cdot \nabla \phi\right) c
+ \mu \left(\nabla \phi \right) \cdot\nabla c\\
0 = \nabla \cdot \nabla \phi + \frac{ze c }{\epsilon}.
If I just expand the derivatives, I end up with a term
\int_{\Omega} u_c \nabla u_p \cdot \nabla v_c dx,
where u_c, u_p are trial functions for c and \phi respectively.
I am currently using
up_n = Function(W)
P = assemble(inner(uc*grad(up_n),grad(vc))*dx
Which works in the sense that I get a problem that I can solve using gmres solver directly. The result deviates from what I get if I let fenics compile the forms automatically, because the problem is non-linear. So, I need to solve the problem using Newton’s method.
But this does not clarify how to assemble the P term. In the linear algebra level, the crossed term should look something like
\int_{\Omega}\Sigma_i U_{c,i} \phi_{c,i} \nabla\Sigma_j U_{p,j} \phi_{p,j} \hat{\phi}_{c,k} dx = \int_{\Omega}\Sigma_{i,j} U_{c,i} U_{p,j} \nabla \phi_{p,j} \phi_{c,i} \hat{\phi}_{c,k} dx
Whenever I try to assemble the problem as:
CG1 = FiniteElement('CG',mesh.ufl_cell(),1)
W = FunctionSpace(mesh,MixedElement([CG1,CG1])
u = TrialFunction(W)
uc,up = split(u)
v = TestFunction(W)
vc,vp = split(v)
def getPForm(u_c,u_p):
PP = u_c*inner(grad(u_p),grad(v_c))*dx
return assemble(PP)
P = getPForm(uc,up)
Gives me the following stack trace:
File "./singlelayer.py", line 446, in main
P = getPForm(c,phi)
File "./singlelayer.py", line 337, in getPForm
return assemble(PP)
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/assembling.py", line 198, in assemble
dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/assembling.py", line 58, in _create_dolfin_form
function_spaces=function_spaces)
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/form.py", line 44, in __init__
mpi_comm=mesh.mpi_comm())
File "/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py", line 47, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py", line 97, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File "/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py", line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File "/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py", line 133, in jit_build
generate=jit_generate)
File "/usr/local/lib/python3.6/dist-packages/dijitso/jit.py", line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params["generator"])
File "/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py", line 66, in jit_generate
prefix=module_name, parameters=parameters, jit=True)
File "/usr/local/lib/python3.6/dist-packages/ffc/compiler.py", line 143, in compile_form
prefix, parameters, jit)
File "/usr/local/lib/python3.6/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
File "/usr/local/lib/python3.6/dist-packages/ffc/analysis.py", line 90, in analyze_ufl_objects
for form in forms)
File "/usr/local/lib/python3.6/dist-packages/ffc/analysis.py", line 90, in <genexpr>
for form in forms)
File "/usr/local/lib/python3.6/dist-packages/ffc/analysis.py", line 174, in _analyze_form
do_apply_restrictions=True)
File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/compute_form_data.py", line 418, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode) # Currently testing how fast this is
File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
check_integrand_arity(itg.integrand(), arguments, complex_mode)
File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
arg_tuples = map_expr_dag(rules, expr, compress=False)
File "/usr/local/lib/python3.6/dist-packages/ufl/corealg/map_dag.py", line 37, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File "/usr/local/lib/python3.6/dist-packages/ufl/corealg/map_dag.py", line 86, in map_expr_dags
r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py", line 63, in product
raise ArityMismatch("Multiplying expressions with overlapping form argument number {0}, argument is {1}.".format(x[0].number(), _afmt(x)))
File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py", line 20, in _afmt
for arg, conj in atuple)
File "/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py", line 20, in <genexpr>
for arg, conj in atuple)
TypeError: 'bool' object is not iterable
Any help is appreciated!