Hello everyone, I am new to FiPy and interested in solving a certain curvature model. I thought using FiPy will make my life easier, but unfortunately I have been struggling with this for quite some time trying different versions of introducing the equations into FiPy.
I decided to try this community, hoping to get some insight on what I am doing wrong. The model includes three variables (two dynamic, one auxiliary) and three coupled equations running on a 2D grid. Please see a summary of the model in the attached PDF file. I am also attaching a minimalistic version of the code. For this version the error message is: "[0] Matrix is missing diagonal entry 2500" (with 2500 being nx*ny). Many thanks in advance, Matan -- To unsubscribe from this group, send email to [email protected] View this message at https://list.nist.gov/fipy --- To unsubscribe from this group and stop receiving emails from it, send an email to [email protected].
CurvatureModelEquations.pdf
Description: Adobe PDF document
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from fipy import CellVariable, Grid2D, GaussianNoiseVariable, DiffusionTerm, TransientTerm, ImplicitSourceTerm, Viewer
from fipy.tools import numerix
# create a 2D solution mesh
nx = ny = 50
dx = dy = 0.1
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)
phi0 = -1.
psi0 = phi0**3 - phi0
sigma0 = 0.001
# create variables
phi = CellVariable(name=r"$\phi$", mesh=mesh, value=phi0, hasOld=True)
psi = CellVariable(name=r"$\psi", mesh=mesh, value=psi0, hasOld=True)
sigma = CellVariable(name=r"$\sigma", mesh=mesh, value=sigma0) #hasOld=True
# Initial condition: initiate phi as a circle of 1. in a field of -1.
r0 = np.sqrt(nx**2 + ny**2)/4 #radius
x0 = nx/2 #position of center of circle
y0 = ny/2
x = numerix.arange(0, nx)
y = numerix.arange(0, ny)
initCond = -1*numerix.ones((y.size, x.size))
mask = (x[numerix.newaxis,:]-x0)**2 + (y[:,numerix.newaxis]-y0)**2 < r0**2
initCond[mask] = 1.
phi[:] = initCond.reshape((int(nx*ny)))
# Parameters
kappa = 0.01
epsilon = 0.2
a0 = 10
# Terms of Equations
Diff1 = DiffusionTerm(coeff=12*kappa*phi*psi, var=phi)
Diff2 = DiffusionTerm(coeff=2*kappa*(3*phi**2 - 1), var=psi)
Diff3 = DiffusionTerm(coeff=(-2*kappa*epsilon**2, 1.), var=psi)
Source1 = ImplicitSourceTerm(coeff=phi**2 - 1., var=phi)
Diff4 = DiffusionTerm(coeff=-epsilon**2, var=phi)
# Equations
eq1 = (TransientTerm(var=phi) == Diff1 + Diff2 + Diff3)
eq2 = (ImplicitSourceTerm(coeff=1., var=psi) == Source1 + Diff4)
eq3 = (TransientTerm(var=sigma) == phi.grad.dot(phi.grad)/2 - a0)
eq = eq1 & eq2 & eq3
# Use exponentially increasing time to accelarate evolution as system relaxes towards the equilibrium state
dexp = -5
elapsed = 0.
duration = 0.5
# solving the equation
while elapsed < duration:
dt = min(100, numerix.exp(dexp))
elapsed += dt
dexp += 0.01
eq.solve(dt=dt)
