# Solving Hyperbolic problem

**URL:** <https://fenicsproject.discourse.group/t/solving-hyperbolic-problem/6262>\
**Category:** Uncategorized\
**Created:** [July 21, 2021, 5:14am UTC](https://fenicsproject.discourse.group/t/solving-hyperbolic-problem/6262 "2021-07-21T05:14:32Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![RAVINA](https://avatars.discourse-cdn.com/v4/letter/r/c4cdca/32.png) [@RAVINA](https://fenicsproject.discourse.group/u/RAVINA)\
**Post date:** [July 21, 2021, 5:14am UTC](https://fenicsproject.discourse.group/t/solving-hyperbolic-problem/6262/1 "2021-07-21T05:14:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [July 21, 2021, 3:22pm UTC](https://fenicsproject.discourse.group/t/solving-hyperbolic-problem/6262/2 "2021-07-21T15:22:47Z")

</div>

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](https://comet-fenics.readthedocs.io/en/latest/demo/elastodynamics/demo_elastodynamics.py.html).

---

<div class="post-metadata">

**Author:** ![Powei](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/powei/32/4660_2.png) [@Powei](https://fenicsproject.discourse.group/u/Powei)\
**Post date:** [July 25, 2021, 6:36pm UTC](https://fenicsproject.discourse.group/t/solving-hyperbolic-problem/6262/3 "2021-07-25T18:36:41Z")

</div>

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

> <https://github.com/pwhuang/FE_in_Fluids_using_FEniCS/blob/master/FE_in_fluids_CH3.ipynb>

I hope this is helpful. Cheers!

---

<div class="post-metadata">

**Author:** ![RAVINA](https://avatars.discourse-cdn.com/v4/letter/r/c4cdca/32.png) [@RAVINA](https://fenicsproject.discourse.group/u/RAVINA)\
**Post date:** [July 26, 2021, 6:35am UTC](https://fenicsproject.discourse.group/t/solving-hyperbolic-problem/6262/4 "2021-07-26T06:35:12Z")

</div>

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.

```auto
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.

```auto
---------------------------------------------------------------------------
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.

---

<div class="post-metadata">

**Author:** ![RAVINA](https://avatars.discourse-cdn.com/v4/letter/r/c4cdca/32.png) [@RAVINA](https://fenicsproject.discourse.group/u/RAVINA)\
**Post date:** [July 26, 2021, 6:36am UTC](https://fenicsproject.discourse.group/t/solving-hyperbolic-problem/6262/5 "2021-07-26T06:36:01Z")

</div>

Thank you for reply.

---

<div class="post-metadata">

**Author:** ![RAVINA](https://avatars.discourse-cdn.com/v4/letter/r/c4cdca/32.png) [@RAVINA](https://fenicsproject.discourse.group/u/RAVINA)\
**Post date:** [July 27, 2021, 5:14am UTC](https://fenicsproject.discourse.group/t/solving-hyperbolic-problem/6262/6 "2021-07-27T05:14:16Z")

</div>

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

Thanks in advance.

---

<div class="post-metadata">

**Author:** ![Yongqiang\_Chen](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/yongqiang_chen/32/7620_2.png) [@Yongqiang\_Chen](https://fenicsproject.discourse.group/u/Yongqiang_Chen)\
**Post date:** [October 3, 2024, 2:38am UTC](https://fenicsproject.discourse.group/t/solving-hyperbolic-problem/6262/7 "2024-10-03T02:38:45Z")

</div>

Hi Nate,

Happy to hear that Fenics can solve any PDEs. Can you give some suggestions to solve the following equations? It has two variables in each equation.  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/56f94906ff0e76fcf439671cb3818f7739b2ef77.png)

THanks,  
Yongqiang