Hello everyone,
I am trying to solve an euler-like problem with dolfin_dg.
The problem converges for order p=0 (density at final time depicted here):
However when I try solving for higher order polynomials the problem does not converge:
Newton iteration 0: r (abs) = 9.062e-02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 2.960e-02 (tol = 1.000e-10) r (rel) = 3.266e-01 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 5.203e-02 (tol = 1.000e-10) r (rel) = 5.742e-01 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 9.504e+00 (tol = 1.000e-10) r (rel) = 1.049e+02 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 3.095e+02 (tol = 1.000e-10) r (rel) = 3.415e+03 (tol = 1.000e-09)
Newton iteration 5: r (abs) = 1.546e+06 (tol = 1.000e-10) r (rel) = 1.706e+07 (tol = 1.000e-09)
Newton iteration 6: r (abs) = 4.473e+11 (tol = 1.000e-10) r (rel) = 4.935e+12 (tol = 1.000e-09)
I tried checking all the boundaries where marked correctly:
As well as using a finer mesh and different timesteps:
But none of these attempts lead to convergence.
Here is my code:
import ufl
from dolfin import *
from dolfin_dg import *
import matplotlib.pyplot as plt
import magpylib as mpl
import meshio
parameters["ghost_mode"] = "shared_facet"
parameters['allow_extrapolation'] = True
# Mesh and function space.
# Dimensions scaled by the radius of the coils
x1, x2 = 0, 15
y1, y2 = 0, 15
N_ele = 20
p = 0
mesh = RectangleMesh(Point(x1,y1), Point(x2,y2),2*N_ele,N_ele,'left/right')
V = VectorFunctionSpace(mesh, 'DG', p, dim=3) #Vector function space for the variables
V1 = FunctionSpace(mesh, 'DG', p) # '' for scalars
du = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)
#Constants of the problem
cs = sqrt(1)
R = 1
# Set up the boundary conditions
n0 = Expression('x[0]> -R and x[0]< R ? cos(0.5*pi*x[0]) : 0', R=R, element = V1.ufl_element())
ux0 = Constant(0)
uz0 = Constant(2)
gDThroat = as_vector((n0, n0*ux0, n0*uz0))
gOut = as_vector((0, 0, 0))
# And intitial conditions
u_vec_0 = Expression(('0.0000001*n0','0.0', '0.0'), n0=n0, element= V.ufl_element())
u_vec = interpolate(u_vec_0,V)
# time step
T = 8
dt_n = 0.005*T
N = int(T/dt_n)
dt = Constant(dt_n)
un = Function(V)
###################################################
#Define Boundaries and boundary conditions
class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[1]<DOLFIN_EPS and between(x[0],(0,1))
bound_mshfunc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) #Facets
bound_mshfunc.set_all(0)
INLET, OUTFLOW = 1, 2
# Mark the whole exterior as outflow
AutoSubDomain(lambda x, on_boundary: on_boundary).mark(bound_mshfunc, OUTFLOW)
# Then overwrite the inflow and outflow markers
inlets = Inlet()
inlets.mark(bound_mshfunc, INLET)
ds = Measure('ds', domain=mesh, subdomain_data=bound_mshfunc)
bcs = [DGDirichletBC(ds(INLET), gDThroat),
DGNeumannBC(ds(OUTFLOW), gOut)]
#Save mesh function to pvd 1 = yes other = no
mesh_func_flag = 0
if mesh_func_flag == 1:
boundary_file = File("boundaries/boundaries.pvd")
boundary_file << bound_mshfunc
#############################################################
# Generate the DGFEM formulation for the operator
class Operator(HyperbolicOperator):
def __init__(self, mesh, V, bcs):
dim = mesh.geometric_dimension()
def F_c(U):
n_i = U[0]
n_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_u/n_i
p = n_i
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_u = as_vector([U[j] for j in range(1, dim+1)])
u = n_u/n_i
c = sqrt(1)
lambdas = [dot(u, n) - c, dot(u, n), dot(u, n) + c]
return lambdas
HyperbolicOperator.__init__(self, mesh, V, bcs, F_c, LocalLaxFriedrichs(alpha))
###################################################
po = Operator(mesh, V, bcs= bcs)
residual = dt*po.generate_fem_formulation(u_vec, v)
residual += dot(u_vec - un, v)*dx
#Write RHS as forcement term
f = as_vector((0,0,0))
residual += - dot(f,v)*dx(mesh)
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()
plot_flag = 1
#Plot the solution for every 0.1*N time-step 1= yes, else = no
for j in range(N):
un.assign(u_vec)
solver.solve(Problem(), u_vec.vector())
if j % int(0.1*N) == 0 and plot_flag == 1:
ni, n_u, n_v = u_vec
plt.colorbar(plot(ni))
plt.show()
Thanks.