gistfile1.txt

"""
This code does the following.

First, obtains an upwind DG1 approximation to div(u*D).
Then, tries to find a BDFM2 flux F such that div F is
equal to this upwind approximation.

we have 

\int_e phi D_1 dx = -\int_e grad phi . u D dx 
                    + \int_{\partial e} phi u.n \tilde{D} ds

where \tilde{D} is the value of D on the upwind face.

Then, if we define F such that

\int_f phi F.n ds = \int_f phi u.n \tilde{D} ds

then

\int_e phi div(F) ds = \int_{\partial e} phi u.n \tilde{D} ds
as required.
"""

from firedrake import *
from numpy import array, pi

mesh = IcosahedralSphereMesh(radius=1.0, refinement_level=5)
# Define global normal
global_normal = Expression(("x[0]/sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])",
"x[1]/sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])",
"x[2]/sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2])"))
mesh.init_cell_orientations(global_normal)
 
# Define function spaces and basis functions
V_dg = FunctionSpace(mesh, "DG", 1)
M = FunctionSpace(mesh, "RT", 1)
 
# advecting velocity
u0 = Expression(('-x[1]','x[0]','0'))
u = Function(M).project(u0)
 
# Mesh-related functions
n = FacetNormal(mesh)
 
# ( dot(v, n) + |dot(v, n)| )/2.0
un = 0.5*(dot(u, n) + abs(dot(u, n)))

dt = 1.0

# D advection equation
phi = TestFunction(V_dg)
D = TrialFunction(V_dg)
a_mass = phi*D*dx
a_int = dot(grad(phi), -u*D)*dx
a_flux = ( dot(jump(phi), un('+')*D('+') - un('-')*D('-')) )*dS

arhs = (a_int + a_flux)

D1 = Function(V_dg)

D0 = Expression("exp(-pow(x[2],2) - pow(x[1],2))")
D = Function(V_dg).interpolate(D0)

t = 0.0
T = 2*pi
k = 0
dumpfreq = 50

D1problem = LinearVariationalProblem(a_mass, action(arhs,D), D1)
D1solver = LinearVariationalSolver(D1problem, solver_parameters={'ksp_type': 'preonly',
                                                          'pc_type': 'lu'})
D1solver.solve()

# Surface Flux equation
RT_ele = FiniteElement("BDFM",triangle,degree=2)
Trace_space = FunctionSpace(mesh,TraceElement(RT_ele))
V1 = FunctionSpace(mesh,RT_ele)
Interior_ele = BrokenElement(FiniteElement("RT",triangle,degree=1))
Interior_space = FunctionSpace(mesh,Interior_ele)

W = MixedFunctionSpace((Trace_space,Interior_space))

gamma, w = TestFunctions(W)

Ft = TrialFunction(V1)
Fs = Function(V1)
 
aFs = (gamma('+')*inner(Ft('+'),n('+'))*dS
       + inner(w,Ft)*dx)
LFs = (gamma('+')*(un('+')*D('+')-un('-')*D('-'))*dS
       + inner(w,u*D)*dx)

Fsproblem = LinearVariationalProblem(aFs, LFs, Fs)
Fssolver = LinearVariationalSolver(Fsproblem)
Fssolver.solve()

file = File('div.pvd')
divFs = Function(V_dg)
solve(a_mass == phi*div(Fs)*dx, divFs)

L2err = assemble((D1-divFs)*(D1-divFs)*dx)
print L2err
assert(L2err<1.0e-12)