I am solving a simple diffusion equation for a mesh with 3 different subdomains. I implemented an initial condition for the concentration for subdomain 1 c=1 and for subdomains 2 and 3 c=0.
The code runs without any errors. What I don’t understand tough is, that there is no flux between the regions with c=1 and c=0. Everything is static.
I would really appreciate any help, because I guess that I did something wrong in my code, because physically I would expect a mixing of the concentrations in the subdomains. In the picture you can see the subdomains. The circles are the subdomains 2 and 3 and everything else is the subdomain 1.
Thanks for your help in advance!
Here is a minimal example of my code:
from fenics import *
import matplotlib.pyplot as plt
D = 5*10**(-10) # Diffusioncoeff. m^2/s
D = Constant(D)
c_ini_f = 1 # initial concentration g/l
# Time variables
T = 2 # final time
num_steps = 20 # number of time steps
dt = T / num_steps # time step size
# Create mesh and define function space
mesh_file = "Netze/GewebeProbe_v1_in_m.h5"
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), mesh_file, "r")
hdf.read(mesh, "/mesh", False)
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
hdf.read(subdomains, "/markers")
hdf.close()
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, P1)
c_f = TrialFunction(V)
dc_f = TestFunction(V)
class InitialCondition(UserExpression): # UserExpression instead of Expression
def __init__(self, subdomains, **kwargs):
super().__init__(**kwargs)
self.subdomains = subdomains
def eval_cell(self, values, x, cell):
if self.subdomains[cell.index] == 1:
values[0] = c_ini_f
elif self.subdomains[cell.index] == 2:
values[0] = 0
elif self.subdomains[cell.index] == 3:
values[0] = 0
c0_f = Function(V)
c_init_sub=InitialCondition(subdomains, degree=0)
c0_f.interpolate(c_init_sub)
Diff_eqn = dc_f*(c_f-c0_f)/dt *dx + D*inner(grad(c_f),grad(dc_f))*dx
# Aufteilen der Gleichung in linke und rechte Seite
l = lhs(Diff_eqn)
r = rhs(Diff_eqn)
t = 0 # Zeit
c_f = Function(V)
##############################################################################
# Lösungsschleife #
##############################################################################
for n in range(num_steps):
# Update current time
t += dt
solve(l == r, c_f,
solver_parameters=dict(linear_solver='mumps'))
assign(c0_f, c_f)
plotKonz = plot(c_f,mode ='color')
plt.colorbar(plotKonz)
plt.show()
Dear Martin,
I ve got same problem with low constants like you in diffusion problem. Firstly, check your time step if you are using dt too large it goes to stable state. it is proximate ~ L^2/(2*D) where L= max dimension.
If not. It is probably because you are multiplying zero by zeros.
Check your dimensions, if it is too low, change units like mm, or micrometres or pm.
I recommend you admen your system.
If the physical constants is too low, I recommend you use “decimal”
from decimal import Decimal
D=float(Decimal(5.E-10))
D = Constant(D)
if nothing works, post one simply geometry similar to your mesh_file. I will check.
Thank you for your answer!
I changed my time step from dt=0.01 s to 0.001 s and 0.0005 s but sadly that didn’t help.
Speaking of the units. I have to admit, that I am not sure, which unit fenics is using for the mesh. When I check with Mesh_Koord = mesh.coordinates() and then look at the biggest entries (which would be my left corners): nx1=np.max(Mesh_Koord) = 10.
Is 10 in mm?
Because if it is, my D was way to small in the beginning.
So I also tried various Diffusioncoefficients, also with what you recommended
:
ranging from D = 1*e-04 to D = 1 but unfortunatly, the result was always the same, that I didn’t see any mixing.
It would be very kind of you, if you could check the problem with my mesh_file, is there a way for me to send my h5-file to you?
I recommend you to check if the circles are connected in domain.
I ve done one problem with same equation with low constants and it is working well.
you have to replot in each step loop
if your D is in m^2/s… your dimensions should be in meters
and your concentration should be in kg/m^3 or grams/litre
It looks like that your doman is 10 meters X 10 meters.
For stable state it will you be (10m)^2/(2*D)=2E+11 seconds or 6337.74777 years.
Maybe your dimensions are in mm and you have to multiply mesh by 1E-3!!
I think the problem might really be, that my circles are not connected to the domain. I will try to figure out how to fix this.
Could you tell me how to achieve this in fenics?
Because I created my mesh the following way:
I created the Geometry in SolidWorks, saved it as a Step-File, defined the physical surfaces in Gmsh and then meshed it and than converted it to a +h5 file for fenics. And the elements are physically connected to each other. I’m not sure what to change in Gmsh so fenics regards the subdomains as physically connected.
Hello,
I have a general question. I always see in examples that, for the definition of a constant, like here with the Diffusioncoefficient, first the value ist defined
D=float(Decimal(5.E-10))
and then it is defined as a fenics-Constant:
D = Constant(D)
I didn’t find an explanation, why you should do that. Is it just better/more efficient coding, or does it changes how the constant is treated while solving?
It changes how the compiled code, the integration kernel is generated. By encapsulating a constant with the Constant class, you can vary this parameter without dolfin having to do just in time compilation
