# Solving Hyperbolic problem

Hello everyone, I want to solve a hyperbolic pde . Can we solve it in fenics? If yes then kindly suggest me some example.

Thank you.

FEniCS is a general toolbox for the approximation of PDE solutions by the finite element method. Anything you can dream up is possible, it just hasnâ€™t been implemented yet.

As an example consider this excellent elastodynamics problem.

I have some examples of solving 1D advection equation, which is an example of a hyperbolic PDE. See:

I hope this is helpful. Cheers!

First of all thanku for the reply.
I am trying to solve the following problem:
(u_{tt},v) + (\nabla{u}, \nabla{v}) + (\alpha u_{t},v) + (\nu u,v) + (f(u),v) \forall v\in V
I am using G-\alpha method.

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

# Form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

mesh = UnitSquareMesh(8,8)
T = 1
num_steps = 20
dt = T / num_steps
V = FunctionSpace(mesh,'CR',1)

v = TestFunction(V)
du = TrialFunction(V)
u = Function(V) # unknown displacement

def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, Constant(0.0) , boundary)

u_D = Expression('sin(pi*x[0])*sin(pi*x[1])*cos(sqrt(2)*pi*t)',degree=2,t=0)
u_t = Expression('-sqrt(2)*pi*sin(pi*x[0])*sin(pi*x[1])*sin(sqrt(2)*pi*t)',degree=2,t=0)
u_old = interpolate(u_D,V)
ut_old = interpolate(u_t,V)
utt_old = Function(V)

alpha_m = Constant(0.5)
alpha_f = Constant(0.5)
gamma   = Constant(0.5)
beta    = Constant(0.25)
alpha = Constant(0.2)
nu = Constant(0.5)

def m(u, v):
return inner(u, v)*dx

def k(u, v):
return inner(grad(u), grad(v))*dx

def c(u, v):
return inner(alpha*u, v)*dx

def d(u,v):
return inner(nu*u, v)*dx

def f(u,v):
return inner(u*u*u, v)*dx

#Update formula for acceleration
def update_utt(u, u_old, ut_old, utt_old, ufl=True):
if ufl:
dt_ = dt
beta_ = beta
else:
dt_ = float(dt)
beta_ = float(beta)
return (u-u_old-dt_*ut_old)/(beta_*dt_**2) - ((1-2*beta_)/2*beta_)*utt_old

# Update formula for velocity
def update_ut(utt, u_old, ut_old, utt_old, ufl=True):
if ufl:
dt_ = dt
gamma_ = gamma
else:
dt_ = float(dt)
gamma_ = float(gamma)
return ut_old + dt_*((1-gamma_)*utt_old + gamma_*utt)

def update_fields(u, u_old, ut_old, utt_old):
"""Update fields at the end of each time step."""

# Get vectors
u_vec, u0_vec  = u.vector(), u_old.vector()
ut0_vec, utt0_vec = ut_old.vector(), utt_old.vector()

# use update functions using vector arguments
utt_vec = update_utt(u_vec, u0_vec, ut0_vec, utt0_vec, ufl=False)
ut_vec = update_ut(utt_vec, u0_vec, ut0_vec, utt0_vec, ufl=False)

# Update (u_old <- u)
ut_old.vector()[:], utt_old.vector()[:] = ut_vec, utt_vec
u_old.vector()[:] = u.vector()

def avg(x_old, x_new, alpha):
return alpha*x_old + (1-alpha)*x_new

utt_new = update_utt(du, u_old, ut_old, utt_old, ufl=True)
ut_new = update_ut(utt_new, u_old, ut_old, utt_old, ufl=True)
F = m(avg(utt_old, utt_new, alpha_m), v) + c(avg(ut_old, ut_new, alpha_f), v) \
+ k(avg(u_old, du, alpha_f), v)  + d(avg(u_old, du, alpha_f), v)+ f(avg(u_old, du, alpha_f),v)

time = np.linspace(0, T, num_steps+1)
for (i, dt) in enumerate(np.diff(time)):

t = time[i+1]
print("Time: ", t)
u_D.t = t
u_t.t = t
solve(F==0,u,bc)

# Update old fields with new quantities
update_fields(u, u_old, ut_old, utt_old)


I am getting the following error.

---------------------------------------------------------------------------
ArityMismatch                             Traceback (most recent call last)
<ipython-input-51-e1482e75ceab> in <module>
6     u_D.t = t
7     u_t.t = t
----> 8     solve(F==0,u,bc)
9
10

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py in solve(*args, **kwargs)
231     # tolerance)
232     elif isinstance(args[0], ufl.classes.Equation):
--> 233         _solve_varproblem(*args, **kwargs)
234
235     # Default case, just call the wrapped C++ solve function

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py in _solve_varproblem(*args, **kwargs)
307             # Create problem
308             problem = NonlinearVariationalProblem(eq.lhs, u, bcs, J,
--> 309                                                   form_compiler_parameters=form_compiler_parameters)
310
311             # Create solver and call solve

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/problem.py in __init__(self, F, u, bcs, J, form_compiler_parameters)
155
156         # Wrap forms
--> 157         F = Form(F, form_compiler_parameters=form_compiler_parameters)
158         if J is not None:
159             J = Form(J, form_compiler_parameters=form_compiler_parameters)

/usr/lib/petsc/lib/python3/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/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py in mpi_jit(*args, **kwargs)
48         # Just call JIT compiler when running in serial
49         if MPI.size(mpi_comm) == 1:
---> 50             return local_jit(*args, **kwargs)
51
52         # Default status (0 == ok, 1 == fail)

/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py in ffc_jit(ufl_form, form_compiler_parameters)
98     p.update(dict(parameters["form_compiler"]))
99     p.update(form_compiler_parameters or {})
--> 100     return ffc.jit(ufl_form, parameters=p)
101
102

/usr/lib/python3/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/lib/python3/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/lib/python3/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/lib/python3/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/lib/python3/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/lib/python3/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/lib/python3/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/lib/python3/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/lib/python3/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/lib/python3/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)
405     preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)
406
--> 407     check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
408
409     # TODO: This member is used by unit tests, change the tests to

/usr/lib/python3/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/lib/python3/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/lib/python3/dist-packages/ufl/corealg/map_dag.py in map_expr_dag(function, expression, compress)
24     Return the result of the final function call.
25     """
---> 26     result, = map_expr_dags(function, [expression], compress=compress)
27     return result
28

/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py in map_expr_dags(function, expressions, compress)
73                 r = handlers[v._ufl_typecode_](v)
74             else:
---> 75                 r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
76
77             # Optionally check if r is in rcache, a memory optimization

/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py in sum(self, o, a, b)
46     def sum(self, o, a, b):
47         if a != b:
---> 48             raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
49         return a
50

ArityMismatch: Adding expressions with non-matching form arguments ('v_1',) vs ().


Problem arising due to non linearity in f. Kindly help me to resolve it.

Thank you in advance.

Thank you for reply.

For linear problem, its working. I am facing problem due to non-linear term. Kindly help me to resolve it.

Thanks in advance.