Piecewise Function Question

sorry, i didnt finish the previous post....

k=0.001
Smax=0.0143
def S():
          return: S_prev+k*(1-S_prev/Smax)*C*dt
S=project(S(), V)
class Sc(Exression):
          def eval(self, values, x):
                        if S(x)>Smax:
                                values[0]=Smax
                        else:
                                values[0]=S(x)
Sc=Sc()

I recieve the error:
File "/usr/lib/python2.6/dist-packages/dolfin/fem/project.py", line 56, in project
    b = assemble(L)
  File "/usr/lib/python2.6/dist-packages/dolfin/fem/assemble.py", line 132, in assemble
    add_values)
  File "/usr/lib/python2.6/dist-packages/dolfin/cpp.py", line 9904, in assemble
    return _cpp.assemble(*args)
Exception: Swig director method error. Error detected when calling 'Expression.eval_data'

It looks like something goes wrong in the eval method. Try evaluating the
expression outside the eval function with dummy arguments (x, values).

You can also acomplish what you want by using the index interface of the
GenericVector:

  S.vector()[S.vector()>Smax] = Smax

Johan

On Tuesday September 20 2011 05:20:51 David Estrella wrote:
> Question #171764 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171764
>
> David Estrella gave more information on the question:
> sorry, i didnt finish the previous post....
>
> k=0.001
> Smax=0.0143
> def S():
> return: S_prev+k*(1-S_prev/Smax)*C*dt
> S=project(S(), V)
> class Sc(Exression):
> def eval(self, values, x):
> if S(x)>Smax:
> values[0]=Smax
> else:
> values[0]=S(x)
> Sc=Sc()
>
> I recieve the error:
> File "/usr/lib/python2.6/dist-packages/dolfin/fem/project.py", line 56, in
> project b = assemble(L)
> File "/usr/lib/python2.6/dist-packages/dolfin/fem/assemble.py", line 132,
> in assemble add_values)
> File "/usr/lib/python2.6/dist-packages/dolfin/cpp.py", line 9904, in
> assemble return _cpp.assemble(*args)
> Exception: Swig director method error. Error detected when calling
> 'Expression.eval_data'

Thanks for the quick answer. I tried using:

def S():
          return: S_prev+k*(1-S_prev/Smax)*C*dt
S=project(S(), V)
S.vector()[S.vector()>Smax] = Smax

but I get this error:

 S.vector()[S.vector()>Smax]=Smax
  File "/usr/lib/python2.6/dist-packages/dolfin/cpp.py", line 2345, in __setitem__
    raise TypeError, "provide a scalar to set single item"
TypeError: provide a scalar to set single item

This is with Smax defined as a constant. Should it be interpolated or projected onto the mesh?

On Tuesday September 20 2011 09:05:57 David Estrella wrote:
> Question #171764 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171764
>
> Status: Answered => Open
>
> David Estrella is still having a problem:
> Thanks for the quick answer. I tried using:
>
> def S():
> return: S_prev+k*(1-S_prev/Smax)*C*dt
> S=project(S(), V)
> S.vector()[S.vector()>Smax] = Smax
>
> but I get this error:
>
> S.vector()[S.vector()>Smax]=Smax
> File "/usr/lib/python2.6/dist-packages/dolfin/cpp.py", line 2345, in
> __setitem__ raise TypeError, "provide a scalar to set single item"
> TypeError: provide a scalar to set single item
>
> This is with Smax defined as a constant. Should it be interpolated or
> projected onto the mesh?

Then pick a scalar instead of the Constant ;)

Or turn the Constant into a scalar:

  Smax = float(Smax)

but that sounds a bit tedious as Smax has to be defined from a scalar.

Johan

Well that wasn't such a complicated solution :)

Thanks a lot for the help!

Thanks Johan Hake, that solved my question.

Hi Johan, I have a related question. I should probably start a new thread for this, but I'm having trouble simplifying my code while still getting the point across. Maybe there is a simple answer to what im looking for.
I am using the ratio of S/Smax in my pde as part of a removal term in my pde, so it looks something like:

if t<dt:
        S=Constant(0)

k=0.001
a = C*w*dx + dt*w*vs*C.dx(0)*dx + k*(1-S/Smax)*C*dt
L = C_prev*w*dx

and I have S here defined after the pde:
def S():
          return: S_prev+k*(1-S_prev/Smax)*C*dt
S=project(S(), V)
 S.vector()[S.vector()>Smax]=Smax

The problem is that when S reaches Smax, the value stops increasing, but it still removes from C at the increased rate. When plotting S it stops at the max value, but then in the end my system is no longer mass conservative because there is some quantity of C removed that is then not accounted for because S is limited at the max.

Sorry that this is a bit confusing to explain. If it doesn't make sense to you, or there is not simple answer based on the given information, I'll start a new thread and try to put together a simple code. Thanks again for the help.

Dave

I am a bit confused and cannot give you a good answer. Maybe you can start a
new thread and try to explain your problem in more details? Maybe also choose
variable names which are not overlapping with function names too?

I am also not sure how you could ever have mass conservation when you enforce
a maximal value on your solution by just reducing some values. You should be
able to enforce a maximal value in your PDE.

Johan

On Wednesday September 21 2011 03:05:48 David Estrella wrote:
> Question #171764 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171764
>
> Status: Solved => Open
>
> David Estrella is still having a problem:
> Hi Johan, I have a related question. I should probably start a new thread
> for this, but I'm having trouble simplifying my code while still getting
> the point across. Maybe there is a simple answer to what im looking for.
> I am using the ratio of S/Smax in my pde as part of a removal term in my
> pde, so it looks something like:
>
> if t<dt:
> S=Constant(0)
>
> k=0.001
> a = C*w*dx + dt*w*vs*C.dx(0)*dx + k*(1-S/Smax)*C*dt
> L = C_prev*w*dx
>
> and I have S here defined after the pde:
> def S():
> return: S_prev+k*(1-S_prev/Smax)*C*dt
> S=project(S(), V)
> S.vector()[S.vector()>Smax]=Smax
>
> The problem is that when S reaches Smax, the value stops increasing, but
> it still removes from C at the increased rate. When plotting S it stops
> at the max value, but then in the end my system is no longer mass
> conservative because there is some quantity of C removed that is then
> not accounted for because S is limited at the max.
>
> Sorry that this is a bit confusing to explain. If it doesn't make sense
> to you, or there is not simple answer based on the given information,
> I'll start a new thread and try to put together a simple code. Thanks
> again for the help.
>
> Dave

Alright, I'm working on making a simplified code now. Sorry for the confusing question, and thanks for the help.

The idea is that I want the removal term to disappear when it has reached the maximum value. I'll move this to a new question.... Thanks again!
Dave