CODEWITHSUNDEEP
numerical-methodspdefipy
alxg
alxg

Reputation: 368

How can i offset the Fipy mesh?

i'm solving some PDE's using Fipy and it would be really handy to translate the mesh so it doesn't start at 0 instead of modifying all my equations.

This is important since my convection term is dependent on x.

Example of my focker-planvk implementation:

from fipy import Viewer,PowerLawConvectionTerm,ConvectionTerm , Variable, FaceVariable, CellVariable, Grid1D, ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer
from fipy.tools import numerix
# from fipy.viewers import viewer
import numpy as np

# plt.close('all')
nx = 150
Lx=15

dx = Lx/nx
mesh = Grid1D(nx=nx,dx=dx)

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

s=0.5
D = s**2/2


timeStepDuration = 0.9 * dx**2 / (2 * D)
steps = 300
x = mesh.cellCenters[0]
t = timeStepDuration * steps

xc = mesh.cellCenters
s0=s/3

Fc=-1*(xc-5)

cl=0.2
tv = Variable(value=0)


# Initial Gaussian distribution
def Radb(x,t, a=1,l=0,s=0.1):
    std=np.sqrt(s**2/(2*np.abs(a+l*t)))
    r1=(x-8)**2
    r2=2*std**2
    return np.exp(-r1/r2) / np.sqrt(std**2*2*np.pi)  


phi.setValue(Radb(x.value,0,a=1,l=0,s=s0))
print('check initial integral:', np.trapz(phi.value,x.value))

Pt=np.zeros((nx,steps))

for step in range(steps - 1):
    # convCoeff=(0.1*x.value[step],)
    tv.setValue(step*timeStepDuration)
    "space/time dependent convection"
    Fc=-1*(xc-(la+cl*tv))

    eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D)-ConvectionTerm(CellVariable(mesh=mesh, value=[Fc]))
    eqI = TransientTerm() == DiffusionTerm(coeff=D)-ConvectionTerm(CellVariable(mesh=mesh, value=[Fc]))
    eqCN = eqX + eqI

    eqCN.solve(var=phi, dt=timeStepDuration)
    Pt[:,step]=  phi.value 
eqI.solve(var=phi, dt=timeStepDuration)
Pt[:,-1]=  phi.value 


if __name__ == '__main__':
     viewer = Viewer(vars=(phi),
                     datamin=0., datamax=1.)
     viewer.plot()

"veryfing conservation"
import matplotlib.pyplot as plt
tp=np.linspace(0,t,steps)
mass=[np.trapz(Pt[:,i],x.value) for i in range(steps)]
fig, ax = plt.subplots()
ax.plot(tp,mass,label=('Total probability'))
ax.axhline(y=1,ls='dashed',color='grey',label=('1'))
ax.axhline(y=0.9999999,ls='dotted',color='red',label=('0.9999999'))
plt.legend()
#%%

"Evolution animation"
from matplotlib.animation import FuncAnimation
fig, ax = plt.subplots()
line, = ax.plot(x, Pt[:,0],color='C0')
def update(i):
    line.set_ydata(Pt[:,i])
    return (line)
anim = FuncAnimation(fig, update, frames=range(steps-1), interval=40)

I would like to be able to define x such that it doesn't start at 0, since the next step is to use some nonlinear convection terms.

Upvotes: 0

Views: 49

Answers (1)

jeguyer
jeguyer

Reputation: 2484

You can translate a mesh by adding a vector to the mesh, e.g.,

mesh = Grid1D(nx=nx,dx=dx) + [[-5.]]
  :
xc = mesh.cellCenters
  :
Fc=-xc

Upvotes: 1

Related Questions