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

[Nils Becker]

hi daniel, thanks for looking into this!


> Hi Nils,
> 
> It seems like the code is converging well as far as I can tell. For
> example, when I modify the sweep loop to be
> 
> 
> print
> print 'step: {0}'.format(tt)
> while res > 1e-9: #Solve.sweeps:
> res = eq.sweep(var=rho, dt=Solve.step,
> solver=Solve.customSolver)
> print res
> 
> 
> every step seems to converge in 5 steps to 1e-10.

ok. it seems that you have not changed the iteration, only added print
functions. so my code seems to be converging, yes?

that's also what i thought. however testing the solution externally, it
seems that the solution is not quite accurate (as time progresses, it
seems to get worse).

of course, this external test could be buggy -- however the fact that it
shows that lhs = rhs quite accurately in the beginning makes me think
that it could be correct.

> The non-linear residual only really tells you how well each step
> individually is converging. It seems like you were only using the very
> first non-linear residual, which uses the old value for the
> calculation, which doesn't tell you anything useful. 

so what sweep returns is basically A(n) * x(n-1) - b(n), yes? as far as
i can tell, i always collect the last one of these values in each step.
they do get small, i also see the convergence to below 1e-9 when i set
the tolerance that low.

> 
> One thing that I'm confused about is that it requires 5 sweeps to
> converge for a fully linear equation. Is it actually linear? I
> couldn't tell from the code.

yes. this also surprised me. the potential phi and the growth rate
lambda_ are of course non-linear functions of space. but all
coefficients are independent of the density, and the density rho should
be really appearing only linearly in the equation. so unless i
constructed the equation or the loop wrong, the system should really be
linear. no idea why it would require sweeps. do linear equations with an
implicit source term require sweeping? am i doing the sweep wrong?

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: Model convection with non-linear adsorption

[Daniel Wheeler]

Hi Martin,

The old values don't appear to be being updated. Is that correct? Don't you
need


...
while t < t_tot:
c_b.updateOld()
c_f.updateOld()
q_b.updateOld()
# Solve iteratively with sweep
while (residual > max_residual_nonlin) & (n_iter_nonlin <
n_max_allowed_iter):
...


Also, you can use an adaptive time step that changes based on whether the
non-linear tolerance is achieved for the given step.

Cheers,

Daniel


On Fri, Aug 5, 2016 at 9:46 AM, Korevaar, Martin <
martin.korev...@kwrwater.nl> wrote:

> Dear mailinglist,
>
>
>
> I am trying to solve a convection problem where the convected substance is
> adsorbed (on active carbon in this case). It is modeled as follows: the
> fluid is divided in a bulk part and in a film layer part; that is the film
> layer surrounding the adsorbing medium (active carbon). To the
> concentration in the bulk applies a convective equation with a sink for the
> part of the solute that is diffusing into the film layer; the sink depends
> on the difference between film and bulk concentration. The film layer thus
> receives the solute from the bulk and loses it to the adsorbing medium and
> is therefore modeled as an equation with only a sink and a source; the
> source depends on the concentration difference between bulk and film, the
> sink depends on the ‘concentration’ difference between the surface and the
> center of active carbon grains. Finally, I also modeled the diffusion of
> the solute from the surface of the active carbon to its center; this source
> term depends on the difference between solution at the surface and center
> of the grain.
>
>
>
> c_b: bulk concentration
>
> c_f: film concentration
>
> q_s: concentration at the surface of the grain (this is not really a
> concentration, but the amount of adsorbed material per amount of adsorbing
> material (active carbon)
>
> q_b: the concentration in the center of the grain
>
>
>
> This then yields the following equations:
>
> dc_b/dt = v dc_b/dx – A1(c_b – c_f)
>
> dc_f/dt = A2 (c_b – c_f) – A3(q_s – q_b)
>
> dq_s/dt = A4(q_s – q_b)
>
>
>
> here A1, A2, A3 and A4 are some mass transfer constants. Furthermore, q_s
> relates to c_f by the Freundlich equation:
>
> q_s = K_Fr * c_f ^n_Fr
>
>
>
> where K_Fr and n_Fr are the (Freundlich) constants. Implementation is
> attached as minimal working example.
>
>
>
> My problem is that this only converges at very small timesteps, that is of
> order of seconds, while I need to simulate order of years. I need larger
> timesteps, but I do not know how to improve this further. The changes in
> time are not so large, so it should be possible. Apparently this is some
> tricky numerical issue. I tried all the different solvers in FiPy and it
> seemed that the most stable solution occurs using LinearLUSolver with zero
> gradient outlet boundary condition at the outlet. I linearized the
> non-linear isotherm which increased my allowed timestep from milliseconds
> to seconds. Can you help me improve my simulation speed? I hope you have
> sufficient information with the implementation and value of the constants
> as attached and given below. If not, please let me know.
>
>
>
> Value for the constants:
>
>
>
> A1 = 0.36
>
> A2 = 0.48
>
> A3 = 0.0148
>
> A4 = 1.e-8
>
> v = 1.6e-3  # m/s
>
> K_Fr = 386
>
> N_Fr = 0.33
>
> Best regards,
>
> M.W. (Martin) Korevaar, MSc
> Scientific researcher - Drinking Water Treatment *|* KWR Watercycle
> Research Institute *|* Groningenhaven 7, P.O. Box 1072, 3430 BB
> Nieuwegein, the Netherlands [image:
> http://www.kwrwater.nl/uploadedImages/mail/gmaps.png]
> 
>  *| T* +31 30 606 9515 *| M *+31 6 11648156 *| E*
> martin.korev...@kwrwater.nl *| W* www.kwrwater.nl* |* Follow KWR on [image:
> http://www.kwrwater.nl/uploadedImages/mail/twitter.png]
> *| *Follow me on  [image:
> http://www.kwrwater.nl/uploadedImages/mail/linkedin.png]
>  *|* Chamber of
> Commerce Utrecht e.o. 27279653.
>
> The KWR office building 
> has won the 'Most Stimulating Environment' award of the Royal Institute of
> Dutch Architects.
>
>
>
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
>


-- 
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

[Daniel Wheeler]

Hi Nils,

It seems like the code is converging well as far as I can tell. For
example, when I modify the sweep loop to be


print
print 'step: {0}'.format(tt)
while res > 1e-9: #Solve.sweeps:
res = eq.sweep(var=rho, dt=Solve.step,
solver=Solve.customSolver)
print res


every step seems to converge in 5 steps to 1e-10.

step: 20.0
0.00192315625797
4.5626520206e-05
1.08935792171e-06
2.60975675747e-08
6.26480480874e-10

step: 21.0
0.00195977695607
4.651749367e-05
1.11103436368e-06
2.66255569394e-08
6.3934852818e-10

step: 22.0
0.00199768415442
4.74384979099e-05
1.13342888153e-06
2.71708361386e-08
6.52634229319e-10

The non-linear residual only really tells you how well each step
individually is converging. It seems like you were only using the very
first non-linear residual, which uses the old value for the
calculation, which doesn't tell you anything useful. Basically "A_new
* x_old -b_new" is not a measure of anything useful. It is only useful
as a quantity for normalization. FiPy uses the old value or the
previous sweep value to calculate the non-linear residual. So the
residual is

   A(n) * x(n-1) - b(n)

where n is the sweep number, and x_(-1) is x_old if the first sweep is
n=0. If only one sweep is executed then the residual isn't a useful
quantity. The linear residual is not returned from the "sweep" method,
which is "A(n) * x(n) - b(n)" and is assumed to be as accurate as the
tolerance of the linear solve. Basically, collecting the first
non-linear residual at each time step is not a useful exercise to
generate a metric for convergence across time steps. In fact, I'm not
even sure how to do that or what it means. I only know how to
understand convergence for a given time step.

One thing that I'm confused about is that it requires 5 sweeps to
converge for a fully linear equation. Is it actually linear? I
couldn't tell from the code.

Cheers,

Daniel




On Fri, Aug 5, 2016 at 8:30 AM, Nils Becker
 wrote:
>>> 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?
>
> if you don't see any obvious errors in building the solution in fipy, i
> could continue poking around if i knew how to evaluate the various terms
> of the solution and their errors at different times. i see methods
> eq_diffprolif.justErrorvector, but i don't quite understand what exactly
> this returns. is it lhs-rhs? what is justResidualVector? can i get those
> for the individual terms on the rhs as well?
>
> cheers, n.
>
> ___
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
>   [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>



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

Model convection with non-linear adsorption

[Korevaar, Martin]

Dear mailinglist,

I am trying to solve a convection problem where the convected substance is 
adsorbed (on active carbon in this case). It is modeled as follows: the fluid 
is divided in a bulk part and in a film layer part; that is the film layer 
surrounding the adsorbing medium (active carbon). To the concentration in the 
bulk applies a convective equation with a sink for the part of the solute that 
is diffusing into the film layer; the sink depends on the difference between 
film and bulk concentration. The film layer thus receives the solute from the 
bulk and loses it to the adsorbing medium and is therefore modeled as an 
equation with only a sink and a source; the source depends on the concentration 
difference between bulk and film, the sink depends on the 'concentration' 
difference between the surface and the center of active carbon grains. Finally, 
I also modeled the diffusion of the solute from the surface of the active 
carbon to its center; this source term depends on the difference between 
solution at the surface and center of the grain.

c_b: bulk concentration
c_f: film concentration
q_s: concentration at the surface of the grain (this is not really a 
concentration, but the amount of adsorbed material per amount of adsorbing 
material (active carbon)
q_b: the concentration in the center of the grain

This then yields the following equations:
dc_b/dt = v dc_b/dx - A1(c_b - c_f)
dc_f/dt = A2 (c_b - c_f) - A3(q_s - q_b)
dq_s/dt = A4(q_s - q_b)

here A1, A2, A3 and A4 are some mass transfer constants. Furthermore, q_s 
relates to c_f by the Freundlich equation:
q_s = K_Fr * c_f ^n_Fr

where K_Fr and n_Fr are the (Freundlich) constants. Implementation is attached 
as minimal working example.

My problem is that this only converges at very small timesteps, that is of 
order of seconds, while I need to simulate order of years. I need larger 
timesteps, but I do not know how to improve this further. The changes in time 
are not so large, so it should be possible. Apparently this is some tricky 
numerical issue. I tried all the different solvers in FiPy and it seemed that 
the most stable solution occurs using LinearLUSolver with zero gradient outlet 
boundary condition at the outlet. I linearized the non-linear isotherm which 
increased my allowed timestep from milliseconds to seconds. Can you help me 
improve my simulation speed? I hope you have sufficient information with the 
implementation and value of the constants as attached and given below. If not, 
please let me know.

Value for the constants:

A1 = 0.36
A2 = 0.48
A3 = 0.0148
A4 = 1.e-8
v = 1.6e-3  # m/s
K_Fr = 386
N_Fr = 0.33
Best regards,
M.W. (Martin) Korevaar, MSc
Scientific researcher - Drinking Water Treatment | KWR Watercycle Research 
Institute | Groningenhaven 7, P.O. Box 1072, 3430 BB Nieuwegein, the 
Netherlands [http://www.kwrwater.nl/uploadedImages/mail/gmaps.png] 

  | T +31 30 606 9515 | M +31 6 11648156 | E 
martin.korev...@kwrwater.nl | W 
www.kwrwater.nl | Follow KWR on 
[http://www.kwrwater.nl/uploadedImages/mail/twitter.png] 
 | Follow me on  
[http://www.kwrwater.nl/uploadedImages/mail/linkedin.png] 
  | Chamber of Commerce Utrecht e.o. 
27279653.

The KWR office building has 
won the 'Most Stimulating Environment' award of the Royal Institute of Dutch 
Architects.



MWE.py
Description: MWE.py
___
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]

>> 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?

if you don't see any obvious errors in building the solution in fipy, i
could continue poking around if i knew how to evaluate the various terms
of the solution and their errors at different times. i see methods
eq_diffprolif.justErrorvector, but i don't quite understand what exactly
this returns. is it lhs-rhs? what is justResidualVector? can i get those
for the individual terms on the rhs as well?

cheers, 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 ]