Python Fipy explicit terms with derivatives and variable diffusion constant
I am trying to solve four coupled PDEs (two scalar density fields and one vector field- but I solve it component wise).
where the second term in first equation is given by
Now, in the Fipy code (the \rho and the vector components are denoted as by m, px and py respectively) I am writing them down as:
mesh = PeriodicGrid2D(dx, dy, nx, ny)
# Variables to use
c = CellVariable(name='c', mesh=mesh, hasOld=1)
m = CellVariable(name='m', mesh=mesh, hasOld=1)
px = CellVariable(name='px', mesh=mesh, hasOld=1)
py = CellVariable(name='py', mesh=mesh, hasOld=1)
x_hat = [1.0, 0.0]
y_hat = [0.0, 1.0]
#-------------------------------------------------------------------------
eq_c = (TransientTerm(var=c) ==
W*(c.grad.dot(x_hat)*m.grad.dot(x_hat) + c.grad.dot(y_hat)*m.grad.dot(y_hat))
- ExponentialConvectionTerm(coeff=x_hat * px*v0 + y_hat * py*v0, var=c)
+ DiffusionTerm(coeff=Da, var=c) + DiffusionTerm(coeff=W*c, var=m) )
eq_m = (TransientTerm(var=m) ==
kb*(c/(c+ch)) - ImplicitSourceTerm(coeff=ku, var=m)
- ExponentialConvectionTerm(coeff=x_hat * px*v0 + y_hat * py*v0, var=m)
+ DiffusionTerm(coeff=Dm, var=m) )
eq_px = (TransientTerm(var=px) ==
DiffusionTerm(coeff=K, var=px) - zeta0*c.grad.dot(x_hat) + zeta*m.grad.dot(x_hat)
+ ImplicitSourceTerm(coeff=((c/cs)-1 -((c/cs)+1)*(px*px+py*py)), var=px) )
eq_py = (TransientTerm(var=py) ==
DiffusionTerm(coeff=K, var=py) - zeta0*c.grad.dot(y_hat) + zeta*m.grad.dot(y_hat)
+ ImplicitSourceTerm(coeff=((c/cs)-1 -((c/cs)+1)*(px*px+py*py)), var=py) )
# Couple them into one big equation
eq = eq_c & eq_m & eq_px & eq_py
Is the Fipy implementation of the second term in c-equation correct? The code does not thow up any error but I am not getting expected results should be an aster as described in this paper. But I get solutions that show aster initially but diffusion takes over and in long term I get homogeneous solutions. I have seen the previous related questions here and here, but not entirely sure how to make use of them. Thanks in advance and sorry for the poor figures!
Here is a minimal working code:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
First version with no polarity dynamics, no equation for n-field
"""
from fipy import (CellVariable, PeriodicGrid2D, TransientTerm, DiffusionTerm, ExplicitDiffusionTerm,
UniformNoiseVariable, Variable, FaceVariable, ImplicitSourceTerm, ExponentialConvectionTerm)
import matplotlib.pyplot as plt
import matplotlib
from matplotlib import cm
matplotlib.rcParams.update({'font.size': 16})
import numpy as np
def run_sim(savename,K,zeta):
#------------ parameters -------------
Da=1.0
Dm=1.0
W=-30.0
kb=1.0/4.0
ku=1.0/8.0
ch=1.0/4.0
cs=1.0/64.0
zeta0=1.0
# zeta=1.0
v0=0.1
#-------------- Mesh --------------------
duration=1
dt=0.5
# Define mesh size and number of points
nx = 60
ny = nx
L = 20
dx = L / nx
dy = dx
mesh = PeriodicGrid2D(dx, dy, nx, ny)
xval = mesh.x
yval = mesh.y
# Variables to use
c = CellVariable(name='c', mesh=mesh, hasOld=1)
m = CellVariable(name='m', mesh=mesh, hasOld=1)
px = CellVariable(name='px', mesh=mesh, hasOld=1)
py = CellVariable(name='py', mesh=mesh, hasOld=1)
# Initial condition + some random noise
c.setValue( 0.1*np.exp(-0.25*(mesh.x-L/2)*(mesh.x-L/2)-0.25*(mesh.y-L/2)*(mesh.y-L/2)) )
m.setValue( UniformNoiseVariable(mesh=mesh, minimum=0.4*(kb/ku), maximum=0.5*(kb/ku)) )
px.setValue( UniformNoiseVariable(mesh=mesh, minimum=0.0, maximum=0.1) )
py.setValue( UniformNoiseVariable(mesh=mesh, minimum=0.0, maximum=0.1) )
x_hat = [1.0, 0.0]
y_hat = [0.0, 1.0]
#-------------------------------------------------------------------------
eq_c = (TransientTerm(var=c) ==
W*(c.grad.dot(x_hat)*m.grad.dot(x_hat) + c.grad.dot(y_hat)*m.grad.dot(y_hat))
- ExponentialConvectionTerm(coeff=x_hat * px*v0 + y_hat * py*v0, var=c)
+ DiffusionTerm(coeff=Da, var=c) + DiffusionTerm(coeff=W*c, var=m) )
eq_m = (TransientTerm(var=m) ==
kb*(c/(c+ch)) - ImplicitSourceTerm(coeff=ku, var=m)
- ExponentialConvectionTerm(coeff=x_hat * px*v0 + y_hat * py*v0, var=m)
+ DiffusionTerm(coeff=Dm, var=m) )
eq_px = (TransientTerm(var=px) ==
DiffusionTerm(coeff=K, var=px) - zeta0*c.grad.dot(x_hat) + zeta*m.grad.dot(x_hat)
+ ImplicitSourceTerm(coeff=((c/cs)-1 -((c/cs)+1)*(px*px+py*py)), var=px) )
eq_py = (TransientTerm(var=py) ==
DiffusionTerm(coeff=K, var=py) - zeta0*c.grad.dot(y_hat) + zeta*m.grad.dot(y_hat)
+ ImplicitSourceTerm(coeff=((c/cs)-1 -((c/cs)+1)*(px*px+py*py)), var=py) )
# Couple them into one big equation
eq = eq_c & eq_m & eq_px & eq_py
elapsed = 0.0
while elapsed < duration:
c.updateOld()
m.updateOld()
px.updateOld()
py.updateOld()
elapsed += dt
res = 1e5
old_res = res * 2
step = 0
while res > 1e-5 and step < 15 and old_res / res > 1.01:
old_res = res
res = eq.sweep(dt=dt)
step += 1
if __name__ == '__main__':
"""
comments
"""
path = 'result/'
def name(K, zeta):
return 'K_{:.2f}_zeta_{:.2f}'.format(K, zeta)
Ks = [1]
zetas = [30.0]
for K in Ks:
for zeta in zetas:
run_sim(name(K, zeta), K=K, zeta=zeta)
Solution
Is the Fipy implementation of the second term in c-equation correct?
No. The term DiffusionTerm(coeff=W*c, var=m) represents 
. There is no reason to run the chain rule to obtain the explicit parts; it is, in fact, undesirable to do so.
Further converting the orientation to a rank-1 CellVariable results in somewhat cleaner equations. FiPy cannot couple scalar c and m with vector p, but the px and py equations weren't coupled, anyway.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
First version with no polarity dynamics, no equation for n-field
"""
from fipy import (CellVariable, Grid2D, PeriodicGrid2D, TransientTerm, DiffusionTerm, ExplicitDiffusionTerm,
UniformNoiseVariable, Variable, FaceVariable, ImplicitSourceTerm, ExponentialConvectionTerm, Viewer)
import matplotlib.pyplot as plt
import matplotlib
from matplotlib import cm
matplotlib.rcParams.update({'font.size': 16})
import numpy as np
def run_sim(savename,K,zeta):
#------------ parameters -------------
Da=1.0
Dm=1.0
W=-30.0
kb=1.0/4.0
ku=1.0/8.0
ch=1.0/4.0
cs=1.0/64.0
zeta0=1.0
# zeta=1.0
v0=0.1
#-------------- Mesh --------------------
duration=1
dt=0.05
# Define mesh size and number of points
nx = 60
ny = nx
L = 20
dx = L / nx
dy = dx
mesh = Grid2D(dx, dy, nx, ny)
xval = mesh.x
yval = mesh.y
# Variables to use
c = CellVariable(name='c', mesh=mesh, hasOld=1)
m = CellVariable(name='m', mesh=mesh, hasOld=1)
p = CellVariable(name='p', mesh=mesh, rank=1, hasOld=True)
# Initial condition + some random noise
c.setValue( 0.1*np.exp(-0.25*(mesh.x-L/2)*(mesh.x-L/2)-0.25*(mesh.y-L/2)*(mesh.y-L/2)) )
m.setValue( UniformNoiseVariable(mesh=mesh, minimum=0.4*(kb/ku), maximum=0.5*(kb/ku)) )
p[0] = UniformNoiseVariable(mesh=mesh, minimum=0.0, maximum=0.1)
p[1] = UniformNoiseVariable(mesh=mesh, minimum=0.0, maximum=0.1)
x_hat = [1.0, 0.0]
y_hat = [0.0, 1.0]
#-------------------------------------------------------------------------
eq_c = (TransientTerm(var=c) ==
- ExponentialConvectionTerm(coeff=p*v0, var=c)
+ DiffusionTerm(coeff=Da, var=c) + DiffusionTerm(coeff=W*c, var=m) )
eq_m = (TransientTerm(var=m) ==
ImplicitSourceTerm(coeff=kb/(c+ch), var=c) - ImplicitSourceTerm(coeff=ku, var=m)
- ExponentialConvectionTerm(coeff=p*v0, var=m)
+ DiffusionTerm(coeff=Dm, var=m) )
eq_p = (TransientTerm(var=p) ==
DiffusionTerm(coeff=[[[K, 0],
[0, K]]], var=p) - zeta0*c.grad + zeta*m.grad
+ ImplicitSourceTerm(coeff=((c/cs)-1 -((c/cs)+1)*p.mag**2), var=p) )
# Couple them into one big equation
eq = eq_c & eq_m
elapsed = 0.0
viewer = Viewer(vars=c)
viewer.plot()
while elapsed < duration:
c.updateOld()
m.updateOld()
p.updateOld()
elapsed += dt
res = 1e5
old_res = res * 2
step = 0
while res > 1e-5 and step < 15 and old_res / res > 1.01:
old_res = res
res = eq.sweep(dt=dt)
eq_p.sweep(dt=dt)
step += 1
viewer.plot()
if __name__ == '__main__':
"""
comments
"""
path = 'result/'
def name(K, zeta):
return 'K_{:.2f}_zeta_{:.2f}'.format(K, zeta)
Ks = [1]
zetas = [30.0]
for K in Ks:
for zeta in zetas:
run_sim(name(K, zeta), K=K, zeta=zeta)
- Activating Anaconda Environment in VsCode
- How to implement the Softmax function in Python?
- Dataframe - Select the minimum set of rows to cover all possible values of each columns
- Non-equi join in polars
- Tkinter and matplotlib - NavigationToolbar2Tk - cursor position precision
- Media Image is not loaded in the output of Django Project
- selenium upload part not working on windows. Path related?
- How to encode UTF8 filename for HTTP headers? (Python, Django)
- Split columns in Csv file
- In the Django admin site, how do I change the display format of time fields?
- Condition in pytest fixture depending on some test-suite results
- How to make a class JSON serializable
- Apache Airflow unit and integration test
- TypeError: 'MultiPolygon' object is not iterable with Plotly Choropleth
- How do I get the current time in milliseconds in Python?
- Python: Decoding base64 encoded strings
- How to fix the error when try to plot data with pandas?
- Python Image - Finding largest branch from image skeleton
- Enumerate is Returning A None Type
- Python finite difference functions?
- How to extract hex color codes from a colormap
- Regex: Get "London, UK" from "Evolution Recruitment (Agency) (London, UK)"
- Visual Studio Code - How to add multiple paths to python path?
- What are the Spark transformations that causes a Shuffle?
- How to find the second lowest lists into a nested list by their second value?
- python-docx: Parse a table to Panda Dataframe
- Deposit and Withdraw on Account
- Why does this code snippet give the error list object is not callable?
- How to find out if argparse argument has been actually specified on command line?
- Print current call stack from a method in code