Hi, 

I`m trying to solve the 1D advection-diffusion-reaction equation dc/dt+u*dc/dx=D*dc2/dx2-kC using Fortan code but I`m still facing some issues. first I solved the advection-diffusion equation without including the source term (reaction) and it works fine. but when including the source term (decay of substence with the fisr order decay -kC)I could not get a correct solution. I have read that I need to use the spliting time mthod or as called also the operator time method, but I do not know how to apply it.
Could any one help to provide a simple code for solving this advection-diffusion-reaction equation ? I would appriciate your help.

Regards,

Kader

kader <fabdelkader@pi.ac.ae> wrote:
 
> I`m trying to solve the 1D advection-diffusion-reaction equation 
> dc/dt+u*dc/dx=D*dc2/dx2-kC using Fortan code but I`m still 
> facing some issues. 

See: https://class.coursera.org/scientificcomp-003

he does it in Matlab, and it also works in Octave. 

Actually, that is in 2D, which makes much nicer pictures.
You can cheat and go directly to lecture 19, 20, or 21.

> first I solved the advection-diffusion equation without 
> including the source term (reaction) and it works fine. 
> but when including the source term (decay of substence with 
> the fisr order decay -kC)I could not get a correct solution. 
> I have read that I need to use the spliting time mthod or 
> as called also the operator time method, but I do not know 
> how to apply it.

Are you doing it with Fourier transforms?

-- glen

On 9/20/2014 5:10 AM, glen herrmannsfeldt wrote:
> kader<fabdelkader@pi.ac.ae>  wrote:
>
>> I`m trying to solve the 1D advection-diffusion-reaction equation
>> dc/dt+u*dc/dx=D*dc2/dx2-kC using Fortan code but I`m still
>> facing some issues.
>
> See: https://class.coursera.org/scientificcomp-003

Sounds interesting but the link does not seem to work.
(The page shows a short *description* of the course but
there seems to be no way to actually see it, also not
after creating an account.)

> he does it in Matlab, and it also works in Octave.
>
> Actually, that is in 2D, which makes much nicer pictures.
> You can cheat and go directly to lecture 19, 20, or 21.

No. There is no list of lecture numbers, let alone
a way to "goto" any of them. :-)

-- 
Jos

Jos Bergervoet <jos.bergervoet@xs4all.nl> wrote:

(snip, I wrote)

>> See: https://class.coursera.org/scientificcomp-003
 
> Sounds interesting but the link does not seem to work.
> (The page shows a short *description* of the course but
> there seems to be no way to actually see it, also not
> after creating an account.)

I will have to see about that.

I can see it because I previously signed up. 

I am not sure about the schedule. Some will let you sign up
anytime, but others won't.

-- glen

Jos Bergervoet <jos.bergervoet@xs4all.nl> wrote:

(snip, I wrote)

>> See: https://class.coursera.org/scientificcomp-003

Seems to be: https://class.coursera.org/scientificcomp-006

I didn't realize that the links go to a specific session of the
course.  The -006 is the current session, that started Aug 26th.

