Solution abnormalities

What is the unknown in your PDE? concentartion?

Yes.

I am pretty much using the same example given in the tutorials except I converted to 3D and vary my f function over space

Are you sure that you set the concentration to 300 on the edges? You mentioned that the maximum value of the solution is 0.

What is the sign of f?

Yes I check the initial conditions that I save in the output and they are set correctly. However, I think the problem with my current approach is that I was solving delC = f(x,y,z) and I want to be solving this equation D*delC = f(x,y,z)*c (where delC is supposed represent that del operator with respect to concentration).

How would I set the equation up in escript? Are there any examples which solve similar problems?

Does "delC = f(x,y,z) " mean "C_{,i} = f_i(x,y,z)" ?
What are D and c in "D*delC = f(x,y,z)*c "?

D is the diffusion coefficient and c is the same as C, sorry capitalization error. I am not sure about the other terminology but here is the full equation:

D * (d2C/dx2 + d2C/dy2 + d2C/dz2) = f(x,y,z)*C(x,y,z)

Were d2C/dx2 represents the second derivative of C with respect to x for example.

I would recommend to use the Helmholtz class, see http://esys.esscc.uq.edu.au/esys13/nightly/epydoc/esys.escript.linearPDEs.Helmholtz-class.html. For your case

from esys.escript.lienarPDEs import Helmholtz
my_pde=Helmholtz(domain)
my_pde.setValue(k=D, omega=f)
C=my_pde.getSolution()

if f>0 no constraints on C are required.

Thank you so much for your help! Is there any way to move the D out of the first grad term, as the diffusion coefficient does not vary with space?

omega*u - grad(k*grad(u)[j])[j] = f

That is what this notation dictates correct? Also if my f is zero then I guess I will need to place some sort of constraint on the solution correct?

