Hello,
I am trying to solve a 1D problem with a boundary condition such the u>=0 always.
The sample problem I am trying to solve is
\begin{equation}
au_x-\nu u_{xx}= s \, \text{in} \, [0,1]
\\
\text{with } s=5e^{-100(x-\frac{1}{8})^2}-5e^{-100(x-\frac{1}{4})^2}
\\
u=0 \, \text{at} \, x=0 \,\text{and}\,x=1
\end{equation}
For the given case a=1 ,\nu=0.01
The solution needs an SUPG scheme, and the solution is
In the next step i need to enforce a penalty on this problem such that u>=0 on the domain.
I am trying to solve thos by adding a penalty term to system such that
\begin{equation}
au_x-\nu u_{xx} - \zeta*min(u,0)=s
\end{equation}
where \zeta is a large number say 10^5.
The corresponding weak form term will be
\int_{\Omega} \zeta*min(u,0)*wd\Omega
However I have tried my best to include this in my code, but nothing seems to work. I am unable to define a function that can be parsed by fenics. Could someone help me on how to define this penalty term?
My code
from __future__ import print_function
from ufl import nabla_div
from ufl import *
from fenics import *
#import pylab as plt
import numpy as np
import matplotlib.pyplot as plt
a = 1.0 # vel coef (1d a vector = scalar)
nu = 0.01 # diff. coef => Pe = 5.0
#Step 1 :Define the unit interval domain
mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, "CG", 1)
#Step 2: Define the BC
def u0_boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, 0.0, u0_boundary)
#Step 3: Define trial and test function
u = TrialFunction(V)
w = TestFunction(V)
# Step 4: Define the SUPG stabilization tau
#h = mesh.hmin() # uniform mesh
h = CellDiameter(mesh)
Pe = a*h/2.0/nu;
beta = 1.0/tanh(Pe) - 1/Pe
nu2 = beta*a*h/2.0;
tau = nu2/(a*a)
# variational forms
s = Expression('5*exp(-100*(x[0]-1./8)*(x[0]-1./8)) - 5*exp(-100*(x[0]-1./4)*(x[0]-1./4))',degree=10)
a_ = w*dot(a, u.dx(0))*dx + inner(w.dx(0), nu*u.dx(0))*dx
L_ = w*s*dx
# stabilization term
Lu_ = dot(a, u.dx(0)) # linear elem => div(nabla_grad(u)) == 0
Pw_ = dot(a, w.dx(0)) # stabal. operator
# stabilized forms
a_ += tau*Pw_*Lu_*dx
L_ += tau*Pw_*s*dx
#############################################################
#define penalty function
def uneg(x):
#try3
val=conditional(ge(x,0),0.0,x)
#try 2
#if x>=0.0:
# val=0.0
#else:
# val=x
#try1
#return (x-abs(x))/2.0 # returns 0 if u is +ve, u if u is negative
return val # returns 0 if u is +ve, u if u is negative
#uneg(-1)
#uneg(1)
# add penalty constraint
penalty=1.0e5
a_-= w*penalty*uneg(u)*dx
#############################################################################
# solve
u = Function(V)
solve(a_ == L_, u, bc)
fig = plt.figure ( )
ax=plt.subplot(111)
plot(u, label='SUPG')
ax.legend ( )
ax.grid ( True )
The output I get is
Calling FFC just-in-time (JIT) compiler, this may take some time.
---------------------------------------------------------------------------
ArityMismatch Traceback (most recent call last)
<ipython-input-86-ce3ff92c8e8e> in <module>
71 # solve
72 u = Function(V)
---> 73 solve(a_ == L_, u, bc)
74 fig = plt.figure ( )
75 ax=plt.subplot(111)
/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py in solve(*args, **kwargs)
218 # tolerance)
219 elif isinstance(args[0], ufl.classes.Equation):
--> 220 _solve_varproblem(*args, **kwargs)
221
222 # Default case, just call the wrapped C++ solve function
/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py in _solve_varproblem(*args, **kwargs)
240 # Create problem
241 problem = LinearVariationalProblem(eq.lhs, eq.rhs, u, bcs,
--> 242 form_compiler_parameters=form_compiler_parameters)
243
244 # Create solver and call solve
/usr/local/lib/python3.6/dist-packages/dolfin/fem/problem.py in __init__(self, a, L, u, bcs, form_compiler_parameters)
54 else:
55 L = Form(L, form_compiler_parameters=form_compiler_parameters)
---> 56 a = Form(a, form_compiler_parameters=form_compiler_parameters)
57
58 # Initialize C++ base class
/usr/local/lib/python3.6/dist-packages/dolfin/fem/form.py in __init__(self, form, **kwargs)
42
43 ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
---> 44 mpi_comm=mesh.mpi_comm())
45 ufc_form = cpp.fem.make_ufc_form(ufc_form[0])
46
/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py in mpi_jit(*args, **kwargs)
45 # Just call JIT compiler when running in serial
46 if MPI.size(mpi_comm) == 1:
---> 47 return local_jit(*args, **kwargs)
48
49 # Default status (0 == ok, 1 == fail)
/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py in ffc_jit(ufl_form, form_compiler_parameters)
95 p.update(dict(parameters["form_compiler"]))
96 p.update(form_compiler_parameters or {})
---> 97 return ffc.jit(ufl_form, parameters=p)
98
99
/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py in jit(ufl_object, parameters, indirect)
215
216 # Inspect cache and generate+build if necessary
--> 217 module = jit_build(ufl_object, module_name, parameters)
218
219 # Raise exception on failure to build or import module
/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py in jit_build(ufl_object, module_name, parameters)
131 name=module_name,
132 params=params,
--> 133 generate=jit_generate)
134 return module
135
/usr/local/lib/python3.6/dist-packages/dijitso/jit.py in jit(jitable, name, params, generate, send, receive, wait)
163 elif generate:
164 # 1) Generate source code
--> 165 header, source, dependencies = generate(jitable, name, signature, params["generator"])
166 # Ensure we got unicode from generate
167 header = as_unicode(header)
/usr/local/lib/python3.6/dist-packages/ffc/jitcompiler.py in jit_generate(ufl_object, module_name, signature, parameters)
64
65 code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
---> 66 prefix=module_name, parameters=parameters, jit=True)
67
68 # Jit compile dependent objects separately,
/usr/local/lib/python3.6/dist-packages/ffc/compiler.py in compile_form(forms, object_names, prefix, parameters, jit)
141 """This function generates UFC code for a given UFL form or list of UFL forms."""
142 return compile_ufl_objects(forms, "form", object_names,
--> 143 prefix, parameters, jit)
144
145
/usr/local/lib/python3.6/dist-packages/ffc/compiler.py in compile_ufl_objects(ufl_objects, kind, object_names, prefix, parameters, jit)
183 # Stage 1: analysis
184 cpu_time = time()
--> 185 analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
186 _print_timing(1, time() - cpu_time)
187
/usr/local/lib/python3.6/dist-packages/ffc/analysis.py in analyze_ufl_objects(ufl_objects, kind, parameters)
88 # Analyze forms
89 form_datas = tuple(_analyze_form(form, parameters)
---> 90 for form in forms)
91
92 # Extract unique elements accross all forms
/usr/local/lib/python3.6/dist-packages/ffc/analysis.py in <genexpr>(.0)
88 # Analyze forms
89 form_datas = tuple(_analyze_form(form, parameters)
---> 90 for form in forms)
91
92 # Extract unique elements accross all forms
/usr/local/lib/python3.6/dist-packages/ffc/analysis.py in _analyze_form(form, parameters)
172 do_apply_geometry_lowering=True,
173 preserve_geometry_types=(Jacobian,),
--> 174 do_apply_restrictions=True)
175 elif r == "tsfc":
176 try:
/usr/local/lib/python3.6/dist-packages/ufl/algorithms/compute_form_data.py in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
416 preprocessed_form = remove_complex_nodes(preprocessed_form)
417
--> 418 check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode) # Currently testing how fast this is
419
420 # TODO: This member is used by unit tests, change the tests to
/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py in check_form_arity(form, arguments, complex_mode)
175 def check_form_arity(form, arguments, complex_mode=False):
176 for itg in form.integrals():
--> 177 check_integrand_arity(itg.integrand(), arguments, complex_mode)
/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py in check_integrand_arity(expr, arguments, complex_mode)
157 key=lambda x: (x.number(), x.part())))
158 rules = ArityChecker(arguments)
--> 159 arg_tuples = map_expr_dag(rules, expr, compress=False)
160 args = tuple(a[0] for a in arg_tuples)
161 if args != arguments:
/usr/local/lib/python3.6/dist-packages/ufl/corealg/map_dag.py in map_expr_dag(function, expression, compress)
35 Return the result of the final function call.
36 """
---> 37 result, = map_expr_dags(function, [expression], compress=compress)
38 return result
39
/usr/local/lib/python3.6/dist-packages/ufl/corealg/map_dag.py in map_expr_dags(function, expressions, compress)
82 # Cache miss: Get transformed operands, then apply transformation
83 if cutoff_types[v._ufl_typecode_]:
---> 84 r = handlers[v._ufl_typecode_](v)
85 else:
86 r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
/usr/local/lib/python3.6/dist-packages/ufl/algorithms/check_arities.py in nonlinear_operator(self, o)
39 for t in traverse_unique_terminals(o):
40 if t._ufl_typecode_ == Argument._ufl_typecode_:
---> 41 raise ArityMismatch("Applying nonlinear operator {0} to expression depending on form argument {1}.".format(o._ufl_class_.__name__, t))
42 return self._et
43
ArityMismatch: Applying nonlinear operator GE to expression depending on form argument v_1.