PDE_1D_MG

Invokes a module, with the statement USE PDE_1D_MG, near the second line of the program unit. The integrator is provided with single or double precision arithmetic, and a generic named interface is provided. We do not recommend using 32-bit floating point arithmetic here. The routine is called within the following loop, and is entered with each value of IDO. The loop continues until a value of IDO results in an exit.

IDO=1

DO

CASE(IDO == 1) {Do required initialization steps}

CASE(IDO == 2) {Save solution, update T0 and TOUT }

IF{Finished with integration} IDO=3

CASE(IDO == 3) EXIT {Normal}

CASE(IDO == 4) EXIT {Due to errors}

CASE(IDO == 5) {Evaluate initial data}

CASE(IDO == 6) {Evaluate differential equations}

CASE(IDO == 7) {Evaluate boundary conditions}

CASE(IDO == 8) {Prepare to solve banded system}

CASE(IDO == 9) {Solve banded system}

CALL PDE_1D_MG (T0, TOUT, IDO, U, &

initial_conditions,&

pde_system_definition,&

boundary_conditions, IOPT)

END DO

The arguments to PDE_1D_MG are required or optional.

Required Arguments

T0—(Input/Output)
This is the value of the independent variable t where the integration of ut begins. It is set to the value TOUT on return.

TOUT—(Input)
This is the value of the independent variable t where the integration of ut ends. Note: Values of T0 < TOUT imply integration in the forward direction, while values of T0 > TOUT imply integration in the backward direction. Either direction is permitted.

IDO—(Input/Output)
This in an integer flag that directs program control and user action. Its value is used for initialization, termination, and for directing user response during reverse communication:

IDO = 1 This value is assigned by the user for the start of a new problem. Internally it causes allocated storage to be reallocated, conforming to the problem size. Various initialization steps are performed.

IDO = 2 This value is assigned by the routine when the integrator has successfully reached the end point, TOUT.

IDO = 3 This value is assigned by the user at the end of a problem. The routine is called by the user with this value. Internally it causes termination steps to be performed.

IDO = 4 This value is assigned by the integrator when a type FATAL or TERMINAL error condition has occurred, and error processing is set NOT to STOP for these types of errors. It is not necessary to make a final call to the integrator with IDO=3 in this case.

Values of IDO = 5,6,7,8,9 are reserved for applications that provide problem information or linear algebra computations using reverse communication. When problem information is provided using reverse communication, the differential equations, boundary conditions and initial data must all be given. The absence of optional subroutine names in the calling sequence directs the routine to use reverse communication. In the module PDE_1D_MG_INT, scalars and arrays for evaluating results are named below. The names are preceded by the prefix “s_pde_1d_mg_” or “d_pde_1d_mg_”, depending on the precision. We use the prefix “?_pde_1d_mg_”, for the appropriate choice.

IDO = 5 This value is assigned by the integrator, requesting data for the initial conditions. Following this evaluation the integrator is re-entered.

(Optional) Update the grid of values in array locations U(NPDE + 1, j) j = 2, …, N. This grid is returned to the user equally spaced, but can be updated as desired, provided the values are increasing.

(Required) Provide initial values for all components of the system at the grid of values U(NPDE + 1, j) j = 1, …, N. If the optional step of updating the initial grid is performed, then the initial values are evaluated at the updated grid.

IDO = 6 This value is assigned by the integrator, requesting data for the differential equations. Following this evaluation the integrator is re-entered. Evaluate the terms of the system of Equation 2. A default value of m = 0 is assumed, but this can be changed to one of the other choices, m = 1 or m = 2 . Use the optional argument IOPT(:) for that purpose. Put the values in the arrays as indicated.

The assign-to equality, a := b, used here and below, is read "the expression b is evaluated and then assigned to the location a."

If any of the functions cannot be evaluated, set pde_1d_mg_ires=3. Otherwise do not change its value.

IDO = 7 This value is assigned by the integrator, requesting data for the boundary conditions, as expressed in Equation 3. Following the evaluation the integrator is re-entered.

The value x∈{xL, xR}, and the logical flag pde_1d_mg_LEFT=.TRUE. for x = xL. It has the value pde_1d_mg_LEFT=.FALSE. for x = xR. If any of the functions cannot be evaluated, set pde_1d_mg_ires=3. Otherwise do not change its value.

IDO = 8 This value is assigned by the integrator, requesting the calling program to prepare for solving a banded linear system of algebraic equations. This value will occur only when the option for “reverse communication solving” is set in the array IOPT(:), with option PDE_1D_MG_REV_COMM_FACTOR_SOLVE. The matrix data for this system is in Band Storage Mode, described in the Reference Material for the IMSL Fortran Numerical Libraries.

IDO = 9 This value is assigned by the integrator , requesting the calling program to solve a linear system with the matrix defined as noted with IDO=8.

U(1:NPDE+1,1:N)—(Input/Output)
This assumed-shape array specifies Input information about the problem size and boundaries. The dimension of the problem is obtained from NPDE + 1 = size(U,1). The number of grid points is obtained by N = size(U,2). Limits for the variable x are assigned as input in array locations, U(NPDE + 1, 1) = xL, U(NPDE  + 1, N) = xR. It is not required to define U(NPDE + 1, j), j = 2, …, N-1. At completion, the array U(1:NPDE,1:N) contains the approximate solution value Ui(xj(TOUT),TOUT) in location U(I,J). The grid value xj(TOUT) is in location U(NPDE+1,J). Normally the grid values are equally spaced as the integration starts. Variable spaced grid values can be provided by defining them as Output from the subroutine initial_conditions or during reverse communication, IDO=5.

Optional Arguments

initial_conditions—(Input)
The name of an external subroutine, written by the user, when using forward communication. If this argument is not used, then reverse communication is used to provide the problem information. The routine gives the initial values for the system at the starting independent variable value T0. This routine can also provide a non-uniform grid at the initial value.

     SUBROUTINE initial_conditions (NPDE,N,U)
        Integer NPDE,N
        REAL(kind(T0)) U(:,)
     END SUBROUTINE

(Optional) Update the grid of values in array locations U(NPDE) + 1, j), j = 2, …, N‑1. This grid is input equally spaced, but can be updated as desired, provided the values are increasing.

(Required) Provide initial values U(:, j), j = 1, …, N for all components of the system at the grid of values U(NPDE) + 1, j), j = 1, …, N. If the optional step of updating the initial grid is performed, then the initial values are evaluated at the updated grid.

pde_system_definition—(Input)
The name of an external subroutine, written by the user, when using forward communication. It gives the differential equation, as expressed in Equation 2.

      SUBROUTINE pde_system_definition (t, x, NPDE, r, dudx, c, q, r, IRES)
          Integer NPDE, IRES
          REAL(kind(T0)) t, x, u(:), dudx(:)
          REAL(kind(T0)) c(:,:), q(:), r(:)
      END SUBROUTINE

Evaluate the terms of the system of equations. A default value of m = 0 is assumed, but this can be changed to one of the other choices m = 1 or 2. Use the optional argument IOPT(:) for that purpose. Put the values in the arrays as indicated.

If any of the functions cannot be evaluated, set IRES=3. Otherwise do not change its value.

