CODEWITHSUNDEEP
pythongeometryfipygmsh
Tiffany Overdick
Tiffany Overdick

Reputation: 11

How to make a diffusion from the center of the circle to the edge on gmsh? [python]

I am trying to simulate the diffusion of a nutrient in a tumor as a function of space and time on python using the gmsh tool. So I need the initial concentration to be in the center, then diffusivity to spread it outward. Here is the complete code, the parameters are not right yet but I'm already trying to get the diffusivity from the inside to the outside and not the other way around.

Here is my code:

from fipy import CellVariable, Gmsh2D, ExponentialConvectionTerm, TransientTerm, DiffusionTerm, ImplicitSourceTerm, Viewer
from fipy.tools import numerix
import numpy as np
import matplotlib.pyplot as plt


# Define some parameters for the creation of the mesh

cellSize = 0.05
radius = 1.


# Define the grid/mesh

mesh = Gmsh2D('''
              cellSize = %(cellSize)g;
              radius = %(radius)g;
              Point(1) = {0, 0, 0, cellSize};
              Point(2) = {-radius, 0, 0, cellSize};
              Point(3) = {0, radius, 0, cellSize};
              Point(4) = {radius, 0, 0, cellSize};
              Point(5) = {0, -radius, 0, cellSize};
              Circle(6) = {2, 1, 3};
              Circle(7) = {3, 1, 4};
              Circle(8) = {4, 1, 5};
              Circle(9) = {5, 1, 2};
              Line Loop(10) = {6, 7, 8, 9};
              Plane Surface(11) = {10};
              ''' % locals())


# Define the model variable and set the boundary conditions

phi = CellVariable(name = "solution variable",
                   mesh = mesh,
                   value = 0.) 
#phi = CellVariable(mesh)

X, Y =  mesh.faceCenters 

dr = np.linalg.norm(mesh.faceCenters, axis=0)
mask = (dr<50) * mesh.exteriorFaces
phi.constrain(1, mask)

mask = (dr>50) * mesh.exteriorFaces
phi.constrain(0, mask)

viewer = None
from fipy import input
if __name__ == '__main__':
  viewer = Viewer(vars=phi, datamin=-1, datamax=1.)
  viewer.plotMesh()

# Define and then solve the equation
D = 1.
eq = TransientTerm() == DiffusionTerm(coeff=D) 

timeStepDuration = 10 * 0.9 * cellSize**2 / (2 * D)
steps = 10
from builtins import range
for step in range(steps):
  eq.solve(var=phi,
           dt=timeStepDuration)
  if viewer is not None:
    viewer.plot()

Thanks a lot !

Upvotes: 1

Views: 185

Answers (1)

jeguyer
jeguyer

Reputation: 2484

You describe an initial condition, but then your code seems to be trying to set a boundary condition.

I'm also confused by your geometry. Your mesh is defined to have a radius of 1, but then you attempt to constrain the value based on whether the radius is less than or greater to 50. Note that all exteriorFaces are (within the resolution of the mesh discretization) at r=1.

If you want an initial nutrient concentration to be at some radius less than the exterior boundary, then, e.g, you'll want to assign cells, not faces:

r = np.linalg.norm(mesh.cellCenters, axis=0)
phi.setValue(1., where=r<0.5)

Whether you constrain the exterior faces is up to you. With no constraint, the natural boundary condition is no-flux.

Also see https://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.anisotropy.html.

Upvotes: 1

Related Questions