Problem with Maximum Removal

After looking into this a little more, I've recognized the problem that the term lang is just using the initial value, whereas the other functions are indeed calculated correctly. Perhaps an easier question to answer is: how can I update the values from the function langmuir to be used in the next time step of the pde?

On Thursday September 22 2011 10:25:57 David Estrella wrote:
> Question #171956 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171956
>
> David Estrella gave more information on the question:
> After looking into this a little more, I've recognized the problem that
> the term lang is just using the initial value, whereas the other
> functions are indeed calculated correctly. Perhaps an easier question
> to answer is: how can I update the values from the function langmuir to
> be used in the next time step of the pde?

Much easier!

  u_prev.vector()[:] = u.vector()

Johan

I haven't looked in detail at what you want to do, but maybe you're asking the wrong question. If you want to have max(a, b) in your variational formulation, use conditional(lt(a,b), b, a) in your form. That will limit the expression at each quadrature point when assembling.

Martin

Den 22. sep. 2011 kl. 19:25 skrev David Estrella <email address hidden>:

> Question #171956 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171956
>
> David Estrella gave more information on the question:
> After looking into this a little more, I've recognized the problem that
> the term lang is just using the initial value, whereas the other
> functions are indeed calculated correctly. Perhaps an easier question
> to answer is: how can I update the values from the function langmuir to
> be used in the next time step of the pde?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

I think this conditional is will solve at least one of the problems I'm having, but I can't quite get it to work. Is it possible to use functions within the conditional? What I would need to do looks something like:

def S():
   return S_prev+(k2*lang)*C*dt
 S=project(S(), V)

Sc=conditional(lt(S,Smax), S, Smax)

a = ..........+ k2*(1-Sc/Smax)*dt*w*C*dx

so when Sc is greater than Smax it is set equal to Smax, otherwise it is evaluated with the defined function.

I am getting a rather long error message that says the conditional is not supported. Have I written this incorrectly?
Thanks for the help,
Dave

sorry, I didn't mean to click problem solved...

Is it maybe possible to use an expression like:

lang=Expression("S>Smax ? 0 : (1- S/Smax)")

instead of the conditional statement? Neither one is working for me, but I'm not sure exactly where is it supported to use functions or expressions.

Dave

I just needed to update my version of fenics, and now the conditional statement works fine.

I still however, have the problem that I need to update the coefficients in a and assemble the stiffness matrix at each timestep. I'm tying solve a function at each time step and use the coefficients at the next, but it appears to be just taking the initial value (for lang).

a =.....+ k2*lang*dt*w*C*dx
C=Function(V)

while t <= T:
            A = assemble(a)
            b = assemble(L)
            bc1.apply(A,b)
            solve(A,C.vector(),b)

            t += dt
            C_prev.assign(C)
           S_prev.assign(S)

         def langmuir():
                 return 1-(S_prev/Smax)
         lang=project(langmuir(), V)

Is there any way to deal with something like this?
Thanks again for the help,
Dave

You must do lang.assign(value) as you do with the other coefficients.

The Form object 'a' keeps a reference to the Function object 'lang'
you use to construct 'a', and when you do 'lang = ...' you create a
new python object which the Form object 'a' has no knowledge of.

Martin

On 25 September 2011 14:45, David Estrella
<email address hidden> wrote:
> Question #171956 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171956
>
> David Estrella gave more information on the question:
> I just needed to update my version of fenics, and now the conditional
> statement works fine.
>
> I still however, have the problem that I need to update the coefficients
> in a and assemble the stiffness matrix at each timestep.   I'm tying
> solve a function at each time step and use the coefficients at the next,
> but it appears to be just taking the initial value (for lang).
>
> a =.....+ k2*lang*dt*w*C*dx
> C=Function(V)
>
> while t <= T:
>            A = assemble(a)
>            b = assemble(L)
>            bc1.apply(A,b)
>            solve(A,C.vector(),b)
>
>            t += dt
>            C_prev.assign(C)
>           S_prev.assign(S)
>
>         def langmuir():
>                 return 1-(S_prev/Smax)
>         lang=project(langmuir(), V)
>
> Is there any way to deal with something like this?
> Thanks again for the help,
> Dave
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Hi Martin, thanks for the answer. This doesn't seem to work for me. Somehow I'm still getting the solution calculated only with the initial value of lang. Should I be somehow using the function rather than the projection?
Dave

On 25 September 2011 19:45, David Estrella
<email address hidden> wrote:
> Question #171956 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171956
>
>    Status: Answered => Open
>
> David Estrella is still having a problem:
> Hi Martin, thanks for the answer.  This doesn't seem to work for me.

You should show us some code, I don't know what you tried to do now.

Martin

 Somehow I'm still getting the solution calculated only with the
initial value of lang.  Should I be somehow using the function rather
than the projection?
> Dave
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Here is the full code, with the recent changes. What I'm tryng to get is the removal term with lang in a to decrease and disappear as removal occurs, but as seen in the final plot I'm still getting the almost instant removal of C as would occur in the first time step.

from dolfin import *

# Create mesh and define function space
mesh = Rectangle(0, 0, 2, 1, 100, 1, "crossed")
V = FunctionSpace(mesh, "CG", 1)

# Define Dirichlet boundary (x = 0), neumann at (x=4)
def Dirichlet_boundary1(x):
 return x[0] < DOLFIN_EPS

dt = 300
T = 9000
t = dt

C1 = Constant(10)
bc1 = DirichletBC(V, C1, Dirichlet_boundary1)

C0 = Constant(0)
C_prev = interpolate(C0, V)