boundary_conditions—(Input)
The name of an external subroutine, written by the user when using forward communication. It gives the boundary conditions, as expressed in Equation 2.

    SUBROUTINE BOUNDARY_CONDITIONS(T,BETA,GAMMA,U,DUDX,NPDE,LEFT,IRES)
            real(kind(1d0)),intent(in)::t
        real(kind(1d0)),intent(out),dimension(:)::BETA, GAMMA
        real(kind(1d0)),intent(in),dimension(:)::U,DUDX
        integer,intent(in)::NPDE
        logical,intent(in)::LEFT
        integer,intent(out)::IRES
    END SUBROUTINE

The value , and the logical flag LEFT=.TRUE. for x = xL. The flag has the value LEFT=.FALSE. for x = xR.

IOPT—(Input)
Derived type array s_options or d_options, used for passing optional data to PDE_1D_MG. See the section Optional Data in the Introduction for an explanation of the derived type and its use. It is necessary to invoke a module, with the statement USE ERROR_OPTION_PACKET, near the second line of the program unit. Examples 2-8 use this optional argument. The choices are as follows:

IOPT(IO) = PDE_1D_MG_CART_COORDINATES
Use the value m = 0 in Equation 2. This is the default.

IOPT(IO) = PDE_1D_MG_CYL_COORDINATES
Use the value m = 1 in Equation 2. The default value is m = 0.

IOPT(IO) = PDE_1D_MG_SPH_COORDINATES
Use the value m = 2 in Equation 2. The default value is m = 0.

IOPT(IO) = ?_OPTIONS(PDE_1D_MG_TIME_SMOOTHING,TAU)
This option resets the value of the parameter  ≥ 0 described above.
The default value is  = 0.

IOPT(IO) = ?_OPTIONS(PDE_1D_MG_SPATIAL_SMOOTHING,KAP)
This option resets the value of the parameter κ ≥ 0, described above.
The default value is κ = 2

IOPT(IO) = ?_OPTIONS(PDE_1D_MG_MONITOR_REGULARIZING,ALPH)
This option resets the value of the parameter α ≥ 0, described above.
The default value is α = 0.01.

IOPT(IO) = ?_OPTIONS(PDE_1D_MG_RELATIVE_TOLERANCE,RTOL)
This option resets the value of the relative accuracy parameter used in DASPG.
The default value is RTOL=1E-2 for single precision and RTOL=1D-4 for double precision.

IOPT(IO) = ?_OPTIONS(PDE_1D_MG_ABSOLUTE_TOLERANCE,ATOL)
This option resets the value of the absolute accuracy parameter used in DASPG. The default value is ATOL=1E-2 for single precision andATOL=1D-4 for double precision.

IOPT(IO) = PDE_1D_MG_MAX_BDF_ORDER
IOPT(IO+1) = MAXBDF
Reset the maximum order for the BDF formulas used in DASPG. The default value is MAXBDF=2. The new value can be any integer between 1 and 5. Some problems will benefit by making this change. We used the default value due to the fact that DASPG may cycle on its selection of order and step-size with orders higher than value 2.

IOPT(IO) = PDE_1D_MG_REV_COMM_FACTOR_SOLVE
The calling program unit will solve the banded linear systems required in the stiff differential-algebraic equation integrator. Values of IDO=8, 9 will occur only when this optional value is used.

IOPT(IO) = PDE_1D_MG_NO_NULLIFY_STACK
To maintain an efficient interface, the routine PDE_1D_MG collapses the subroutine call stack with CALL_E1PSH(“NULLIFY_STACK”). This implies that the overhead of maintaining the stack will be eliminated, which may be important with reverse communication. It does not eliminate error processing. However, precise information of which routines have errors will not be displayed. To see the full call chain, this option should be used. Following completion of the integration, stacking is turned back on with CALL_E1POP(“NULLIFY_STACK”).

FORTRAN 90 Interface

Generic: CALL PDE_1D_MG (T0, TOUT, IDO, [, …])

Specific: The specific interface names are S_PDE_1D_MG and D_PDE_1D_MG.

Description

The equation

ut = f(u, x, t), xL < x < xR, t > t0,

is approximated at N time-dependent grid values

Using the total differential

transforms the differential equation to

Using central divided differences for the factor ux leads to the system of ordinary differential equations in implicit form

The terms Ui, Fi respectively represent the approximate solution to the partial differential equation and the value of f(u,x,t) at the point (x,t) = (xi,(t),t). The truncation error is second-order in the space variable, x. The above ordinary differential equations are underdetermined, so additional equations are added for the variation of the time-dependent grid points. It is necessary to discuss these equations, since they contain parameters that can be adjusted by the user. Often it will be necessary to modify these parameters to solve a difficult problem. For this purpose the following quantities are defined:

Note: The three-tiered equal sign, used here and below, is read "a≡b, or a and b are exactly the same object or value."

The values ni are the so-called point concentration of the grid, and κ ≥ 0 denotes a spatial smoothing parameter. Now the grid points are defined implicitly so that

where  ≥ 1 is a time-smoothing parameter. Choosing very large results in a fixed grid. Increasing the value of from its default avoids the error condition where grid lines cross. The divisors are

The value κ determines the level of clustering or spatial smoothing of the grid points. Decreasing κ from its default decrease the amount of spatial smoothing. The parameters Mi approximate arc length and help determine the shape of the grid or xi-distribution. The parameter prevents the grid movement from adjusting immediately to new values of the Mi, thereby avoiding oscillations in the grid that cause large relative errors. This is important when applied to solutions with steep gradients.

The discrete form of the differential equation and the smoothing equations are combined to yield the implicit system of differential equations.

This is frequently a stiff differential-algebraic system. It is solved using the integrator DASPG and its subroutines, including D2SPG. These are documented in this chapter. Note that DASPG is restricted to use within PDE_1D_MG until the routine exits with the flag IDO = 3. If DASPG is needed during the evaluations of the differential equations or boundary conditions, use of a second processor and inter-process communication is required. The only options for DASPG set by PDE_1D_MG are the Maximum BDF Order, and the absolute and relative error values, ATOL and RTOL. Users may set other options using the Options Manager. This is described in routine DASPG and generally in Chapter 11 of this manual.

Remarks on the Examples

Due to its importance and the complexity of its interface, this subroutine is presented with several examples. Many of the program features are exercised. The problems complete without any change to the optional arguments, except where these changes are required to describe or to solve the problem.

In many applications the solution to a PDE is used as an auxiliary variable, perhaps as part of a larger design or simulation process. The truncation error of the approximate solution is commensurate with piece-wise linear interpolation on the grid of values, at each output point.

Example 1 - Electrodynamics Model

This example is from Blom and Zegeling (1994). The system is

We make the connection between the model problem statement and the example:

The boundary conditions are

Rationale: Example 1

This is a non-linear problem with sharply changing conditions near t = 0. The default settings of integration parameters allow the problem to be solved. The use of PDE_1D_MG with forward communication requires three subroutines provided by the user to describe the initial conditions, differential equations, and boundary conditions.

program PDE_EX1

! Electrodynamics Model:

USE PDE_1d_mg_int

IMPLICIT NONE

INTEGER, PARAMETER :: NPDE=2, N=51, NFRAMES=5

INTEGER I, IDO

! Define array space for the solution.

real(kind(1d0)) U(NPDE+1,N), T0, TOUT