Unlike some courses, for this one all homework and quizzes are due
on the last day (perfect for procrastinators) which seems to
be Oct 31, 11:59PM PDT.  (If you set the right time zone, the
assignments will do it right. In any case PDT is UTC-0700.

You will find that the instructor doesn't follow the discussion
forum all that closely, but if you already know about
Advection-Diffusion that might not be a problem.

You can decide from your background how closely, or not, to follow
the lecture, quiz, and homework schedule.  There is no cerificate,
so it is only for personal gratification and learning.

For one of the problems, you do a 2D advection-diffusion equation
and watch it move in color. 

Much of the lecture discussion is on stability of solutions
to differential equations.  Good stuff to know.

-- glen

On 9/20/2014 10:48 AM, glen herrmannsfeldt wrote:
> Jos Bergervoet<jos.bergervoet@xs4all.nl>  wrote:
>
> (snip, I wrote)
>
>>> See: https://class.coursera.org/scientificcomp-003
>
> Seems to be: https://class.coursera.org/scientificcomp-006

That page contains a text titled "Announcements", which
might be a correct title, but there is no actual entry
link to the course. From the links in the left-hand
pane a tried a few:

https://class.coursera.org/scientificcomp-006/wiki/WelcometotheCourse
gives again only some description of the course. But then:
https://class.coursera.org/scientificcomp-006/lecture
finally shows some entry points, only it seems to be
unrelated to the numbering you used (in your advice to
go to lecture 19, 20, or 21.) Looking, for instance, for
the text "Advection-diffusion", I find some sections that
might be appropriate, like this one:

https://class.coursera.org/scientificcomp-006/lecture/45

I'm not sure whether it is the lecture you actually meant,
but it clearly is a video lecture. (I'll leave it to those
readers here whit lots of time to further inspect the
contents!)

  ..
> Much of the lecture discussion is on stability of solutions
> to differential equations.  Good stuff to know.

It seems to contain a lot about the equations. Does he
also treat the solvers? (Stuff like all those different
conjugate-gradient implementations, pre-conditioners,
multi-grid approach, etc. I don't think I see it listed
on the pages I scanned through..)

-- 
Jos

Jos Bergervoet <jos.bergervoet@xs4all.nl> wrote:

>> Seems to be: https://class.coursera.org/scientificcomp-006
 
> That page contains a text titled "Announcements", which
> might be a correct title, but there is no actual entry
> link to the course. From the links in the left-hand
> pane a tried a few:
 
> https://class.coursera.org/scientificcomp-006/wiki/WelcometotheCourse
> gives again only some description of the course. But then:
> https://class.coursera.org/scientificcomp-006/lecture
> finally shows some entry points, only it seems to be
> unrelated to the numbering you used (in your advice to
> go to lecture 19, 20, or 21.) Looking, for instance, for
> the text "Advection-diffusion", I find some sections that
> might be appropriate, like this one:

It seems that the lectures get added week by week.

But the lecture notes are available:

http://spark-public.s3.amazonaws.com/scientificcomp/Lecture%20Note%20Packet/LectureNotePacketSciComputing.pdf

(Maybe only if you are signed up).

It says that the lecture notes are for classroom and personal research
use only. 
 
(snip)

> It seems to contain a lot about the equations. Does he
> also treat the solvers? (Stuff like all those different
> conjugate-gradient implementations, pre-conditioners,
> multi-grid approach, etc. I don't think I see it listed
> on the pages I scanned through..)

It is close to the beginning of solvers, but yes.

Well, for Advection-Diffusion you have to solve both the time
dependence and space dependence. 

For time dependence (initial value) he goes through Euler, Runga-Kutta,
and Adams, and mostly uses the RK that is part of Matlab and Octave.
It is pretty easy, you write the function and pass it to the appropriate
routine to call.

The spatial solution is done with Fourier transforms. 

Then to schemes like Leap-Frog (2,2) for combining the time step
and spatial solution for a stable result. 

As always, bigger time step for faster solution, smaller for more
accurate solution, but you have to properly consider the time step
and spatial mesh size.

But much of the discussion is in the direction that the examples
and problems are going to use.

Oh, the way the assignments work is that you write and run the
program, and generate the values of certain time/space points.
You then put those in, and, if close enough, you get the points.

If you use Octave instead of Matlab, sometimes the nubmers are
slightly different, but usually close enough.

I would solve it in Octave, then go through and translate 
to Fortran.  The debugging will be a lot faster.


-- glen

Although we didn't include a source term, which appears to be your primary =
difficulty, Chapter 12 of my book has a parallel Burgers equation solver th=
at could be a good starting point for you. See the Resources tab at http://=
www.cambridge.org/rouson.  If you email me, I would be glad to answer quest=
ions related to the Burgers equation solver.

The Cray, Intel, and GNU Fortran compilers are all capable of compiling it.=
  (With GNU, you'll need the GCC 5.0, the current development branch. If yo=
u want to do parallel runs, you'll also need OpenCoarrays from http://openc=
oarrays.org. But you can do serial runs without OpenCoarrays.)

Some aspects of a related heat equation solver are described in the modern =
Fortran videos at http://jolts.stanford.edu/math.=20

Damian

On 9/20/2014 12:12 PM, glen herrmannsfeldt wrote:
> Jos Bergervoet<jos.bergervoet@xs4all.nl>  wrote:
  ...
>> unrelated to the numbering you used (in your advice to
>> go to lecture 19, 20, or 21.) Looking, for instance, for
>> the text "Advection-diffusion", I find some sections that
>> might be appropriate, like this one:
>
> It seems that the lectures get added week by week.
>
> But the lecture notes are available:
>
> http://spark-public.s3.amazonaws.com/scientificcomp/Lecture%20Note%20Packet/LectureNotePacketSciComputing.pdf

Yes, I would say that's the one to start with!

> It says that the lecture notes are for classroom and personal research
> use only.

That's exactly what we need then, it seems. (It saves
you hundreds of hours of watching tedious videos! You
can clearly see here why the written word is the main
foundation of our civilization.)

