So I am very new to Fenics and Finite Element.
I am trying to solve a 1-D equation of the form:
h_t+q_x=0
q_t+(ghh0.5+qq/h)_x=0, with DG Finite Element method and this is my code below:
parameters["ghost_mode"] = "shared_facet"
mesh = UnitSquareMesh(8,8)
P1 = FiniteElement('DG', triangle, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
U=TrialFunction(V)
W=TestFunction(V)
h,q=split(U)
v1,v2=split(W)
U_old = Function(V)
h_old, q_old=split(U_old)
beta = Constant(0.5)
end_time=1
delta_t=0.1
#initial conditions
h0 = Constant(3.3)
q0=Expression("sin(2*pi*x[0])", degree=1, t=0) #
h_old=h0
q_old=q0
#Boundary Conditions
q_0 = Constant(0.0)
h_0 = Expression("sin(2*pi*(x[0] ))", degree=1)
bc_vector = as_vector([h_0, q_0])
bcs_h = DirichletBC(V.sub(0), h_0, "on_boundary", method='geometric')
bcs_q = DirichletBC(V.sub(1), q_0, "on_boundary", method='geometric')
bcs = [bcs_h, bcs_q]
y= Constant((1.0,1.0))
s0=Constant(1.0)
sf=Constant(0.0)
g = Constant(9.18)
time_steps = int(end_time/delta_t)
n = FacetNormal(mesh)
#fluxes
Fh_u = q
Fq_u=g*h*h*0.5+q*q/h
alpha_h=Constant(1.0)
alpha_q=ufl.Max(abs(q/h-sqrt(g*h)), abs(q/h+sqrt(g*h)))
flux_q = avg(dot(y*(g*h*h*0.5+q*q/h),n)) + 0.5*alpha_q*jump(q)
flux_h = avg(dot(y*q,n)) + 0.5*alpha_h*jump(h)
#Variational formulation
F_h= Constant(1/delta_t)*(h-h_old)*v1*dx+Fh_u*grad(v1)[0]*dx
F_h += dot(jump(v1), flux_h)*dS + dot(v1, h_0*n)*ds.
F_q=Constant(1/delta_t)*(q-q_old)*v2*dx+Fq_u*grad(v2)[0]*dx
F_q += dot(jump(v2), flux_q)*dS + dot(v2, q_0*n)*ds
F=F_h+F_q
Unew=Function(V)
t = 0
for n in range(time_steps):
solve(F==0, Unew, bcs)
U_old.assign(Unew)
h_old,q_new=split(Uold)
t += delta_t
my code is giving me some errors and I can not seem to figure out how to fix it.
-
ufl.log.UFLException: Discontinuous type Argument must be restricted. Error message if I run it
-
If I should change the n in the flux definitions to n(‘+’), I get
ufl.log.UFLException: Cannot restrict an expression twice. -
I get ufl.log.UFLException: Dot product requires non-scalar arguments, got arguments with ranks 0 and 1. for the last term in the variational formula
Any help on how to go about this will be much appreciated and FENICS tutorials on the shallow water equation as well. Thank you in advance.