real(kind(1d0)) :: ZERO=0D0, ONE=1D0, &

DELTA_T=10D0, TEND=4D0

EXTERNAL IC_01, PDE_01, BC_01

! Start loop to integrate and write solution values.

IDO=1

DO

SELECT CASE (IDO)

! Define values that determine limits.

CASE (1)

T0=ZERO

TOUT=1D-3

U(NPDE+1,1)=ZERO;U(NPDE+1,N)=ONE

OPEN(FILE='PDE_ex01.out',UNIT=7)

WRITE(7, "(3I5, 4F10.5)") NPDE, N, NFRAMES,&

U(NPDE+1,1), U(NPDE+1,N), T0, TEND

! Update to the next output point.

! Write solution and check for final point.

CASE (2)

WRITE(7,"(F10.5)")TOUT

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

T0=TOUT;TOUT=TOUT*DELTA_T

IF(T0 >= TEND) IDO=3

TOUT=MIN(TOUT, TEND)

! All completed. Solver is shut down.

CASE (3)

CLOSE(UNIT=7)

EXIT

END SELECT

! Forward communication is used for the problem data.

CALL PDE_1D_MG (T0, TOUT, IDO, U,&

initial_conditions= IC_01,&

PDE_system_definition= PDE_01,&

boundary_conditions= BC_01)

END DO

END

SUBROUTINE IC_01(NPDE, NPTS, U)

! This is the initial data for Example 1.

IMPLICIT NONE

INTEGER NPDE, NPTS

REAL(KIND(1D0)) U(NPDE+1,NPTS)

U(1,:)=1D0;U(2,:)=0D0

END SUBROUTINE

SUBROUTINE PDE_01(T, X, NPDE, U, DUDX, C, Q, R, IRES)

! This is the differential equation for Example 1.

IMPLICIT NONE

INTEGER NPDE, IRES

REAL(KIND(1D0)) T, X, U(NPDE), DUDX(NPDE),&

C(NPDE,NPDE), Q(NPDE), R(NPDE)

REAL(KIND(1D0)) :: EPS=0.143D0, P=0.1743D0,&

ETA=17.19D0, Z, TWO=2D0, THREE=3D0

C=0D0;C(1,1)=1D0;C(2,2)=1D0

R=P*DUDX;R(1)=R(1)*EPS

Z=ETA*(U(1)-U(2))/THREE

Q(1)=EXP(Z)-EXP(-TWO*Z)

Q(2)=-Q(1)

END SUBROUTINE

SUBROUTINE BC_01(T, BTA, GAMA, U, DUDX, NPDE, LEFT, IRES)

! These are the boundary conditions for Example 1.

IMPLICIT NONE

INTEGER NPDE, IRES

LOGICAL LEFT

REAL(KIND(1D0)) T, BTA(NPDE), GAMA(NPDE),&

U(NPDE), DUDX(NPDE)

IF(LEFT) THEN

BTA(1)=1D0;BTA(2)=0D0

GAMA(1)=0D0;GAMA(2)=U(2)

ELSE

BTA(1)=0D0;BTA(2)=1D0

GAMA(1)=U(1)-1D0;GAMA(2)=0D0

END IF

END SUBROUTINE

Example 2 - Inviscid Flow on a Plate

This example is a first order system from Pennington and Berzins, (1994). The equations are

Following elimination of w, there remain NPDE = 2 differential equations. The variable t is not time, but a second space variable. The integration goes from t = 0 to t = 5. It is necessary to truncate the variable x at a finite value, say xmax = xR = 25. In terms of the integrator, the system is defined by letting m = 0 and

The boundary conditions are satisfied by

We use N = 10 + 51 = 61 grid points and output the solution at steps of Δt = 0.1.

Rationale: Example 2

This is a non-linear boundary layer problem with sharply changing conditions near t = 0. The problem statement was modified so that boundary conditions are continuous near t = 0. Without this change the underlying integration software, DASPG, cannot solve the problem. The continuous blending function u = exp(-20t) is arbitrary and artfully chosen. This is a mathematical change to the problem, required because of the stated discontinuity at t = 0. Reverse communication is used for the problem data. No additional user-written subroutines are required when using reverse communication. We also have chosen 10 of the initial grid points to be concentrated near xL = 0, anticipating rapid change in the solution near that point. Optional changes are made to use a pure absolute error tolerance and non-zero time-smoothing.

program PDE_1D_MG_EX02

! Inviscid Flow Over a Plate

USE PDE_1d_mg_int

USE ERROR_OPTION_PACKET

IMPLICIT NONE

INTEGER, PARAMETER :: NPDE=2, N1=10, N2=51, N=N1+N2

INTEGER I, IDO, NFRAMES

! Define array space for the solution.

real(kind(1d0)) U(NPDE+1,N), T0, TOUT, DX1, DX2, DIF

real(kind(1d0)) :: ZERO=0D0, ONE=1D0, DELTA_T=1D-1,&

TEND=5D0, XMAX=25D0

real(kind(1d0)) :: U0=1D0, U1=0D0, TDELTA=1D-1, TOL=1D-2

TYPE(D_OPTIONS) IOPT(3)

! Start loop to integrate and record solution values.

IDO=1

DO

SELECT CASE (IDO)

! Define values that determine limits and options.

CASE (1)

T0=ZERO

TOUT=DELTA_T

U(NPDE+1,1)=ZERO;U(NPDE+1,N)=XMAX

OPEN(FILE='PDE_ex02.out',UNIT=7)

NFRAMES=NINT((TEND+DELTA_T)/DELTA_T)

WRITE(7, "(3I5, 4D14.5)") NPDE, N, NFRAMES,&

U(NPDE+1,1), U(NPDE+1,N), T0, TEND

DX1=XMAX/N2;DX2=DX1/N1

IOPT(1)=D_OPTIONS(PDE_1D_MG_RELATIVE_TOLERANCE,ZERO)

IOPT(2)=D_OPTIONS(PDE_1D_MG_ABSOLUTE_TOLERANCE,TOL)

IOPT(3)=D_OPTIONS(PDE_1D_MG_TIME_SMOOTHING,1D-3)

! Update to the next output point.

! Write solution and check for final point.

CASE (2)

T0=TOUT

IF(T0 <= TEND) THEN

WRITE(7,"(F10.5)")TOUT

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

TOUT=MIN(TOUT+DELTA_T,TEND)

IF(T0 == TEND)IDO=3

END IF

! All completed. Solver is shut down.

CASE (3)

CLOSE(UNIT=7)

EXIT

! Define initial data values.

CASE (5)

U(:NPDE,:)=ZERO;U(1,:)=ONE

DO I=1,N1

U(NPDE+1,I)=(I-1)*DX2

END DO

DO I=N1+1,N

U(NPDE+1,I)=(I-N1)*DX1

END DO

WRITE(7,"(F10.5)")T0

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

! Define differential equations.

CASE (6)

D_PDE_1D_MG_C=ZERO

D_PDE_1D_MG_C(1,1)=ONE

D_PDE_1D_MG_C(2,1)=D_PDE_1D_MG_U(1)

D_PDE_1D_MG_R(1)=-D_PDE_1D_MG_U(2)

D_PDE_1D_MG_R(2)= D_PDE_1D_MG_DUDX(1)

