> Could you possible post your code so I can play around with it? I > can't think offhand what could be the issue or exactly what the error > is that you're calculating, it seems like the normalized residual or > something.
hi daniel, thanks for your reply. the way i was calculating the residual
is to export the solution on a grid, import that into mathematica,
interpolate in space and time and have mathematica calculate
||lhs-rhs|| / ||lhs|| where ||...|| is the L2 norm. what i saw was that
the relative error was low in the beginning but increased over time,
especially in places where the solution starts to grow exponentially.
i am sure this relative error can be gotten from fipy using
residualVector or some other methods of the terms but i have not been
able to understand the docs well enough to replicate that.
so i am posting pretty much the (streamlined but otherwise unchanged)
example script i am having trouble with below. h5py is needed only for
the export, you can comment that out if you dont export. Lambda3 is the
source term coefficient -- by increasing it's values the solution
becomes more rapidly growing and the errors also increase.
i hope the indentation survives emailing....
cheers, n.
from __future__ import division
import fipy as fp
from fipy import numerix as nx
# the module for the bare numpy functions
# this is essentially numpy!
NX = nx.NUMERIX
from time import sleep
# this is where we are going to write results
import os
datadir = os.path.join(
os.path.dirname(os.path.realpath(__file__)),
'pde_test_data/'
)
datadir = os.path.realpath(datadir)
try: os.mkdir(datadir)
except OSError: pass
# utility
def nsq(v):
return nx.dot(v, v)
class Domain(object):
upper = 300.
lower = 0
size = upper - lower
class Mesh(Domain):
def __init__(self):
delta = 1.5
Lx, Ly = (self.size,) * 2
dx, dy = (delta,) * 2
nx, ny = (self.size // delta,) * 2
self.mesh = fp.Grid2D(dx=dx, dy=dx, nx=nx, ny=ny, Lx=Lx, Ly=Ly)
# from Mma, the last iteration of the potential was this:
class Phi(object):
scale = 100.
x1, y1 = [100., 100]
x2, y2 = [180., 180]
s1 = 4.
s2 = 12.
dp = 3.
strength = 30.
sharpness_a = 120.
# the following should act on the position field x
# this should be a cell varable in fipy.
@classmethod
def itp(cls, z):
t = nx.tanh(z/cls.sharpness_a**2) + 1
return 0.5 * t
@classmethod
def scalephi(cls, x, y):
[d1, d2] = [(x-cls.x1)**2 + (y-cls.y1)**2,
(x-cls.x2)**2 + (y-cls.y2)**2]
[f2, f1] = [cls.itp(d1 - d2), cls.itp(d2 - d1)]
p1 = 0.5 * cls.s1 * d1 / cls.scale**2
p2 = 0.5 * cls.s2 * d2 / cls.scale**2 - cls.dp
return f1 * p1 + f2 * p2
@classmethod
def __call__(cls, x, y):
return cls.strength * cls.scalephi(x, y)
class Lambda3(object):
def __init__(self, curvature=None, direction=None, start=None):
if curvature is None: curvature = 0.2 / Domain.size**2
if direction is None: direction = nx.array([1, 1.])
direction /= nx.sqrt(nsq(nx.array(direction)))
if start is None: start = nx.array([0.,0])
max_pos = 0.5 * nx.sqrt(2) * Domain.size
max_ht = 1./2 * max_pos**2 * curvature # for 0 at the origin
max_ht = 0.1 * max_ht # tame the beast
(self.curvature, self.direction, self.max_pos, self.max_ht,
self.start) = (
curvature, direction, max_pos, max_ht, start)
def __call__(self, x, y):
distance2 = (NX.dot(self.direction,
(nx.array([x, y]).T - self.start).T) - self.max_pos)**2
lam = self.max_ht - 1./2 * self.curvature * distance2
return lam
class InitialRho(object):
'''Gaussian-like initial condition'''
c = [[100], [100]]
covm = NX.eye(2) * 30**2 # bare numpy
stm = NX.linalg.inv(covm)
Z = (2. * NX.linalg.det(covm))**(1./2)
@classmethod
def __call__(cls, x, y):
p = NX.array([x,y])
return 1. / cls.Z * NX.exp(-1./2 * NX.diag(
cls.stm.dot(p - cls.c).T.dot(p - cls.c)))
mesh = Mesh().mesh
x, y = mesh.cellCenters
x1d, y1d = x[:mesh.nx], y[::mesh.ny] # yes!
# diffusion
Dtensor = nx.eye(2) * 10.
# constant diff coeff fields
Dtensor = fp.CellVariable(mesh, 'dtensor', Dtensor)
# make the input fields
phi = fp.CellVariable(name=r'$\phi$', mesh=mesh, value=Phi()(x, y))
v = - phi.grad # on the cells, not the faces.
lambda_ = fp.CellVariable(name=r'$\lambda$', mesh=mesh,
value=Lambda3(curvature=1./Domain.size**2, direction=[1.,-1],
start=[0.,Domain.size])(x, y))
# terms
def diff(DD):
# Correction actually needed? in any case it seems to work.
# list of diffusion coefficients of different orders
return fp.DiffusionTermCorrection(coeff=[DD,])
def drift(v):
return - fp.PowerLawConvectionTerm(coeff=v)
def prolif(lambda_):
# implicit = multiplies the var
return fp.ImplicitSourceTerm(coeff=lambda_)
def lhs():
return fp.TransientTerm()
# construct the density
rho = fp.CellVariable(name=r'$\rho$', mesh=mesh, hasOld=True)
# boundary conditions
border = mesh.exteriorFaces
bc = 'dirichlet'
bc = 'no-flux'
# Dirichlet bc
if bc == 'dirichlet':
rho.constrain(0., where=border)
# no-flux: set convection and diffusion 0
elif bc == 'no-flux':
pass # no flux is standard b.c. in fipy.
else: raise ValueError("invalid bc type")
eq_diffprolif = lhs() == diff(Dtensor) + drift(v) + prolif(lambda_)
# initial condition
ic = fp.CellVariable(name=r'$\rho_0$', mesh=mesh, value=InitialRho()(x,y))
variants = { 'diffprolif': (eq_diffprolif, Dtensor, phi, lambda_, ic) }
# set up visualization for it
viewer = fp.Viewer(rho, datamax = 3. * ic.max(), datamin = -0.1 * ic.max())
# solver parameters
class Solve(object):
sampling = 5
step = sampling / 5
# T = 10
T = 80
sample_t = NX.arange(0, T + sampling / 2., sampling)
sweeps = 0.01
# customSolver = fp.LinearCGSSolver()
customSolver = None
assert NX.isclose(0, NX.modf(sampling / step)[0])
# do it for all the equations of interest:
snapshots, residuals = {}, {}
eq, _, _, _, ic_ = variants['diffprolif']
# initialize
rho.setValue(ic_.value)
snapshots = []
residuals = []
# first
snapshots.append(rho.value.copy())
viewer.plot()
for (st, stnext) in zip(Solve.sample_t, Solve.sample_t[1:]):
for tt in nx.arange(st, stnext, Solve.step):
rho.updateOld()
res = 1.
while res > Solve.sweeps:
res = eq.sweep(var=rho, dt=Solve.step,
solver=Solve.customSolver)
residuals.append(res)
snapshots.append(rho.value.copy())
sleep(0.1)
viewer.plot()
# Now finally we have to export this. Use hdf5 as from Mma to reproduce
the output.
# sampling for export
sampling_delta = 1.5
ex_x1, ex_y1 = NX.arange(1., 280, sampling_delta), NX.arange(1., 280,
sampling_delta)
ex_x, ex_y = NX.meshgrid(ex_x1, ex_y1)
# interpolate the intermediates at these values!
def sample_field(cellvar, value=None, xs=None, ys=None, normalize=False):
from scipy.interpolate import RectBivariateSpline
val = cellvar.value if value is None else value
f = val.reshape(cellvar.mesh.nx, -1)
x1 = cellvar.mesh.x[:cellvar.mesh.nx]
y1 = cellvar.mesh.y[::cellvar.mesh.ny] # sic.
field = RectBivariateSpline(x1, y1, f,
bbox=[Domain.lower, Domain.upper] * 2, s=0)
xs = ex_x1 if xs is None else xs
ys = ex_y1 if ys is None else ys
sf = field(xs, ys)
area_el = (xs[1]-xs[0]) * (ys[1]-ys[0])
int_ = sf.sum() * area_el
return (sf / int_, int_) if normalize else sf
# write to HDF5 file
import h5py as h5
outfile = os.path.join(datadir, 'twohump.h5')
eq, Dtens, ph, lam, ic_ = variants['diffprolif']
Dt = Dtens[:,:,0] # get back constant array from cell variable
with h5.File(name=outfile, mode='w') as f:
f['phi'] = sample_field(ph)
f['lambda'] = sample_field(lam)
dNs = [sample_field(rho, value=snap, normalize=True)
for snap in snapshots]
f['density'] = NX.array([dn[0] for dn in dNs])
f['Nt'] = NX.array([dn[1] for dn in dNs])
f['diffcoeff'] = NX.diag(Dt)
f['timeline'] = Solve.sample_t
assert f['density'].shape[0] == f['timeline'].shape[0]
f['xs'], f['ys'] = ex_x1, ex_y1
0xF0A5C638.asc
Description: application/pgp-keys
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
