Hi there, I want to conduct a simple calculation for a circular structure at equilibrium by using level set method with lsmlib. The circular structure is shrinking at a normal velocity of Constant1*(1+kappa), where kappa is the curvature (1/radius), against a resistance, Constant2
v_n = C3*(Constant2 - Constant1*(1+kappa)) Constant1 and Constant2 are selected to ensure the whole structure is under equilibrium analytically. Since it's a circular structure with a uniform radius, the curvature should be isotropic and the structure should always keep the circular shape. Please run the test scripts anisotropy_test1.py. From the results, it can be seen some parts of the structure is shrinking and other are growing as shown in the plot attached. This indicate an anisotropy behavior in curvature calculation, which is done by phi.faceGrad.divergence. I re-ran the same simulation with a higher mesh resolution as shown in anisotropy_test2.py. However, similar results were observed. My question is how this anisotropy effect can be eliminated? As I am studying the effect of anisotropy at this moment, my results suffer from this error anisotropy. Thank you. Best, Yunbo
from fipy import DistanceVariable, CellVariable, Viewer, Grid2D
from pvd.pylsmlib import computeExtensionFields2d, computeDistanceFunction2d
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import figure, show
from matplotlib.pylab import rcParams
from numpy import reshape
L = 100.0 #
N = 100
dL = L / N
dt = 0.5
radius = 20.0
Decrease_0 = 0.544
Growth_factor = 0.548
R = 8.314
numberofStep = 20000
m = Grid2D(nx=N, ny=N, Lx=L, Ly=L)
phi = DistanceVariable(mesh=m, value=1.0, name='distance function')
vel = CellVariable(mesh=m, rank=1, elementshape=(1,), name='velocity')
extVel = CellVariable(mesh=m, name=r'extension velocity')
phi[:] = (m.x - L/2.)**2 + (m.y - L/2.)**2 - (radius)**2
phi[:] = computeDistanceFunction2d(phi.value, nx=m.nx, ny=m.ny, dx=m.dx, dy=m.dy)
for step in range(numberofStep):
Decrease = Decrease_0 * (1 + (phi.faceGrad.divergence)*0.1)
Net_growth = (Growth_factor-Decrease)
vel[:] = Net_growth
extVel[:] = vel[:]
tmpphi, extVel[:] = computeExtensionFields2d(phi.value, vel.value, nx=m.nx, ny=m.ny, dx=m.dx, dy=m.dy)
if step%100==0:
#plot the zero levelset
X = m.x[0:N]
Y = X
X1, Y1 = np.meshgrid(X, Y)
values = phi.value
values = values.reshape((N, N))
fig = plt.figure(1, figsize=(8,8))
ax = fig.add_subplot(111)
ax.contour(X1, Y1, values, 1,
linewidths=3, color='blue')
plt.savefig('../../plot_%.3d.pdf' % step)
plt.clf()
print step
phi[:] = phi[:] - dt * extVel[:]
from fipy import DistanceVariable, CellVariable, Viewer, Grid2D
from pvd.pylsmlib import computeExtensionFields2d, computeDistanceFunction2d
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import figure, show
from matplotlib.pylab import rcParams
from numpy import reshape
L = 100.0 #
N = 400
dL = L / N
dt = 0.5/16.
radius = 20.0
Decrease_0 = 0.544
Growth_factor = 0.5465
R = 8.314
numberofStep = 20000
m = Grid2D(nx=N, ny=N, Lx=L, Ly=L)
phi = DistanceVariable(mesh=m, value=1.0, name='distance function')
vel = CellVariable(mesh=m, rank=1, elementshape=(1,), name='velocity')
extVel = CellVariable(mesh=m, name=r'extension velocity')
phi[:] = (m.x - L/2.)**2 + (m.y - L/2.)**2 - (radius)**2
phi[:] = computeDistanceFunction2d(phi.value, nx=m.nx, ny=m.ny, dx=m.dx, dy=m.dy)
for step in range(numberofStep):
Decrease = Decrease_0 * (1 + (phi.faceGrad.divergence)*0.1)
Net_growth = (Growth_factor-Decrease)
vel[:] = Net_growth
extVel[:] = vel[:]
tmpphi, extVel[:] = computeExtensionFields2d(phi.value, vel.value, nx=m.nx, ny=m.ny, dx=m.dx, dy=m.dy)
if step%100==0:
#plot the zero levelset
X = m.x[0:N]
Y = X
X1, Y1 = np.meshgrid(X, Y)
values = phi.value
values = values.reshape((N, N))
fig = plt.figure(1, figsize=(8,8))
ax = fig.add_subplot(111)
ax.contour(X1, Y1, values, 1,
linewidths=3, color='blue')
plt.savefig('../../plot_%.3d.pdf' % step)
plt.clf()
print step
phi[:] = phi[:] - dt * extVel[:]
plot_900.pdf
Description: Adobe PDF document
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