D_PDE_1D_MG_Q(1)= ZERO

D_PDE_1D_MG_Q(2)= &

D_PDE_1D_MG_U(2)*D_PDE_1D_MG_DUDX(1)

! Define boundary conditions.

CASE (7)

D_PDE_1D_MG_BETA=ZERO

IF(PDE_1D_MG_LEFT) THEN

DIF=EXP(-20D0*D_PDE_1D_MG_T)

! Blend the left boundary value down to zero.

D_PDE_1D_MG_GAMMA=(/D_PDE_1D_MG_U(1)-DIF,D_PDE_1D_MG_U(2)/)

ELSE

D_PDE_1D_MG_GAMMA=(/D_PDE_1D_MG_U(1)-ONE,D_PDE_1D_MG_DUDX(2)/)

END IF

END SELECT

! Reverse communication is used for the problem data.

CALL PDE_1D_MG (T0, TOUT, IDO, U, IOPT=IOPT)

END DO

end program

Example 3 - Population Dynamics

This example is from Pennington and Berzins (1994). The system is

This is a notable problem because it involves the unknown

across the entire domain. The software can solve the problem by introducing two dependent algebraic equations:

This leads to the modified system

In the interface to the evaluation of the differential equation and boundary conditions, it is necessary to evaluate the integrals, which are computed with the values of u(x,t) on the grid. The integrals are approximated using the trapezoid rule, commensurate with the truncation error in the integrator.

Rationale: Example 3

This is a non-linear integro-differential problem involving non-local conditions for the differential equation and boundary conditions. Access to evaluation of these conditions is provided using reverse communication. It is not possible to solve this problem with forward communication, given the current subroutine interface. Optional changes are made to use an absolute error tolerance and non-zero time-smoothing. The time-smoothing value  = 1 prevents grid lines from crossing.

program PDE_1D_MG_EX03

! Population Dynamics Model.

USE PDE_1d_mg_int

USE ERROR_OPTION_PACKET

IMPLICIT NONE

INTEGER, PARAMETER :: NPDE=1, N=101

INTEGER IDO, I, NFRAMES

! Define array space for the solution.

real(kind(1d0)) U(NPDE+1,N), MID(N-1), T0, TOUT, V_1, V_2

real(kind(1d0)) :: ZERO=0D0, HALF=5D-1, ONE=1D0,&

TWO=2D0, FOUR=4D0, DELTA_T=1D-1,TEND=5D0, A=5D0

TYPE(D_OPTIONS) IOPT(3)

! Start loop to integrate and record solution values.

IDO=1

DO

SELECT CASE (IDO)

! Define values that determine limits.

CASE (1)

T0=ZERO

TOUT=DELTA_T

U(NPDE+1,1)=ZERO;U(NPDE+1,N)=A

OPEN(FILE='PDE_ex03.out',UNIT=7)

NFRAMES=NINT((TEND+DELTA_T)/DELTA_T)

WRITE(7, "(3I5, 4D14.5)") NPDE, N, NFRAMES,&

U(NPDE+1,1), U(NPDE+1,N), T0, TEND

IOPT(1)=D_OPTIONS(PDE_1D_MG_RELATIVE_TOLERANCE,ZERO)

IOPT(2)=D_OPTIONS(PDE_1D_MG_ABSOLUTE_TOLERANCE,1D-2)

IOPT(3)=D_OPTIONS(PDE_1D_MG_TIME_SMOOTHING,1D0)

! Update to the next output point.

! Write solution and check for final point.

CASE (2)

T0=TOUT

IF(T0 <= TEND) THEN

WRITE(7,"(F10.5)")TOUT

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

TOUT=MIN(TOUT+DELTA_T,TEND)

IF(T0 == TEND)IDO=3

END IF

! All completed. Solver is shut down.

CASE (3)

CLOSE(UNIT=7)

EXIT

! Define initial data values.

CASE (5)

U(1,:)=EXP(-U(2,:))/(TWO-EXP(-A))

WRITE(7,"(F10.5)")T0

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

! Define differential equations.

CASE (6)

D_PDE_1D_MG_C(1,1)=ONE

D_PDE_1D_MG_R(1)=-D_PDE_1D_MG_U(1)

! Evaluate the approximate integral, for this t.

V_1=HALF*SUM((U(1,1:N-1)+U(1,2:N))*&

(U(2,2:N) - U(2,1:N-1)))

D_PDE_1D_MG_Q(1)=V_1*D_PDE_1D_MG_U(1)

! Define boundary conditions.

CASE (7)

IF(PDE_1D_MG_LEFT) THEN

! Evaluate the approximate integral, for this t.

! A second integral is needed at the edge.

V_1=HALF*SUM((U(1,1:N-1)+U(1,2:N))*&

(U(2,2:N) - U(2,1:N-1)))

MID=HALF*(U(2,2:N)+U(2,1:N-1))

V_2=HALF*SUM(MID*EXP(-MID)*&

(U(1,1:N-1)+U(1,2:N))*(U(2,2:N)-U(2,1:N-1)))

D_PDE_1D_MG_BETA=ZERO

D_PDE_1D_MG_GAMMA=G(ONE,D_PDE_1D_MG_T)*V_1*V_2/(V_1+ONE)**2-&

D_PDE_1D_MG_U

ELSE

D_PDE_1D_MG_BETA=ZERO

D_PDE_1D_MG_GAMMA=D_PDE_1D_MG_DUDX(1)

END IF

END SELECT

! Reverse communication is used for the problem data.

CALL PDE_1D_MG (T0, TOUT, IDO, U, IOPT=IOPT)

END DO

CONTAINS

FUNCTION G(z,t)

IMPLICIT NONE

REAL(KIND(1d0)) Z, T, G

G=FOUR*Z*(TWO-TWO*EXP(-A)+EXP(-T))**2

G=G/((ONE-EXP(-A))*(ONE-(ONE+TWO*A)*&

EXP(-TWO*A))*(1-EXP(-A)+EXP(-T)))

END FUNCTION

end program

Example 4 - A Model in Cylindrical Coordinates

This example is from Blom and Zegeling (1994). The system models a reactor-diffusion problem:

The axial direction z is treated as a time coordinate. The radius r is treated as the single space variable.

Rationale: Example 4

This is a non-linear problem in cylindrical coordinates. Our example illustrates assigning m = 1 in Equation 2. We provide an optional argument that resets this value from its default, m = 0. Reverse communication is used to interface with the problem data.

program PDE_1D_MG_EX04

! Reactor-Diffusion problem in cylindrical coordinates.

USE pde_1d_mg_int

USE error_option_packet

IMPLICIT NONE

INTEGER, PARAMETER :: NPDE=1, N=41

INTEGER IDO, I, NFRAMES

! Define array space for the solution.

real(kind(1d0)) T(NPDE+1,N), Z0, ZOUT

real(kind(1d0)) :: ZERO=0D0, ONE=1D0, DELTA_Z=1D-1,&

ZEND=1D0, ZMAX=1D0, BTA=1D-4, GAMA=1D0, EPS=1D-1

TYPE(D_OPTIONS) IOPT(1)

! Start loop to integrate and record solution values.

IDO=1

DO

SELECT CASE (IDO)

! Define values that determine limits.

CASE (1)

Z0=ZERO

ZOUT=DELTA_Z

