Using the solution of a function from a previous time step

On 08/02/11 13:41, David Estrella wrote:
> New question #166719 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/166719
>
> Hi Everyone,
>
> I'm a bit of a newcomer to FEniCS, and I'm working on an advection-dispersion-reaction groundwater solute transport model. In this model there is a sorption term (S()) that calculates the total amount of removed solute during a time step. What I need to do is save this value as S_prev and use it in the next time step so that the removed amount accumulates over time. It is clear to me how to do this with the solution of the variational problem (C), but I am having trouble doing the same thing with a simple user defined function. I believe that I am not correctly saving S_prev as a value but as the function S().

Correct, you could try doing something of the form

   S_current = project(S(), V)
   S_prev.assign(S_current)

inside your loop.

--
Marie

> Then the system is unable to calculate the solution. Here is my code:
>
> mesh = Rectangle(0, 0, 2, 1, 100, 1, 'crossed')
> V = FunctionSpace(mesh, "CG", 1)
>
> # Define Dirichlet boundary
> def Dirichlet_boundary1(x):
> return x[0]< DOLFIN_EPS
>
> dt = 60
> T = 14400
> t = dt
>
> C1 = Constant(1.5)
> bc1 = DirichletBC(V, C1, Dirichlet_boundary1)
>
> C0 = Constant(0)
> C0.t = 0
> C_prev = interpolate(C0, V)
> C_prev2 = interpolate(C_prev, V)
>
> S0 = Constant(0)
> S0.t = 0
> S_prev = interpolate(S0, V)
>
>
> # Define variational problem
>
> C = TrialFunction(V)
> w = TestFunction(V)
> S = Function
> f = Constant(0) #Flux() # Source
> alphdisp=Constant(0.0051)
> q = Constant(8.19e-5) #Darcy Velocity
> De = Constant(2.807e-7)
> D = Expression('alphdisp*vs+De') # Dispersion Coefficient
> #D = Constant(0)
> k = Constant(5.53e-5) # Constant removal term
> #k = Constant(0)
>
>
> n = Constant(0.34) # Porosity
> vs = Constant(3.35e-4) # Seepage Velocity
> Smax = Constant(1.999e-3) # Maximum possible Sorption
> pb = Constant(2531) # Density kg/m3
> den = 1658
>
> def katt():
> return (1-S_prev/(Smax+n/den*k*C_prev))
>
> def S():
> return S_prev+(n/den*(katt()+k)*C_prev*dt)
>
>
> a = C*w*dx + dt*w*C.dx(0)*vs*dx + dt*D*inner(grad(C),grad(w))*dx +k*dt*C*w*dx + (1-S()/(Smax+(n/den*k*C_prev)))*dt*C*w*dx
> L = (C_prev + dt*f)*w*dx
>
> C = Function(V)
>
> A = assemble(a)
>
> # Save solution in VTK format
> file = File("transportsolution.pvd");
>
> while t<= T:
>
> b = assemble(L)
> C0.t = t
> bc1.apply(A,b)
> solve(A,C.vector(),b)
>
> t += dt
> S_prev=S()
> C_prev2.assign(C_prev)
> C_prev.assign(C)
> file<< (C)
>
>
> # Plot solution
> plot(C, interactive=True)
> plot(S(), interactive=True)
>
>
> The program begins to run but then slows down before giving this error
>
> File "trans.py", line 90, in<module>
> S_prev=S()
> File "trans.py", line 69, in S
> return S_prev+(n/den*(katt()+k)*C_prev*dt)
> File "trans.py", line 66, in katt
> return (1-S_prev/(Smax+n/den*k*C_prev))
> File "/usr/lib/python2.6/dist-packages/ufl/exproperators.py", line 134, in _rsub
> return Sum(o, -self)
> File "/usr/lib/python2.6/dist-packages/ufl/exproperators.py", line 168, in _neg
> return -1*self
> File "/usr/lib/python2.6/dist-packages/ufl/exproperators.py", line 110, in _rmul
> return _mult(o, self)
> File "/usr/lib/python2.6/dist-packages/ufl/exproperators.py", line 90, in _mult
> p = Product(a, b)
> File "/usr/lib/python2.6/dist-packages/ufl/algebra.py", line 231, in __new__
> self._init(*operands)
> File "/usr/lib/python2.6/dist-packages/ufl/algebra.py", line 243, in _init
> self._repr = "Product(%s)" % ", ".join(repr(o) for o in self._operands)
> MemoryError
>
> Thanks in advance for the help.
>
> Best,
> David Estrella
>
>

Thanks for the answer, and sorry for the late response. I was away last week without my computer.

I tried to implement this into the code, but it still doesnt seem to be projecting the values correctly. S_current and S_prev are just returned as 0. Here are the updated parts of the code I used:
...
def S():
 return S_prev+(n/den*(katt()+k)*C_prev*dt)
S_current = project(S(), V)
...
t += dt
 S_prev.assign(S_current)
 C_prev2.assign(C_prev)
 C_prev.assign(C)
 file << (C)

Perhaps I am not properly understanding how 'project' works. Thanks again in advance.

I think I may have made this question more complicated than it is, so I will try to rephrase in a simpler way. All I really need to do is store the value of a user defined function and assign this value to be used during the next time step. All of the values of this function are known, but when I try to store the value it just stores the function rather than the value. From what I understand project works only for the solution of the pde, which I believe is why it did not work above.

On 08/23/11 11:05, David Estrella wrote:
> Question #166719 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/166719
>
> David Estrella gave more information on the question:
> I think I may have made this question more complicated than it is, so I
> will try to rephrase in a simpler way. All I really need to do is store
> the value of a user defined function and assign this value to be used
> during the next time step. All of the values of this function are
> known, but when I try to store the value it just stores the function
> rather than the value. From what I understand project works only for
> the solution of the pde, which I believe is why it did not work above.
>

Could you try really stripping down your original code to provide
a very simple (but runnable!) example of the idea of what you want to do?
(A single simple Python function, a single dolfin Function, and a simple
loop
would probably suffice.) That would very much lower the effort level
involved for taking a look at it.

Also, project operates on dolfin Function (which the solution of a pde
typically
is an example of).

--
Marie

Thanks, Marie, for the quick response. I'm going to try to reply in two parts, first with just the simplest explanation, and again with a more stripped down code. I'm working on stripping it down, but I'm not sure if it will still be too much, so first this is the simplest way that might be enough to avoid looking more deeply into the code.

I would like to take a user-defined function:

def S():
   return S_prev+(k*C_prev*dt)

all of these values are known, with C_prev being the solution of the pde from the previous time step stored already. I basically need to calculate this as a value at each point in the domain and save the values as S_prev to be used in the next time step. The idea is that this value accumulates over time. Maybe I need to use an expression rather than defning a function.

S_prev=S()
only saves the function and doesnt work in the more complicated model that I have

S_current=project(S(), V)
returns 0

I would also like to be able to store these values like the pde solution for analysis after the program runs

I hope this is enough, but I will still work on stripping the whole code down. Sorry to make things so complicated.

Thanks again,
Dave

I have got it working now! Thank you very much!

Thanks Marie Rognes, that solved my question.

Great!

--
Marie

On 23. aug. 2011, at 17:55, David Estrella <email address hidden> wrote:

> Question #166719 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/166719
>
> Status: Answered => Solved
>
> David Estrella confirmed that the question is solved:
> I have got it working now! Thank you very much!
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.