If D does not vary in space it does not matter if you write D in front of the first grad or between the two grads.
by default the boundary condition "n.grad(u)=0" is applied (n is the surface normal".

Correct, if f is zero you need to add a constraint.

NB: the notation in "grad(k*grad(u)[j])[j]" is not very clear. It actually means "div(k*grad(u))" or "(ku_{,i}_{,i}". Again
div(k*grad(u))=k * laplace(u) if k is constant.

So I looked at the example problems using Helmholtz equation and wanted to make sure I am understanding how to create boundary conditions. If I want to set f to zero, I can call the function .setValue(f = 0). I already know how to set the values for omega and k. However, again I have a box where the values at the edge should be set to some value greater than 0. To accomplish this do I need to set the g or alpha parameter? I am assuming that I need to set G and assume that the alpha is zero... although maybe I misunderstand?

Conversely if f = 0, I will need to apply a constraint. In this case the parameters for constraints are r and q, but I am not quite sure what they mean or how to set them. I am assuming I must make r equal to a small number or zero so that the simulation cant be negative for a given positive q?

yes, by default alpha and g are assumed to be zero.

and

yes, if f =0 you need to set constraints. r sets the value of the solution and q sets the location where constraints are applied.
The following code snippet sets the solution to v[i] on the faces x[i]=0 and x[i]=L[i] (i=0,1,2):

x=domain.getX()
my_q=0
my_r=0
for i in range(DIM):
    m = whereZero( x[i] ) + whereZero( x[i] - L[i] )
    my_q=my_q+m
    my_r = (1-m) * my_r + m * v[i]

pde.setValue(r=my_r, q=my_q)

Ok thanks for the explnantion. I think I undertsand how this works now. You need to use q to define where to define the boundary conditions are applied, and use can use R to define the value of the solution at those points. I used a my own script to define the edges of the box where I want to apply constraints:

for x in xrange(0, 21):
    for y in xrange(0, 21):
        i1 = pointToIndex(x,y,0)
        i2 = pointToIndex(x,y,20)
        gammaD.setValueOfDataPoint(i1, 1)
        gammaD.setValueOfDataPoint(i2, 1)

for x in xrange(0, 21):
    for z in xrange(0, 21):
        i1 = pointToIndex(x,0,z)
        i2 = pointToIndex(x,20, z)
        gammaD.setValueOfDataPoint(i1, 1)
        gammaD.setValueOfDataPoint(i2, 1)

for y in xrange(0, 21):
    for z in xrange(0, 21):
        i1 = pointToIndex(0,y,z)
        i2 = pointToIndex(20,y, z)
        gammaD.setValueOfDataPoint(i1, 1)
        gammaD.setValueOfDataPoint(i2, 1)

def pointToIndex(x, y, z):
    return int(x + (21)*y + (21)*(21)*z)

Where in this case the length of my brick was 20 by 20 by 20. I then used the following code to define r and q:

q = gammaD
r = gammaD*100

Which sets the values at the border of my solution to 100. I then used this code to generate a random omega defined throughout space:

f0 = Scalar(0, d.getX().getFunctionSpace())
alpha = 2
for i in xrange(21*21*21):
    f0.setValueOfDataPoint(i, numpy.random.randint(-1,2)*alpha)

And then I graph and visualize using the visit cli. Thanks for your help! I did want to note one interesting feature of this system is that it appears to become unstable if the w function is allowed to vary wildly over the larger ranges. For example in this code if alpha is allowed to increase above 5, the solution then becomes negative over the range, which is not possible.

Again thank you for all of your help!

Ok so I have to add one more thing to this problem in case anyone else ever asks for a similar application. The source functions which are independent of concentration need to be defined in a separate function which is set to f. Only decay terms dependent on concentration should be included in the omega definition. Again, thanks for all of your help, I really appreciate it!

I would like to emphasize that in the way you are using setValueOfDataPoint to set gammaD you make assumptions about the ordering of the data points. This is dangerous and will not work if you try to run your program under MPI or use the optimize option in the Brick function. Therefore we strongly discourage people to use setValueOfDataPoint. As I showed #12 you don't need to use setValueOfDataPoint for setting gammaD.

In fact, things go wrong when you set f0 as f0 has in fact 20*20*20 *8 data points (not 21*21*21) as f0 is hold on the quadrature points of elements. So only an 1/8 of the array representing f0 is filled by random numbers the rest by zero. This can be easily be fixed by replacing "xrange(21*21*21)" by "xrange(f0.getNumDataPoints()".

NB: I am not so sure why you expect C to be positive regardless of what the value for alpha is.

The problem I had with 12 is that it produced summations where the different edges of the box overlapped. To avoid this I wrote the other script. I am confused about the number of data points though, when you ask the structure to print itself, it says it has 9261 data points, which is 21X21X21, and when use the command you mentioned above it also returns this number. Are you saying that the setValueofDataPoint() only sets one of the eight values associated with this data point, and thus should not be used in this manner? How would you recommend that I set individual values in the functions?

Ops - missed a step. omega is assumed to be presented at numerical integration points but your f0 is created at grid points - so your 9261 data points is correct. There is an interpolation of f0 happening behind the scene. Sorry - all good for f0.

For q it is only important that the value is positive the actual value is not important. For r this is a different story but has r is contant you can just say "pde.setValue(r=100)".

I see, so a positive q defines the locations to apply the boundary condition whose value is specified by r, which essentially allow you to specify the Dirichlet boundary conditions. If I want to specify a Neumann boundary condition then I need to use alpha and g parameters which allow you to control the flux. For this you said the value of alpha matters, so should I specify a constant value of alpha over the edge domain and zero otherwise, and then define a value for g which will control the actual value of the flux at that point? Or do i need to set alpha to zero and then define g only on the edges of the box?

alpha and g have the property "FunctionOnBoundary" so they are defined on the boundary only and they are set on the integration points with each surface element. If you want top set the normal flux to 1 on all boundaries (by default alpha and g =0):

pde.setValue(g=1, alpha=0)

if you want to set different values on different faces:

g=Scalar(0., FunctionOnBoundary(d))
g.setTaggedValue("right", 1.)
pde.setValue(g=1)

will set the normal flux to 1 on flur the "right" face (x[0]=0).