Hi, I have problem wherein in know the final amount of substance (mass) in an object. I trying to set this final amount as a initial condition (concentration) on a boundary and observe how it diffuses through the model. Therefore, ideally when I integrate the concentration distributed across the mesh at the end of the simulation the amount should match the initial amount set at the boundary.
With a simple box mesh and a simple diffusion equation, the solution has an error of <0.5% with a relatively coarse mesh. I’ve implemented the initial condition as:
This works well. However, when I try to the implement the same idea for a more complex cylindrical boundary with a coupled PDE solving for mixed element I get an error of almost 9% for a very refine mesh. I’m not sure if I am overlooking something. I have implemented it like so:
As you have not provided a minimal code example (should be possible to run without adding additional code), I am not sure I follow what you want to obtain.
This is my interpretation of what you would like:
Given a mass M, you would like to distribute the mass along a boundary (as density).
Here is my minimal code example of what I think you would like to achieve:
from dolfin import (Function, FunctionSpace, Constant, assemble,
ds, dx, Point, DirichletBC, UnitSquareMesh)
for N in [10,100,500]:
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, "CG", 1)
u_init = Function(V)
mass = 3
boundary_area = assemble(Constant(1)*ds(domain=mesh))
density = mass/boundary_area
bc = DirichletBC(V, density, "on_boundary")
bc.apply(u_init.vector())
print("Mass:", assemble(u_init*ds), assemble(u_init*dx))
As you can observe, the mass distributed along the surface matches what I gave as an initial mass,
while the mass in the domain, as the width of the volume integral is going to zero, is tending to zero.
Thank you for getting back to me. Sorry for not providing MWE but my code has no errors and the only issue I have is of how to determine the value of density. Your interpretation is spot on but currently I have the total mass within the system so I need assemble(u*dx)) to equal to mass (concentration * volume). Mass is the experimental value that I have. I am not sure how to then determine the value of density? Such that setting the density as the initial condition (u_init) on the boundary would yeild a solution (u) which would result in assemble(u*dx)) = mass. Its a diffusion study.
What I have done is to divide the Area by the number of nodes and multiply the ‘Area per node’ by the number of nodes on the boundary. Then I’ve defined the initial condition to be the mass / area at the boundary. We are reducing a 3D problem to a 2D study.
Area per node would only be correct on meshes with uniform cell and facet size (structured meshes).
How about this approach?
from dolfin import (Function, FunctionSpace, Constant, assemble,
ds, dx, Point, DirichletBC, UnitSquareMesh)
for N in [10,100,500]:
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, "CG", 1)
mass = 3
boundary_area = assemble(Constant(1)*ds(domain=mesh))
# density = mass/boundary_area
u_n = Function(V)
bc = DirichletBC(V, 1, "on_boundary")
bc.apply(u_n.vector())
density = mass/assemble(u_n*dx)
bc_density = DirichletBC(V, density, "on_boundary")
u_init = Function(V)
bc_density.apply(u_init.vector())
print("Mass:", assemble(u_init*dx))
As you can observe here, the integral of u_init is constant equal to the mass, independent of the mesh size, due to the scaling of \frac{1}{\int u_n \mathrm{d} x}.