T(NPDE+1,1)=ZERO;T(NPDE+1,N)=ZMAX

OPEN(FILE='PDE_ex04.out',UNIT=7)

NFRAMES=NINT((ZEND+DELTA_Z)/DELTA_Z)

WRITE(7, "(3I5, 4D14.5)") NPDE, N, NFRAMES,&

T(NPDE+1,1), T(NPDE+1,N), Z0, ZEND

IOPT(1)=PDE_1D_MG_CYL_COORDINATES

! Update to the next output point.

! Write solution and check for final point.

CASE (2)

IF(Z0 <= ZEND) THEN

WRITE(7,"(F10.5)")ZOUT

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")T(I,:)

END DO

ZOUT=MIN(ZOUT+DELTA_Z,ZEND)

IF(Z0 == ZEND)IDO=3

END IF

! All completed. Solver is shut down.

CASE (3)

CLOSE(UNIT=7)

EXIT

! Define initial data values.

CASE (5)

T(1,:)=ZERO

WRITE(7,"(F10.5)")Z0

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")T(I,:)

END DO

! Define differential equations.

CASE (6)

D_PDE_1D_MG_C(1,1)=ONE

D_PDE_1D_MG_R(1)=BTA*D_PDE_1D_MG_DUDX(1)

D_PDE_1D_MG_Q(1)= -GAMA*EXP(D_PDE_1D_MG_U(1)/&

(ONE+EPS*D_PDE_1D_MG_U(1)))

! Define boundary conditions.

CASE (7)

IF(PDE_1D_MG_LEFT) THEN

D_PDE_1D_MG_BETA=ONE; D_PDE_1D_MG_GAMMA=ZERO

ELSE

D_PDE_1D_MG_BETA=ZERO; D_PDE_1D_MG_GAMMA=D_PDE_1D_MG_U(1)

END IF

END SELECT

! Reverse communication is used for the problem data.

! The optional derived type changes the internal model

! to use cylindrical coordinates.

CALL PDE_1D_MG (Z0, ZOUT, IDO, T, IOPT=IOPT)

END DO

end program

Example 5 - A Flame Propagation Model

This example is presented more fully in Verwer, et al., (1989). The system is a normalized problem relating mass density u(x,t) and temperature v(x,t):

Rationale: Example 5

This is a non-linear problem. The example shows the model steps for replacing the banded solver in the software with one of the user’s choice. Reverse communication is used for the interface to the problem data and the linear solver. Following the computation of the matrix factorization in DL2CRB, we declare the system to be singular when the reciprocal of the condition number is smaller than the working precision. This choice is not suitable for all problems. Attention must be given to detecting a singularity when this option is used.

program PDE_1D_MG_EX05

! Flame propagation model

USE pde_1d_mg_int

USE ERROR_OPTION_PACKET

USE Numerical_Libraries, ONLY :&

dl2crb, dlfsrb

IMPLICIT NONE

INTEGER, PARAMETER :: NPDE=2, N=40, NEQ=(NPDE+1)*N

INTEGER I, IDO, NFRAMES, IPVT(NEQ)

! Define array space for the solution.

real(kind(1d0)) U(NPDE+1,N), T0, TOUT

! Define work space for the banded solver.

real(kind(1d0)) WORK(NEQ), RCOND

real(kind(1d0)) :: ZERO=0D0, ONE=1D0, DELTA_T=1D-4,&

TEND=6D-3, XMAX=1D0, BTA=4D0, GAMA=3.52D6

TYPE(D_OPTIONS) IOPT(1)

! Start loop to integrate and record solution values.

IDO=1

DO

SELECT CASE (IDO)

! Define values that determine limits.

CASE (1)

T0=ZERO

TOUT=DELTA_T

U(NPDE+1,1)=ZERO; U(NPDE+1,N)=XMAX

OPEN(FILE='PDE_ex05.out',UNIT=7)

NFRAMES=NINT((TEND+DELTA_T)/DELTA_T)

WRITE(7, "(3I5, 4D14.5)") NPDE, N, NFRAMES,&

U(NPDE+1,1), U(NPDE+1,N), T0, TEND

IOPT(1)=PDE_1D_MG_REV_COMM_FACTOR_SOLVE

! Update to the next output point.

! Write solution and check for final point.

CASE (2)

T0=TOUT

IF(T0 <= TEND) THEN

WRITE(7,"(F10.5)")TOUT

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

TOUT=MIN(TOUT+DELTA_T,TEND)

IF(T0 == TEND)IDO=3

END IF

! All completed. Solver is shut down.

CASE (3)

CLOSE(UNIT=7)

EXIT

! Define initial data values.

CASE (5)

U(1,:)=ONE; U(2,:)=2D-1

WRITE(7,"(F10.5)")T0

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

! Define differential equations.

CASE (6)

D_PDE_1D_MG_C=ZERO

D_PDE_1D_MG_C(1,1)=ONE; D_PDE_1D_MG_C(2,2)=ONE

D_PDE_1D_MG_R=D_PDE_1D_MG_DUDX

D_PDE_1D_MG_Q(1)= D_PDE_1D_MG_U(1)*F(D_PDE_1D_MG_U(2))

D_PDE_1D_MG_Q(2)= -D_PDE_1D_MG_Q(1)

! Define boundary conditions.

CASE (7)

IF(PDE_1D_MG_LEFT) THEN

D_PDE_1D_MG_BETA=ZERO;D_PDE_1D_MG_GAMMA=D_PDE_1D_MG_DUDX

ELSE

D_PDE_1D_MG_BETA(1)=ONE

D_PDE_1D_MG_GAMMA(1)=ZERO

D_PDE_1D_MG_BETA(2)=ZERO

IF(D_PDE_1D_MG_T >= 2D-4) THEN

D_PDE_1D_MG_GAMMA(2)=12D-1

ELSE

D_PDE_1D_MG_GAMMA(2)=2D-1+5D3*D_PDE_1D_MG_T

END IF

D_PDE_1D_MG_GAMMA(2)=D_PDE_1D_MG_GAMMA(2)-&

D_PDE_1D_MG_U(2)

END IF

CASE(8)

! Factor the banded matrix. This is the same solver used

! internally but that is not required. A user can substitute

! one of their own.

call dl2crb (neq, d_pde_1d_mg_a, pde_1d_mg_lda, &

pde_1d_mg_iband, pde_1d_mg_iband, d_pde_1d_mg_a, &

pde_1d_mg_lda, ipvt, rcond, work)

IF(rcond <= EPSILON(ONE)) pde_1d_mg_panic_flag = 1

CASE(9)

! Solve using the factored banded matrix.

call dlfsrb(neq, d_pde_1d_mg_a, pde_1d_mg_lda, &

pde_1d_mg_iband, pde_1d_mg_iband, ipvt, &

d_pde_1d_mg_rhs, 1, d_pde_1d_mg_sol)

END SELECT

! Reverse communication is used for the problem data.

CALL PDE_1D_MG (T0, TOUT, IDO, U, IOPT=IOPT)

END DO

CONTAINS

FUNCTION F(Z)

IMPLICIT NONE

REAL(KIND(1D0)) Z, F

F=GAMA*EXP(-BTA/Z)

END FUNCTION

end program

Example 6 - A ‘Hot Spot’ Model

