diff --git a/pydemo/scalar/scalar.py b/pydemo/scalar/scalar.py new file mode 100644 index 0000000000000000000000000000000000000000..46f60bab979b27e5f44f0f9e35a019bd07b4f0bb --- /dev/null +++ b/pydemo/scalar/scalar.py @@ -0,0 +1,27 @@ +from ufl import * +from dune.ufl import Space +class Transport1D: + name = 'transport' + + def F_c(U): + return [ [U[0], 0] ] + + def outflowValue(u): + return u + boundaryValue = {1: outflowValue} + + def maxLambda(U,n): + return abs(n[0]) + def velocity(U): + return as_vector( [1,0] ) + def physical(U): + return 1 + def jump(U,V): + return U[0]-V[0] + +space = Space(2,1) +def RiemanProblem(UL,UR): + x = SpatialCoordinate(space.cell()) + return conditional(x[0]<-0.5,as_vector(UL),as_vector(UR)) +def Constant(): + return as_vector( [0.1] ) diff --git a/pydemo/scalar/shock_tube.py b/pydemo/scalar/shock_tube.py new file mode 100644 index 0000000000000000000000000000000000000000..5e5bafcfe1e43fa7906c92864d324fcd3637aceb --- /dev/null +++ b/pydemo/scalar/shock_tube.py @@ -0,0 +1,50 @@ +import time +from dune.grid import structuredGrid, cartesianDomain +from dune.fem import parameter +import dune.create as create +from dune.models.elliptic.formfiles import loadModels +from dune.femdg import createFemDGSolver + +from scalar import RiemanProblem +from scalar import Transport1D as Model + +parameter.append("parameter") +parameter.append({"fem.verboserank": -1}) + +x0,x1,N = [-1, 0], [1, 0.1], [20, 5] +grid = structuredGrid(x0,x1,N) +# grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N)) +dimR = 1 +t = 0 +saveStep = 0.01 +saveTime = saveStep +count = 0 +endTime = 1. + +def useODESolver(): + global count, t, saveTime + initial = RiemanProblem( [1], [0] ) + + spaceName = "dgonb" + polOrder = 2 + space = create.space(spaceName, grid, order=polOrder, dimrange=dimR) + u_h = space.interpolate(initial, name='u_h') + u_h_n = u_h.copy(name="previous") + operator = createFemDGSolver( Model, space ) + + start = time.time() + grid.writeVTK(Model.name, celldata=[u_h], number=count) + while t < endTime: + u_h_n.assign(u_h) + operator(u_h_n, u_h) + dt = operator.deltaT() + t += dt + if t > saveTime: + print('dt = ', dt, 'time = ',t, 'count = ',count ) + count += 1 + grid.writeVTK(Model.name, celldata=[u_h], number=count) + saveTime += saveStep + grid.writeVTK(Model.name, celldata=[u_h], number=count) + print("time loop:",time.time()-start) + +useODESolver()