>> It seems to contain a lot about the equations. Does he
>> also treat the solvers? (Stuff like all those different
>> conjugate-gradient implementations, pre-conditioners,
>> multi-grid approach, etc. I don't think I see it listed
>> on the pages I scanned through..)
>
> It is close to the beginning of solvers, but yes.

Section 2.3 at least mentions them. But for a real
treatment I think "Conjugate Gradient Method Without
the Agonizing Pain" still rules supremely:
http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf

> Well, for Advection-Diffusion you have to solve both the time
> dependence and space dependence.

Actually Navier-Stokes would be the most interesting
example of such a problem (of Millennium-Prize fame!
Nathan should give it a dedicated section, I'd say..)

> For time dependence (initial value) he goes through Euler, Runga-Kutta,
> and Adams, and mostly uses the RK that is part of Matlab and Octave.

Which is useful in giving examples of how those Octave
routines are called. Saves you hours of reading through
tedious help files.

> It is pretty easy, you write the function and pass it to the appropriate
> routine to call.

I'm not sure, however, if we are really coming closer
to answering OP's question. (But then again, I forgot
what it was..)

> ...
> I would solve it in Octave, then go through and translate
> to Fortran.  The debugging will be a lot faster.

Well, debugging is only needed if you make errors! So
with today's structured programming techniques this
hardly seems a factor of importance..

-- 
Jos

On Sunday, September 21, 2014 2:28:57 AM UTC+10, Jos Bergervoet wrote:
> On 9/20/2014 12:12 PM, glen herrmannsfeldt wrote:
>
> > I would solve it in Octave, then go through and translate
> > to Fortran.  The debugging will be a lot faster.
> 
> Well, debugging is only needed if you make errors! So
> with today's structured programming techniques this
> hardly seems a factor of importance..

Debugging is always needed, particularly in Fortran,
for errors can creep in in many ways.
At the very least, runs with test data are needed in order to check
that the correct results are obtained.

Then there are compiler errors, that may need extensive
work to track down and identify.

We used to have a customer who, whenever he discovered a programming
error in his program, always proclaimed:
"Can't understand it! Never made a mistake before in my life!"

Op zaterdag 20 september 2014 04:08:34 UTC+2 schreef kader:

> Could any one help to provide a simple code for solving this advection-diffusion-reaction equation ? I would appriciate your help.
> 

You have not described the difficulties you experience, but have you thought of reducing the timestep? If the decay coefficient is too large, it will lead to instabilities (i.e. the concentration will grow in size with each time).

Alternatively, you can try an implicit method:

(C(t+1) - C(t))/Dt + udC(t+1)/dx = Dd2C(t+1)/dx2 - kC(t+1)

evaluating all the terms at the new time level. Then the method is stable, although you have to solve a set of linear equations now, but that is  straightforward enough.

Regards,

Arjen

Arjen Markus <arjen.markus895@gmail.com> wrote:

(snip)

> You have not described the difficulties you experience, but 
> have you thought of reducing the timestep? If the decay 
> coefficient is too large, it will lead to instabilities 
> (i.e. the concentration will grow in size with each time).

The previously mentioned Coursera course describes the 
interaction between time step, space step, and instability.
 
> Alternatively, you can try an implicit method:
 
> (C(t+1) - C(t))/Dt + udC(t+1)/dx = Dd2C(t+1)/dx2 - kC(t+1)
 
> evaluating all the terms at the new time level. 
> Then the method is stable, although you have to solve a set 
> of linear equations now, but that is  straightforward enough.

The lectures use a Fourier transform based method for the
spatial part. Works well in 2D or even 3D.

-- glen

Op maandag 22 september 2014 23:15:52 UTC+2 schreef glen herrmannsfeldt:

> 
> The previously mentioned Coursera course describes the 
> interaction between time step, space step, and instability.
> 
> 
> > Alternatively, you can try an implicit method:>  
> 
> > (C(t+1) - C(t))/Dt + udC(t+1)/dx = Dd2C(t+1)/dx2 - kC(t+1) 
> 
> > evaluating all the terms at the new time level. 
> > Then the method is stable, although you have to solve a set 
> > of linear equations now, but that is  straightforward enough.
> 
> The lectures use a Fourier transform based method for the
> spatial part. Works well in 2D or even 3D.
> 

