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"))