Hi Onno, If you check out the QG demo, it uses an RK scheme with DG for pv. I've been intending to supply a demo for DG transient advection fora while, maybe over the summer I'll find time. Cheers cjc On Friday, 17 June 2016, Onno Bokhove <O.Bokhove@leeds.ac.uk> wrote:
Hi FD's
Trying to get my feet wet. Simple query.
Trying to make a DG conservation law/system code,
starting with linear advection, backward Euler in steps.
How do I make the code snippet below forward Euler
or Crank-Nicolson? I used all the old time step terms in L
but that does not work. (Now commented out -thet factor terms)
Also, why does the nonlinear solver not work?
I think it is because the printing can't be done in DG0
so there needs to assign between the DG solution of the
prime unknown.
Would be nice to have a standard DG example for
U_t + div F = 0
forward Euler or RK3.
Thanks for the reply as then I can keep playing in the plane!
Onno
P.S. Code below
tend = 4.211
timestep = tend/40
V = FunctionSpace(tmesh, "DG", 0)
W = VectorFunctionSpace(tmesh, "CG", 1)
# Advection velocity
velocity = Expression(("1.0", "0.0"))
u0 = project(velocity, W)
#Inflow BC
inflow = Expression("(x[0] < 0.0) ? 1.0 : 0")
inflow2 = Expression("(x[0] < 0.5) ? 1.0 : 0")
inflow3 = Expression("(x[0] < 1.0) ? 1.0 : 0")
#Inflow Variable
D0 = Function(V)
D0.interpolate(inflow)
m1 = Function(V)
m1.interpolate(inflow2)
m2 = Function(V)
m2.interpolate(inflow3)
#Initial Condition
ic = Expression("(x[0] < .5) ? 1.0 : 0")
ic2 = Expression("(x[0] < .5) ? 1.0 : 0")
ic3 = Expression("(x[0] < .5) ? 1.0 : 0")
#Set up element normal
n = FacetNormal(tmesh)
#
un = 0.5*(dot(u0, n) + abs(dot(u0, n)))
# Previous timestep for BE timestep
D_ = Function(V)
m1_ = Function(V)
m2_ = Function(V)
D_.interpolate(ic)
m1_.interpolate(ic2)
m2_.interpolate(ic3)
D = TrialFunction(V)
phi = TestFunction(V)
#
# Setting up Bilinear form
#
thet = 0;
a1 = ((D/timestep)*phi - (1-thet)*D*dot(u0, grad(phi)))*dx # (1-thet) interior term at new time step
a2 = (1-thet)*dot(jump(phi), un('+')*D('+') - un('-')*D('-'))*dS # Internal flux
a3 = (1-thet)*dot(phi, un*D)*(ds(112)) # outflow at top wall
a = a1 + a2 + a3
#
L1 = (D_/timestep)*phi*dx - D0*phi*dot(u0, n)*(ds(114)) # previous timestep & given inflow at bottom wall
# L2 = thet*(D_*dot(u0, grad(phi)))*dx # thet times interior term at old time stepL
# L3 = -thet*dot(jump(phi), un('+')*D_('+') - un('-')*D_('-'))*dS # Internal flux
# L4 = -thet*dot(phi, un*D_)*(ds(112)) # outflow at top wall
L = L1
F = a-L1
out = Function(V)
outfile = File("./Results/swe2d.pvd")
# Write IC to file
# outfile << D_
D_.assign(out)
outfile.write(out)
# outfile.wrote(project(out, W, name="D"))
t = 0.0
while (t <= tend):
#Solve problem
solve(a == L, out)
#solve(F == 0, out)
#Update previous timestep
D_.assign(out)
#Update time
t+=timestep
# outfile << out
outfile.write(out)
# outfile.write(project(out, W, name="D"))
-- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916