Hello everyone, I am using dolfin_dg to perform some simulations of modified euler equations. basically what I want to solve is the euler equations with a modified pressure and a polytropic relation substituting the energy equation.
\frac{d \vec q_j}{dt} + \vec\nabla_i \textit{F}_{i,j} = 0
where
q_j = (n,n m u_1, n m u_2)^T\\
F_{i,j}=(n\vec u, m n \vec( u \otimes\vec u)+ T n I_{2 \times 2})^T
with T a constant.
I tried to supply the flux as daughter class inheriting from the HyperbolicOperator class:
class PlasmaOperator(HyperbolicOperator):
def __init__(self, mesh, V, bcs):
dim = mesh.geometric_dimension()
def F_c(U):
n_i = U[0]
n_m_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_m_u/(m_i*n_i)
p = n_i*electron_temp
inertia = n_i*ufl.outer(u, u) + p*Identity(dim)
res = ufl.as_tensor([n_i*u,
*[inertia[d, :] for d in range(dim)]])
return res
def alpha(U, n):
n_i = U[0]
n_m_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_m_u/(m_i*n_i+DOLFIN_EPS)
c = sqrt(T/m_i)
lambdas = [dot(u, n) - c, dot(u, n), dot(u, n) + c]
return lambdas
HyperbolicOperator.__init__(self, mesh, V, bcs, F_c, LocalLaxFriedrichs(alpha))
The boundary conditions I am trying to implement are basically an Dirichlet in the left side of the box and free outflow anywhere, while the initial condition should be zero density and velocity everywhere, this generates problems as you get 0/0, anyway the zero density is not a strong constraint as I could change it. Here’s how I implement them:
gDThroat = as_vector((n0, m_i*n0*u0, m_i*n0*v0))
class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] < DOLFIN_EPS
class Outlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] > DOLFIN_EPS
bound_mshfunc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) #Facets
bound_mshfunc.set_all(0)
INLET, OUTFLOW = 1, 2
inlets = Inlet()
inlets.mark(bound_mshfunc, INLET)
outlet = Outlet()
outlet.mark(bound_mshfunc, OUTFLOW)
ds = Measure('ds', domain=mesh, subdomain_data=bound_mshfunc)
bcs = [DGDirichletBC(ds(INLET), gDThroat),
DGDirichletBC(ds(OUTFLOW),gD2)]
However when I solve the problem the residual from Newton’s solver increases steadily and diverges.
Here is my full code:
#%%
import matplotlib.pyplot as plt
import ufl
from dolfin import *
from dolfin_dg import *
parameters["ghost_mode"] = "shared_facet"
# Mesh and function space.
N_ele = 10
p = 2
mesh = UnitSquareMesh(N_ele,N_ele,'left/right')
V = VectorFunctionSpace(mesh, 'DG', p, dim=3)
du = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)
#Constants
electron_temp = 20 #T
m_i = 1
e = 1.6
n0 = Constant(1)
u0 = Constant(1)
v0 = Constant(0.0)
gDThroat = as_vector((n0, m_i*n0*u0, m_i*n0*v0))
gD2 = as_vector((0.1, 0, 0))
u_vec_0 = Constant((n0, m_i*n0*u0, m_i*n0*v0))
u_vec = interpolate(u_vec_0,V)
T = 0.2
h = 1/N_ele
CFL = 0.6
dt_n = CFL*h/(p+1)**2
N = int(T/dt_n)
dt = Constant(dt_n)
un = Function(V)
###################################################
class PlasmaOperator(HyperbolicOperator):
def __init__(self, mesh, V, bcs):
dim = mesh.geometric_dimension()
def F_c(U):
n_i = U[0]
n_m_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_m_u/(m_i*n_i)
p = n_i*electron_temp
inertia = n_i*ufl.outer(u, u) + p*Identity(dim)
res = ufl.as_tensor([n_i*u,
*[inertia[d, :] for d in range(dim)]])
return res
def alpha(U, n):
n_i = U[0]
n_m_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_m_u/(m_i*n_i+DOLFIN_EPS)
c = sqrt(electron_temp/m_i)
lambdas = [dot(u, n) - c, dot(u, n), dot(u, n) + c]
return lambdas
HyperbolicOperator.__init__(self, mesh, V, bcs, F_c, LocalLaxFriedrichs(alpha))
###################################################
class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] < DOLFIN_EPS
class Outlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] > DOLFIN_EPS
bound_mshfunc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) #Facets
bound_mshfunc.set_all(0)
INLET, OUTFLOW = 1, 2
inlets = Inlet()
inlets.mark(bound_mshfunc, INLET)
outlet = Outlet()
outlet.mark(bound_mshfunc, OUTFLOW)
ds = Measure('ds', domain=mesh, subdomain_data=bound_mshfunc)
bcs = [DGDirichletBC(ds(INLET), gDThroat),
DGDirichletBC(ds(OUTFLOW),gD2)]
po = PlasmaOperator(mesh, V, bcs= bcs)
residual = po.generate_fem_formulation(u_vec, v)
residual += dot(u_vec - un, v)*dx
J = derivative(residual, u_vec, du)
class Problem(NonlinearProblem):
def F(self, b, x):
assemble(residual, tensor=b)
def J(self, A, x):
assemble(J, tensor=A)
solver = NewtonSolver()
for j in range(N):
un.assign(u_vec)
solver.solve(Problem(), u_vec.vector())
Thanks in advance.
