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