This example is presented more fully in Verwer, et al., (1989). The system is a normalized problem relating the temperature u(x,t), of a reactant in a chemical system. The formula for h(z) is equivalent to their example.

Rationale: Example 6

This is a non-linear problem. The output shows a case where a rapidly changing front, or hot-spot, develops after a considerable way into the integration. This causes rapid change to the grid. An option sets the maximum order BDF formula from its default value of 2 to the theoretical stable maximum value of 5.

USE pde_1d_mg_int

USE error_option_packet

IMPLICIT NONE

INTEGER, PARAMETER :: NPDE=1, N=80

INTEGER I, IDO, NFRAMES

! Define array space for the solution.

real(kind(1d0)) U(NPDE+1,N), T0, TOUT

real(kind(1d0)) :: ZERO=0D0, ONE=1D0, DELTA_T=1D-2,&

TEND=29D-2, XMAX=1D0, A=1D0, DELTA=2D1, R=5D0

TYPE(D_OPTIONS) IOPT(2)

! Start loop to integrate and record solution values.

IDO=1

DO

SELECT CASE (IDO)

! Define values that determine limits.

CASE (1)

T0=ZERO

TOUT=DELTA_T

U(NPDE+1,1)=ZERO; U(NPDE+1,N)=XMAX

OPEN(FILE='PDE_ex06.out',UNIT=7)

NFRAMES=(TEND+DELTA_T)/DELTA_T

WRITE(7, "(3I5, 4D14.5)") NPDE, N, NFRAMES,&

U(NPDE+1,1), U(NPDE+1,N), T0, TEND

! Illustrate allowing the BDF order to increase

! to its maximum allowed value.

IOPT(1)=PDE_1D_MG_MAX_BDF_ORDER

IOPT(2)=5

! Update to the next output point.

! Write solution and check for final point.

CASE (2)

T0=TOUT

IF(T0 <= TEND) THEN

WRITE(7,"(F10.5)")TOUT

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

TOUT=MIN(TOUT+DELTA_T,TEND)

IF(T0 == TEND)IDO=3

END IF

! All completed. Solver is shut down.

CASE (3)

CLOSE(UNIT=7)

EXIT

! Define initial data values.

CASE (5)

U(1,:)=ONE

WRITE(7,"(F10.5)")T0

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

! Define differential equations.

CASE (6)

D_PDE_1D_MG_C=ONE

D_PDE_1D_MG_R=D_PDE_1D_MG_DUDX

D_PDE_1D_MG_Q= - H(D_PDE_1D_MG_U(1))

! Define boundary conditions.

CASE (7)

IF(PDE_1D_MG_LEFT) THEN

D_PDE_1D_MG_BETA=ZERO

D_PDE_1D_MG_GAMMA=D_PDE_1D_MG_DUDX

ELSE

D_PDE_1D_MG_BETA=ZERO

D_PDE_1D_MG_GAMMA=D_PDE_1D_MG_U(1)-ONE

END IF

END SELECT

! Reverse communication is used for the problem data.

CALL PDE_1D_MG (T0, TOUT, IDO, U, IOPT=IOPT)

END DO

CONTAINS

FUNCTION H(Z)

real(kind(1d0)) Z, H

H=(R/(A*DELTA))*(ONE+A-Z)*EXP(-DELTA*(ONE/Z-ONE))

END FUNCTION

end program

Example 7 - Traveling Waves

This example is presented more fully in Verwer, et al., (1989). The system is a normalized problem relating the interaction of two waves, u(x,t) and v(x,t) moving in opposite directions. The waves meet and reduce in amplitude, due to the non-linear terms in the equation. Then they separate and travel onward, with reduced amplitude.

Rationale: Example 7

This is a non-linear system of first order equations.

program PDE_1D_MG_EX07

! Traveling Waves

USE pde_1d_mg_int

USE error_option_packet

IMPLICIT NONE

INTEGER, PARAMETER :: NPDE=2, N=50

INTEGER I, IDO, NFRAMES

! Define array space for the solution.

real(kind(1d0)) U(NPDE+1,N), TEMP(N), T0, TOUT

real(kind(1d0)) :: ZERO=0D0, HALF=5D-1, &

ONE=1D0, DELTA_T=5D-2,TEND=5D-1, PI

TYPE(D_OPTIONS) IOPT(5)

! Start loop to integrate and record solution values.

IDO=1

DO

SELECT CASE (IDO)

! Define values that determine limits.

CASE (1)

T0=ZERO

TOUT=DELTA_T

U(NPDE+1,1)=-HALF; U(NPDE+1,N)=HALF

OPEN(FILE='PDE_ex07.out',UNIT=7)

NFRAMES=(TEND+DELTA_T)/DELTA_T

WRITE(7, "(3I5, 4D14.5)") NPDE, N, NFRAMES,&

U(NPDE+1,1), U(NPDE+1,N), T0, TEND

IOPT(1)=D_OPTIONS(PDE_1D_MG_TIME_SMOOTHING,1D-3)

IOPT(2)=D_OPTIONS(PDE_1D_MG_RELATIVE_TOLERANCE,ZERO)

IOPT(3)=D_OPTIONS(PDE_1D_MG_ABSOLUTE_TOLERANCE,1D-3)

IOPT(4)=PDE_1D_MG_MAX_BDF_ORDER

IOPT(5)=3

! Update to the next output point.

! Write solution and check for final point.

CASE (2)

T0=TOUT

IF(T0 <= TEND) THEN

WRITE(7,"(F10.5)")TOUT

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

TOUT=MIN(TOUT+DELTA_T,TEND)

IF(T0 == TEND)IDO=3

END IF

! All completed. Solver is shut down.

CASE (3)

CLOSE(UNIT=7)

EXIT

! Define initial data values.

CASE (5)

TEMP=U(3,:)

U(1,:)=PULSE(TEMP); U(2,:)=U(1,:)

WHERE (TEMP < -3D-1 .or. TEMP > -1D-1) U(1,:)=ZERO

WHERE (TEMP < 1D-1 .or. TEMP > 3D-1) U(2,:)=ZERO

WRITE(7,"(F10.5)")T0

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

! Define differential equations.

CASE (6)

D_PDE_1D_MG_C=ZERO

D_PDE_1D_MG_C(1,1)=ONE; D_PDE_1D_MG_C(2,2)=ONE

D_PDE_1D_MG_R=D_PDE_1D_MG_U

D_PDE_1D_MG_R(1)=-D_PDE_1D_MG_R(1)

D_PDE_1D_MG_Q(1)= 100D0*D_PDE_1D_MG_U(1)*D_PDE_1D_MG_U(2)

D_PDE_1D_MG_Q(2)= D_PDE_1D_MG_Q(1)

! Define boundary conditions.

CASE (7)

D_PDE_1D_MG_BETA=ZERO;D_PDE_1D_MG_GAMMA=D_PDE_1D_MG_U

END SELECT

! Reverse communication is used for the problem data.

CALL PDE_1D_MG (T0, TOUT, IDO, U, IOPT=IOPT)

END DO

CONTAINS

FUNCTION PULSE(Z)

real(kind(1d0)) Z(:), PULSE(SIZE(Z))

PI=ACOS(-ONE)

PULSE=HALF*(ONE+COS(10D0*PI*Z))

