I solve diffusion equation, but it don’t work. My code:
from dolfin import *
mesh = UnitSquareMesh(40, 40)
F = FunctionSpace(mesh, 'DG', 1)
Rho = FunctionSpace(mesh, 'DG', 1)
V = VectorFunctionSpace(mesh, 'P', 1)
T = 10.0 # final time
num_steps = 10000 # number of time steps
dt = T / num_steps # time step size
k= Constant(dt)
u_n = Constant((0,0))
def boundary(x, on_boundary):
return on_boundary
f_D = Expression('0', degree = 1)
bc = DirichletBC(F, f_D, boundary)
bcf = [bc]
f_expr_init = Expression('(x[0] >= 0.4 && x[0] <= 0.6 && x[1] >= 0.4 && x[1] <= 0.6) ? 1.0 : 0.0', degree=1)
f_n = Function(F)
f_n = project(f_expr_init, F)
rho_1 = Constant(1021)
rho_2 = Constant(1117)
def rho_f(ff):
return project(rho_1*rho_2/(ff*rho_2 + (1 - ff)*rho_1), Rho)
rho = rho_f(f_n)
f = TrialFunction(F)
fi = TestFunction(F)
D = Constant(1e-10)
F4 = rho * dot((f - f_n) / k, fi)*dx \
+ div(rho * f * u_n)*fi*dx \
- D * inner(rho * grad(f), grad(fi))*dx
a4 = lhs(F4)
L4 = rhs(F4)
A4 = assemble(a4)
[bc.apply(A4) for bc in bcf]
t = 0
f_ = Function(F)
#num_steps
for n in range(5000):
t += dt
b4 = assemble(L4)
[bc.apply(b4) for bc in bcf]
solve(A4, f_.vector(), b4)
rho = rho_f(f_)
f_n.assign(f_)
max_f = np.abs(np.array(f_n.vector())).max()
print('t = %.3f: max_f = %.10f' % (t, max_f))
Max f (concentration) for each iteration > 1.
Please be more precise with your code. What do you imply by it not working?
As far as I can tell, the main issue is that you are trying to use a DG 0 function space for your state variable. This does not make sense, as the gradient of a DG 0 space is always 0.
See for instance: Bitbucket for an example of how to write a DG variational formulation
However, I see that there is no gradient due to the use of such elements. If i will not use DG, what elements should I use to set the square in the center with a value of one, and the rest of the field to zero? Sorry for my silly quetions.
If you use a DG-1 space, you need to modify your variational form such that the degrees of freedom on each cell is coupled, as shown in the demo.
If you use a CG-1 space, you can get an approximate initial condition in the CG 1 space by interpolating the expression, i.e.:
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(40, 40)
F = FunctionSpace(mesh, 'CG', 1)
Rho = FunctionSpace(mesh, 'DG', 1)
V = VectorFunctionSpace(mesh, 'P', 1)
T = 10.0 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
k= Constant(dt)
u_n = Constant((0,0))
def boundary(x, on_boundary):
return on_boundary
f_D = Expression('0', degree = 1)
bc = DirichletBC(F, f_D, boundary)
bcf = [bc]
f_expr_init = Expression('(x[0] >= 0.4 && x[0] <= 0.6 && x[1] >= 0.4 && x[1] <= 0.6) ? 1.0 : 0.0', degree=1)
f_n = Function(F)
f_n = interpolate(f_expr_init, F)
rho_1 = Constant(1021)
rho_2 = Constant(1117)
def rho_f(ff):
return project(rho_1*rho_2/(ff*rho_2 + (1 - ff)*rho_1), Rho)
rho = rho_f(f_n)
f = TrialFunction(F)
fi = TestFunction(F)
D = Constant(1e-2)
F4 = rho * dot((f - f_n) / k, fi)*dx \
+ div(rho * f * u_n)*fi*dx \
+ D * inner(rho * grad(f), grad(fi))*dx
a4 = lhs(F4)
L4 = rhs(F4)
A4 = assemble(a4)
[bc.apply(A4) for bc in bcf]
out = File("u.pvd")
t = 0
f_ = Function(F)
out << f_n, 0
#num_steps0
for n in range(5000):
t += dt
b4 = assemble(L4)
[bc.apply(b4) for bc in bcf]
solve(A4, f_.vector(), b4)
rho = rho_f(f_)
f_n.assign(f_)
out << f_n, t
max_f = np.abs(np.array(f_n.vector())).max()
print('t = %.3f: max_f = %.10f' % (t, max_f))
1 Like
I replaced boundary conditions, but i did’t see it:
def boundary(x, on_boundary):
if on_boundary:
if x[0] <= 0.2 and x[1] <= 0.2:
return True
else:
return False
else:
return False
f_D = Expression('1', degree = 1)
bc = DirichletBC(F, f_D, boundary)
bcf = [bc]
#Expression('(x[0] >= 0.4 && x[0] <= 0.6 && x[1] >= 0.4 && x[1] <= 0.6) ? 1.0 : 0.0', degree=1)
f_expr_init = Expression('0', degree = 1)
f_n = Function(F)
f_n = interpolate(f_expr_init, F)
For task:
My f always 0.0000
I cannot reproduce this with the following code. Please have a careful look at your implementation before reposting.
from dolfin import *
import numpy as np
mesh = UnitSquareMesh(40, 40)
F = FunctionSpace(mesh, 'CG', 1)
Rho = FunctionSpace(mesh, 'DG', 1)
V = VectorFunctionSpace(mesh, 'P', 1)
T = 10.0 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
k= Constant(dt)
u_n = Constant((0,0))
def boundary(x, on_boundary):
return on_boundary
# f_D = Expression('0', degree = 1)
# bc = DirichletBC(F, f_D, boundary)
# bcf = [bc]
# f_expr_init = Expression('(x[0] >= 0.4 && x[0] <= 0.6 && x[1] >= 0.4 && x[1] <= 0.6) ? 1.0 : 0.0', degree=1)
# f_n = Function(F)
# f_n = interpolate(f_expr_init, F)
def boundary(x, on_boundary):
if on_boundary:
if x[0] <= 0.2 and x[1] <= 0.2:
return True
else:
return False
else:
return False
f_D = Expression('1', degree = 1)
bc = DirichletBC(F, f_D, boundary)
bcf = [bc]
#Expression('(x[0] >= 0.4 && x[0] <= 0.6 && x[1] >= 0.4 && x[1] <= 0.6) ? 1.0 : 0.0', degree=1)
f_expr_init = Expression('0', degree = 1)
f_n = Function(F)
f_n = interpolate(f_expr_init, F)
rho_1 = Constant(1021)
rho_2 = Constant(1117)
def rho_f(ff):
return project(rho_1*rho_2/(ff*rho_2 + (1 - ff)*rho_1), Rho)
rho = rho_f(f_n)
f = TrialFunction(F)
fi = TestFunction(F)
D = Constant(1e-2)
F4 = rho * dot((f - f_n) / k, fi)*dx \
+ div(rho * f * u_n)*fi*dx \
+ D * inner(rho * grad(f), grad(fi))*dx
a4 = lhs(F4)
L4 = rhs(F4)
A4 = assemble(a4)
[bc.apply(A4) for bc in bcf]
out = File("u.pvd")
t = 0
f_ = Function(F)
out << f_n, 0
#num_steps0
for n in range(5000):
t += dt
b4 = assemble(L4)
[bc.apply(b4) for bc in bcf]
solve(A4, f_.vector(), b4)
rho = rho_f(f_)
f_n.assign(f_)
out << f_n, t
max_f = np.abs(np.array(f_n.vector())).max()
print('t = %.3f: max_f = %.10f' % (t, max_f))
