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.
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.
THanks,
Yongqiang