END FUNCTION

end program

Example 8 - Black-Scholes

The value of a European “call option”, , with exercise price e and expiration date T, satisfies the “asset-or-nothing payoff ” . Prior to expiration is estimated by the Black-Scholes differential equation

The parameters in the model are the risk-free interest rate, r, and the stock volatility, σ. The boundary conditions are and . This development is described in Wilmott, et al. (1995), pages 41-57. There are explicit solutions for this equation based on the Normal Curve of Probability. The normal curve, and the solution itself, can be efficiently computed with the IMSL function ANORDF, IMSL (1994), page 186. With numerical integration the equation itself or the payoff can be readily changed to include other formulas, , and corresponding boundary conditions. We use

e = 100, r = 0.08, T-t = 0.25, σ2 = 0.04, sL = 0, and sR = 150

Rationale: Example 8

This is a linear problem but with initial conditions that are discontinuous. It is necessary to use a positive time-smoothing value to prevent grid lines from crossing. We have used an absolute tolerance of 10-3. In $US, this is one-tenth of a cent.

program PDE_1D_MG_EX08

! Black-Scholes call price

USE pde_1d_mg_int

USE error_option_packet

IMPLICIT NONE

INTEGER, PARAMETER :: NPDE=1, N=100

INTEGER I, IDO, NFRAMES

! Define array space for the solution.

real(kind(1d0)) U(NPDE+1,N), T0, TOUT, SIGSQ, XVAL

real(kind(1d0)) :: ZERO=0D0, HALF=5D-1, ONE=1D0, &

DELTA_T=25D-3, TEND=25D-2, XMAX=150, SIGMA=2D-1, &

R=8D-2, E=100D0

TYPE(D_OPTIONS) IOPT(5)

! Start loop to integrate and record solution values.

IDO=1

DO

SELECT CASE (IDO)

! Define values that determine limits.

CASE (1)

T0=ZERO

TOUT=DELTA_T

U(NPDE+1,1)=ZERO; U(NPDE+1,N)=XMAX

OPEN(FILE='PDE_ex08.out',UNIT=7)

NFRAMES=NINT((TEND+DELTA_T)/DELTA_T)

WRITE(7, "(3I5, 4D14.5)") NPDE, N, NFRAMES,&

U(NPDE+1,1), U(NPDE+1,N), T0, TEND

SIGSQ=SIGMA**2

! Illustrate allowing the BDF order to increase

! to its maximum allowed value.

IOPT(1)=PDE_1D_MG_MAX_BDF_ORDER

IOPT(2)=5

IOPT(3)=D_OPTIONS(PDE_1D_MG_TIME_SMOOTHING,5D-3)

IOPT(4)=D_OPTIONS(PDE_1D_MG_RELATIVE_TOLERANCE,ZERO)

IOPT(5)=D_OPTIONS(PDE_1D_MG_ABSOLUTE_TOLERANCE,1D-2)

! Update to the next output point.

! Write solution and check for final point.

CASE (2)

T0=TOUT

IF(T0 <= TEND) THEN

WRITE(7,"(F10.5)")TOUT

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

TOUT=MIN(TOUT+DELTA_T,TEND)

IF(T0 == TEND)IDO=3

END IF

! All completed. Solver is shut down.

CASE (3)

CLOSE(UNIT=7)

EXIT

! Define initial data values.

CASE (5)

U(1,:)=MAX(U(NPDE+1,:)-E,ZERO) ! Vanilla European Call

U(1,:)=U(NPDE+1,:) ! Asset-or-nothing Call

WHERE(U(1,:) <= E) U(1,:)=ZERO ! on these two lines

WRITE(7,"(F10.5)")T0

DO I=1,NPDE+1

WRITE(7,"(4E15.5)")U(I,:)

END DO

! Define differential equations.

CASE (6)

XVAL=D_PDE_1D_MG_X

D_PDE_1D_MG_C=ONE

D_PDE_1D_MG_R=D_PDE_1D_MG_DUDX*XVAL**2*SIGSQ*HALF

D_PDE_1D_MG_Q=-(R-SIGSQ)*XVAL*D_PDE_1D_MG_DUDX+R*D_PDE_1D_MG_U

! Define boundary conditions.

CASE (7)

IF(PDE_1D_MG_LEFT) THEN

D_PDE_1D_MG_BETA=ZERO

D_PDE_1D_MG_GAMMA=D_PDE_1D_MG_U

ELSE

D_PDE_1D_MG_BETA=ZERO

D_PDE_1D_MG_GAMMA=D_PDE_1D_MG_DUDX(1)-ONE

END IF

END SELECT

! Reverse communication is used for the problem data.

CALL PDE_1D_MG (T0, TOUT, IDO, U, IOPT=IOPT)

END DO

end program

Example 9 - Electrodynamics, Parameters Studied with MPI

Note: For a detailed description of MPI Requirements see Dense Matrix Parallelism Using MPI in Chapter 10 of this manual.

This example, described above in Example 1, is from Blom and Zegeling (1994). The system parameters ɛ, ρ, and η, are varied, using uniform random numbers. The intervals studied are 0.1 ≤ ɛ ≤ 0.2, 0.1 ≤ ρ ≤ 0.2, and 10 ≤ η ≤ 20. Using N = 21 grid values and other program options, the elapsed time, parameter values, and the value are sent to the root node. This information is written on a file. The final summary includes the minimum value of and the maximum and average time per integration, per node.

Rationale: Example 9

This is a non-linear simulation problem. Using at least two integrating processors and MPI allows more values of the parameters to be studied in a given time than with a single processor. This code is valuable as a study guide when an application needs to estimate timing and other output parameters. The simulation time is controlled at the root node. An integration is started, after receiving results, within the first SIM_TIME seconds. The elapsed time will be longer than SIM_TIME by the slowest processor’s time for its last integration.

program PDE_1D_MG_EX09

! Electrodynamics Model, parameter study.

USE PDE_1d_mg_int

USE MPI_SETUP_INT

USE RAND_INT

USE SHOW_INT

IMPLICIT NONE

INCLUDE "mpif.h"

INTEGER, PARAMETER :: NPDE=2, N=21

INTEGER I, IDO, IERROR, CONTINUE, STATUS(MPI_STATUS_SIZE)

INTEGER, ALLOCATABLE :: COUNTS(:)

! Define array space for the solution.

real(kind(1d0)) :: U(NPDE+1,N), T0, TOUT

real(kind(1d0)) :: ZERO=0D0, ONE=1D0,DELTA_T=10D0, TEND=4D0

! SIM_TIME is the number of seconds to run the simulation.

real(kind(1d0)) :: EPS, P, ETA, Z, TWO=2D0, THREE=3D0, SIM_TIME=60D0

real(kind(1d0)) :: TIMES, TIMEE, TIMEL, TIME, TIME_SIM, V_MIN, &

DATA(5)

real(kind(1d0)), ALLOCATABLE :: AV_TIME(:), MAX_TIME(:)

TYPE(D_OPTIONS) IOPT(4), SHOW_IOPT(2)

TYPE(S_OPTIONS) SHOW_INTOPT(2)

MP_NPROCS=MP_SETUP(1)

MPI_NODE_PRIORITY=(/(I-1,I=1,MP_NPROCS)/)

