Calculation of the EMP product term 'p'


INTERFACE:

   subroutine findp_bisection(numc, cc, derivative, dt, accuracy, pi)
DESCRIPTION:

Auxiliary subroutine for finding the Extended Modified Patankar product term $ p$ with the bisection technique.

This subroutine solves the non-linear problem


$\displaystyle {\vec{c}^{n + 1} = \vec{c}^n + \Delta t \: \vec{f}(t^n,\vec{c}^n)\prod\limits_{j \in J^n} {\frac{c_j^{n + 1} }{c_j^n}}}$
    with $\displaystyle J^n = \left\{ {i:f_i (t^n,\vec{c}^n) < 0,i \in \{1,...,I\}} \right\}$ (261)

using the fact that it can be reduced to the problem of finding the root of a polynomial in one unknown $ p: = \prod\limits_{j \in J^n}{{c_j^{n + 1}}/{c_j^n}}$:

$\displaystyle g(p) = \prod\limits_{j \in J^n} {\left( {1 + \frac{\Delta t \: f_j(t^n,\vec{c}^n)}{c_j^n }p} \right)} - p = 0,$ (262)

with

$\displaystyle J^n = \left\{ {i:f_i (t^n,\vec{c}^n) < 0,i \in \{1,...,I\}} \right\},$ (263)

Additionally, it makes use of the the positivity requirement $ \vec{c}^{n+1}_i>0 \forall i$, which imposes restriction

$\displaystyle p \in \left( 0, \min \left( {1,\mathop {\min }\limits_{j \in J^n}...
...( { - \frac{c_j^n }{\Delta t \: f_j (t^n,\vec{c}^n)}} \right)} \right) \right).$ (264)

It has been proved that there exists exactly one $ p$ for which the above is true, see Bruggeman et al. (2006). The resulting problem is solved using the bisection scheme, which is guaranteed to converge.


USES:

   implicit none
INPUT PARAMETERS:
   integer, intent(in)   :: numc
   REALTYPE, intent(in)  :: cc(1:numc), derivative(1:numc)
   REALTYPE, intent(in)  :: dt, accuracy
  !OUTPUT PARAMETER:
   REALTYPE, intent(out) :: pi
REVISION HISTORY:
   Original author(s): Jorn Bruggeman
LOCAL VARIABLES:
   REALTYPE :: pileft, piright, fnow
   REALTYPE :: relderivative(1:numc)
   integer  :: iter, i, potnegcount

Karsten Bolding 2012-12-28