Interesting - I am not familiar with this type of methods. (In most cases the equations I work with have some non-linear reaction term.) I should have a look at this course :). Thanks.

Regards,

Arjen

Arjen Markus <arjen.markus895@gmail.com> wrote:

(snip, I wrote)
>> The lectures use a Fourier transform based method for the
>> spatial part. Works well in 2D or even 3D.
 
> Interesting - I am not familiar with this type of methods. 
> (In most cases the equations I work with have some 
> non-linear reaction term.) I should have a look at this 
> course :). Thanks.

http://www.thefouriertransform.com/applications/differentialequations.php

(Except much easier to use omega instead of 2*pi*f.)

Fourier transforms convert a differential equation into an
algebraic equation. Solve that, inverse transform, and you have
the solution.

The example above, but with omega (w):

  y''(t)-y(t)=-g(t)

Fourier transform both sides, and remember that the Fourier transform
of the derivative of a function is -iw times the transform 
of the function.

    F{y'}=iw F{y} = iwY

So:  -w**2*Y(w)-Y(w) = -G(w)

solve for Y(w):    Y(w)=-G(w)/(1+w**2). 


Especially convenient for boundary value problems where the function
goes to zero at the boundary, or for periodic boundary conditions.

-- glen

Op dinsdag 23 september 2014 10:34:50 UTC+2 schreef glen herrmannsfeldt:

> 
> Especially convenient for boundary value problems where the function
> goes to zero at the boundary, or for periodic boundary conditions.
> 
> -- glen

Thanks for this clarification. It has been a while since I last used Fourier transforms or Laplace transforms in this way. Also, I was thinking of Fourier series instead. 

Regards,

Arjen

Thank you guys for your kind reply and guidence, I really appriciate your help

Actually I`m using finite difference method to solve the equation and using the Crank Nicolson scheme for the distretisation.
I`m puting here the methodology of the code but not sure if is it 100 percent correct specially for the time splitting approach

first is to solve the Advection diffusion equation only and get the system equation thatcan be solved using TDMA method. 

The total time step for the whole ADR equation is DT=DT1+DT2
where DT1 is the time step to solve the AD equation only from time (n) to time (n+1/2)and DT2 is the time step to solve the Reaction equation from time (n+1/2) to time (n+1). the solution of time (n+1/2) is considered as an initial solution for the reaction equation.

!.....START TIME LOOP 
!
      TIME=0.
      DT=DT1+DT2
!
      DO 100 NT=1,NTMAX
      TIME=NT*DT1 

!.....SOLVE EQUATION SYSTEM for the time step DT1 to find the concentration at time 1/2

   Ap*C1(i,n+1/2)+AE*C1(i+1,n+1/2)+AW*C1(i-1,n+1/2)=Ap*C1(i,n)+AE*C1(i+1,n)+AW*C1(i-1,n) 
!
      CALL TDMA

!      TIME=TIME+DT2
!  
!       DO i=2,NM1
!       AA=(1./DT1)+(k/2.)
!       BB=(1./DT1)-(k/2.)
!       C(i,n+1)=BB/AA*C1(i,n+1/2)
       END DO

100 continue

could you please advise if I`m wrong somewhere ? 

Regards,

kader

"kader"  wrote

> Actually I`m using finite difference method to solve the equation and 
> using the Crank Nicolson scheme for the distretisation.
> I`m puting here the methodology of the code but not sure if is it 100 
> percent correct specially for the time splitting approach

> first is to solve the Advection diffusion equation only and get the system 
> equation thatcan be solved using TDMA method.

> !.....SOLVE EQUATION SYSTEM for the time step DT1 to find the 
> concentration at time 1/2

> 
> Ap*C1(i,n+1/2)+AE*C1(i+1,n+1/2)+AW*C1(i-1,n+1/2)=Ap*C1(i,n)+AE*C1(i+1,n)+AW*C1(i-1,n)

>     CALL TDMA

By TDMA do you mean "TriDiagonal Matrix Algorithm"? I recently looked at a 
different problem in another Fortran related newsgroup that was improved by
replacing Gaussian Elimination to solve a linear system by the TDMA method.

TDMA takes arguments

sub diagonal
main diagonal
super diagonal
right hand side
solution vector
[size of vectors]

diagonals are indexed from row 1, even if it is missing.

-- e