Re: reaction-advection-diffusion; problems making the solution converge

[Nils Becker]

> I'm not sure if using ||lhs|| is a good way to normalize the residual.
> Surely that could be very small. Is ||lhs - rhs|| getting smaller
> without normalization?

ok, i was not very precise in my statement. actually i was doing this:

1. make a plot of lhs pointwise.
2. make a plot of lhs - rhs pointwise, on the same scale.

then compare the size of the deviations lhs-rhs. what i see is that the
absolute deviations lhs-rhs do increase over time and that they reach a
level over 10% of the lhs, at those points in the phase space where the
lhs is largest in magnitude

sorry for the confusion.

n.


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 ]

Re: reaction-advection-diffusion; problems making the solution converge

[Daniel Wheeler]

On Wed, Aug 3, 2016 at 9:27 AM, Nils Becker
 wrote:
>> 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'm not sure if using ||lhs|| is a good way to normalize the residual.
Surely that could be very small. Is ||lhs - rhs|| getting smaller
without normalization?

-- 
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]

Re: reaction-advection-diffusion; problems making the solution converge

[Nils Becker]

> 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 

Re: Applying Bounding Limits to PDE Solution Variables

[Daniel Wheeler]

On Tue, Aug 2, 2016 at 2:02 PM, Campbell, Ian
 wrote:
>
> I've tried the mean-square and RMS values, but they're both orders of 
> magnitude away from the residual returned by '.sweep'. As such, can you 
> please tell me the formula used by '.sweep' to compute the residual, so that 
> I may use the same one?

I think it's the L2 norm of "A * x - b" without any normalization,
where x is the variable value before the linear solve. See,

https://github.com/usnistgov/fipy/blob/develop/fipy/solvers/solver.py#L135

-- 
Daniel Wheeler
___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]