! If NP_NPROCS=1, the program stops. Change

! MPI_ROOT_WORKS=.TRUE. if MP_NPROCS=1.

MPI_ROOT_WORKS=.FALSE.

IF(.NOT. MPI_ROOT_WORKS .and. MP_NPROCS == 1) STOP

ALLOCATE(AV_TIME(MP_NPROCS), MAX_TIME(MP_NPROCS), COUNTS(MP_NPROCS))

! Get time start for simulation timing.

TIME=MPI_WTIME()

IF(MP_RANK == 0) OPEN(FILE='PDE_ex09.out',UNIT=7)

SIMULATE: DO

! Pick random parameter values.

EPS=1D-1*(ONE+rand(EPS))

P=1D-1*(ONE+rand(P))

ETA=10D0*(ONE+rand(ETA))

! Start loop to integrate and communicate solution times.

IDO=1

! Get time start for each new problem.

DO

IF(.NOT. MPI_ROOT_WORKS .and. MP_RANK == 0) EXIT

SELECT CASE (IDO)

! Define values that determine limits.

CASE (1)

T0=ZERO

TOUT=1D-3

U(NPDE+1,1)=ZERO;U(NPDE+1,N)=ONE

IOPT(1)=PDE_1D_MG_MAX_BDF_ORDER

IOPT(2)=5

IOPT(3)=D_OPTIONS(PDE_1D_MG_RELATIVE_TOLERANCE,1D-2)

IOPT(4)=D_OPTIONS(PDE_1D_MG_ABSOLUTE_TOLERANCE,1D-2)

TIMES=MPI_WTIME()

! Update to the next output point.

! Write solution and check for final point.

CASE (2)

T0=TOUT;TOUT=TOUT*DELTA_T

IF(T0 >= TEND) IDO=3

TOUT=MIN(TOUT, TEND)

! All completed. Solver is shut down.

CASE (3)

TIMEE=MPI_WTIME()

EXIT

! Define initial data values.

CASE (5)

U(1,:)=1D0;U(2,:)=0D0

! Define differential equations.

CASE (6)

D_PDE_1D_MG_C=0D0;D_PDE_1D_MG_C(1,1)=1D0;D_PDE_1D_MG_C(2,2)=1D0

D_PDE_1D_MG_R=P*D_PDE_1D_MG_DUDX

D_PDE_1D_MG_R(1)=D_PDE_1D_MG_R(1)*EPS

Z=ETA*(D_PDE_1D_MG_U(1)-D_PDE_1D_MG_U(2))/THREE

D_PDE_1D_MG_Q(1)=EXP(Z)-EXP(-TWO*Z)

D_PDE_1D_MG_Q(2)=-D_PDE_1D_MG_Q(1)

! Define boundary conditions.

CASE (7)

IF(PDE_1D_MG_LEFT) THEN

D_PDE_1D_MG_BETA(1)=1D0;D_PDE_1D_MG_BETA(2)=0D0

D_PDE_1D_MG_GAMMA(1)=0D0;D_PDE_1D_MG_GAMMA(2)=D_PDE_1D_MG_U(2)

ELSE

D_PDE_1D_MG_BETA(1)=0D0;D_PDE_1D_MG_BETA(2)=1D0

D_PDE_1D_MG_GAMMA(1)=D_PDE_1D_MG_U(1)- &

1D0;D_PDE_1D_MG_GAMMA(2)=0D0

END IF

END SELECT

! Reverse communication is used for the problem data.

CALL PDE_1D_MG (T0, TOUT, IDO, U)

END DO

TIMEL=TIMEE-TIMES

DATA=(/EPS, P, ETA, U(2,N), TIMEL/)

IF(MP_RANK > 0) THEN

! Send parameters and time to the root.

CALL MPI_SEND(DATA, 5, MPI_DOUBLE_PRECISION,0, MP_RANK, &

MP_LIBRARY_WORLD, IERROR)

! Receive back a "go/stop" flag.

CALL MPI_RECV(CONTINUE, 1, MPI_INTEGER, 0, MPI_ANY_TAG, &

MP_LIBRARY_WORLD, STATUS, IERROR)

! If root notes that time is up, it sends node a quit flag.

IF(CONTINUE == 0) EXIT SIMULATE

ELSE

! If root is working, record its result and then stand ready

! for other nodes to send.

IF(MPI_ROOT_WORKS) WRITE(7,*) MP_RANK, DATA

! If all nodes have reported, then quit.

IF(COUNT(MPI_NODE_PRIORITY >= 0) == 0) EXIT SIMULATE

! See if time is up. Some nodes still must report.

IF(MPI_WTIME()-TIME >= SIM_TIME) THEN

CONTINUE=0

ELSE

CONTINUE=1

END IF

! Root receives simulation data and finds which node sent it.

IF(MP_NPROCS > 1) THEN

CALL MPI_RECV(DATA, 5, MPI_DOUBLE_PRECISION, &

MPI_ANY_SOURCE, MPI_ANY_TAG, MP_LIBRARY_WORLD, &

STATUS, IERROR)

WRITE(7,*) STATUS(MPI_SOURCE), DATA

! If time at the root has elapsed, nodes receive signal to stop.

! Send the reporting node the "go/stop" flag.

! Mark if a node has been stopped.

CALL MPI_SEND(CONTINUE, 1, MPI_INTEGER, &

STATUS(MPI_SOURCE), &0, MP_LIBRARY_WORLD, IERROR)

IF (CONTINUE == 0) MPI_NODE_PRIORITY(STATUS(MPI_SOURCE)+1)&

=- MPI_NODE_PRIORITY(STATUS(MPI_SOURCE)+1)-1

END IF

IF (CONTINUE == 0) MPI_NODE_PRIORITY(1)=-1

END IF

END DO SIMULATE

IF(MP_RANK == 0) THEN

ENDFILE(UNIT=7);REWIND(UNIT=7)

! Read the data. Find extremes and averages.

MAX_TIME=ZERO;AV_TIME=ZERO;COUNTS=0;V_MIN=HUGE(ONE)

DO

READ(7,*, END=10) I, DATA

COUNTS(I+1)=COUNTS(I+1)+1

AV_TIME(I+1)=AV_TIME(I+1)+DATA(5)

IF(MAX_TIME(I+1) < DATA(5)) MAX_TIME(I+1)=DATA(5)

V_MIN=MIN(V_MIN, DATA(4))

END DO

10 CONTINUE

CLOSE(UNIT=7)

! Set printing Index to match node numbering.

SHOW_IOPT(1)= SHOW_STARTING_INDEX_IS

SHOW_IOPT(2)=0

SHOW_INTOPT(1)=SHOW_STARTING_INDEX_IS

SHOW_INTOPT(2)=0

CALL SHOW(MAX_TIME,"Maximum Integration Time, per process:",IOPT=SHOW_IOPT)

AV_TIME=AV_TIME/MAX(1,COUNTS)

CALL SHOW(AV_TIME,"Average Integration Time, per process:",IOPT=SHOW_IOPT)

CALL SHOW(COUNTS,"Number of Integrations",IOPT=SHOW_INTOPT)

WRITE(*,"(1x,A,F6.3)") "Minimum value for v(x,t),at x=1,t=4: ",V_MIN

END IF

MP_NPROCS=MP_SETUP("Final")

end program