S=Constant(0)
S_prev = interpolate(S, V)

# Define variational problem

C = TrialFunction(V)
w = TestFunction(V)
v=Constant(5e-3) #velocity
S0max = Constant(1e-2) # Maximum possible Sorption
Smax=S0max
k1=Constant(5e-4)
k2=Constant(5e-2)

if t<=dt:
 lang=Constant(1)

a = C*w*dx + dt*v*w*C.dx(0)*dx + k1*dt*w*C*dx+ k2*lang_prev*dt*w*C*dx
L = C_prev*w*dx

C = Function(V)

# Save solution in VTK format
file = File("Ssolution.pvd");

while t <= T:
 A = assemble(a)
 b = assemble(L)
 bc1.apply(A,b)
 solve(A,C.vector(),b)

 t += dt
 C_prev.assign(C)
 S_prev.assign(S)
 lang_prev.assign(lang)

 def langmuir():
  return 1-(S_prev/Smax)
 lang=project(langmuir(), V)

 def S():
  return S_prev+((2.12e-4)*(k2*lang)*C*dt)
 S=project(S(), V)
 file<<(S)

 b = assemble(L, tensor=b)

plot(lang, title="Langmuir coefficient")
plot(S, title="S")
plot(C, title="C")
interactive()

On 25 September 2011 20:50, David Estrella
<email address hidden> wrote:
> Question #171956 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171956
>
>    Status: Answered => Open
>
> David Estrella is still having a problem:
> Here is the full code, with the recent changes.  What I'm tryng to get
> is the removal term with lang in a to decrease and disappear as removal
> occurs, but as seen in the final plot I'm still getting the almost
> instant removal of C as would occur in the first time step.
>
> from dolfin import *
>
> # Create mesh and define function space
> mesh = Rectangle(0, 0, 2, 1, 100, 1, "crossed")
> V = FunctionSpace(mesh, "CG", 1)
>
> # Define Dirichlet boundary (x = 0), neumann at (x=4)
> def Dirichlet_boundary1(x):
>        return x[0] < DOLFIN_EPS
>
> dt = 300
> T = 9000
> t = dt
>
> C1 = Constant(10)
> bc1 = DirichletBC(V, C1, Dirichlet_boundary1)
>
> C0 = Constant(0)
> C_prev = interpolate(C0, V)
>
> S=Constant(0)
> S_prev = interpolate(S, V)
>
> # Define variational problem
>
> C = TrialFunction(V)
> w = TestFunction(V)
> v=Constant(5e-3) #velocity
> S0max = Constant(1e-2) # Maximum possible Sorption
> Smax=S0max
> k1=Constant(5e-4)
> k2=Constant(5e-2)
>
> if t<=dt:
>        lang=Constant(1)
>
> a = C*w*dx + dt*v*w*C.dx(0)*dx + k1*dt*w*C*dx+ k2*lang_prev*dt*w*C*dx

This is not a runnable code, you're using lang_prev without defining
it. Also you didn't use assign on lang in your loop as I suggested in
a previous mail.

Martin

Hi Martin, I'm sorry about the non-running code. I was under the impression that something like this was what you had suggested with the lang_prev.assign(lang). If thats not what you meant, I'm afraid I didn't understand the suggestion. Thanks again for the help and quick responses. It is greatly appreciated.

Dave

"""This demo program solves

1D transient transport column problem
check page 54 of laboratory notebook for equations
"""

__author__ = "D-Munnie"
__date__ = "2011-07-15"
__copyright__ = ""
__license__ = "GNU LGPL Version 2.1"

# Begin model

from dolfin import *

# Create mesh and define function space
mesh = Rectangle(0, 0, 2, 1, 100, 1, "crossed")
V = FunctionSpace(mesh, "CG", 1)

# Define Dirichlet boundary (x = 0), neumann at (x=4)
def Dirichlet_boundary1(x):
 return x[0] < DOLFIN_EPS

dt = 300
T = 9000
t = dt

C1 = Constant(10)
bc1 = DirichletBC(V, C1, Dirichlet_boundary1)

C0 = Constant(0)
C_prev = interpolate(C0, V)

S=Constant(0)
S_prev = interpolate(S, V)

# Define variational problem

C = TrialFunction(V)
w = TestFunction(V)
v=Constant(5e-3) #velocity
S0max = Constant(1e-2) # Maximum possible Sorption
k1=Constant(5e-4)
k2=Constant(5e-2)

if t<=dt:
 lang=Constant(1)
 lang_prev=interpolate(lang, V)

def langmuir():
 return 1-(S_prev/S0max)
lang=project(langmuir(), V)
#
a = C*w*dx + dt*v*w*C.dx(0)*dx + k1*dt*w*C*dx+ k2*lang_prev*dt*w*C*dx
L = C_prev*w*dx

C = Function(V)

# Save solution in VTK format
file = File("Ssolution.pvd");

while t <= T:
 A = assemble(a)
 b = assemble(L)
 bc1.apply(A,b)
 solve(A,C.vector(),b)

 t += dt
 C_prev.assign(C)
 S_prev.assign(S)
 lang_prev.assign(lang)

 def S():
  return S_prev+((2.12e-4)*(k2*lang)*C*dt)
 S=project(S(), V)
 file<<(S)

 b = assemble(L, tensor=b)

plot(lang, title="Langmuir coefficient")
plot(S, title="S")
plot(C, title="C")
interactive()

Hi Martin, Thanks for the help! I got it to work. I got rid of the python functions used project in assign as follows:
Sc_prev.assign(project(Sc,V)).
Dave

Thanks Martin Sandve Alnæs, that solved my question.