# Computer Program 1 - Demonstration Code - Fortran by rrk61112

VIEWS: 10 PAGES: 8

• pg 1
```									AUG. 25, 2009                                      ATMS 502 – CS 505 – CSE 566                                       Jewett
c
c ... Parameters and input .........................
Computer Program 1 – Demonstration Code – Fortran              c
c = 1.
c                                                                       dx = 0.1
c ATMS 502 / CS 505 / CSE 566 -- Fall, 2009
c Demo for pgm1: Linear and nonlinear advection                            print*,'Enter courant number:'
c Brian Jewett, University of Illinois                                     read*,courant
c                                                                          dt = courant*dx/c
c Compile with: "ncargf77 p1demo.f"
c You'll need the following in your home directory                         print*,'Number of steps for one complete loop =',
".soft" file (minus the Fortran comments & spaces here):               >      ( dx*real(nx) / c / dt )
c
c     # NCAR graphics                                                      print*,'Enter number of time steps:'
c     PATH += \$NCARG_ROOT/bin
c                                                                          print*,'Enter plot interval, in steps:'
implicit none
c                                                                          print*,'Enter L for linear, N for nonlinear
parameter (nx = 25, maxstep = 200)                                    linear = .true.
else
c ... Arrays and other variables                                             linear = .false.
endif
real s1(0:nx+1),s2(0:nx+1),strue(0:nx+1)                  c
real strace(maxstep),c,dt,dx,courant,smax                 c ... Open the NCAR Graphics package and set colors.
integer i,n,nstep,nplot                                   c ... To invert black/white colors, comment out calls to
character*15 name                                         c
logical linear                                                   call opngks
call gscr(1,1,0.,0.,0.)
c ... Variables for run history                                  c
c ... X-axis label
integer worksize                                          c
parameter (worksize = 2*nx*maxstep+nx+maxstep)                   call anotat('I','S',0,0,0,0)
real history(nx,maxstep)
real workarray(worksize)

Page 1
c                                                               c . . . . Plot fields when needed
c ... Set default Y plot bounds.                                         if (n .eq. nstep) then
c                                                                          call plot1d(s2(1:nx),nx,n,smax,1,strue(1:nx),
call agsetf('Y/MINIMUM.',-1.2)                                >              'Final solution',name)
call agsetf('Y/MAXIMUM.', 1.2)                                    else if (mod(n,nplot).eq.0) then
c                                                                          call plot1d(s2(1:nx),nx,n,smax,0,strue(1:nx),
c ... Set and plot the initial condition                             >              'Solution',name)
c                                                                        endif
call ic(s1,dx,nx)
call stats(s1,nx,0,smax)                                 c . . . . Check if problem out of bounds
call plot1d(s1(1:nx),nx,0,smax,0,strue(1:nx),                     if (smax.gt.1.5) then
>            'Initial condition',name)                                print*,'Stop - solution blowing up at step',n
c                                                                          call plot1d(s2(1:nx),nx,n,smax,0,strue(1:nx),
c ... Copy the initial condition to the "strue" array.               >               'Solution blowing up',name)
c                                                                          nstep = n
strue = s1                                                          go to 100
c                                                                        endif
c ... Integrate ......................................
c                                                               c . . . . end of time loop from n=1 to nstep
do n = 1,nstep                                                  enddo
c
c . . . . Set boundary conditions                               c ... Run complete - do final plots ......................
call bc(s1,nx)                                         c
100    continue
c . . . . Compute values at next step
call integrate(s1,s2,c,dx,dt,nx,linear)                c ... Plot Smax(t)
call agsetf('Y/MINIMUM.',0.)
c . . . . Do array update at end of time step.                         call agsetf('Y/MAXIMUM.',1.5)
call update(s1,s2,nx)                                         call anotat('N','Max value',0,0,0,0)
call plot1d(strace,nstep,0,smax,0,strue,
c . . . . Copy latest solution s2(1:nx) to history array.            >            'Smax vs. time',name)
history(1:nx,n) = s2(1:nx)
c ... Plot history (time-space) surface s(x,t)
c . . . . Stats                                                        call sfc(history,nx,nstep,dt*real(nstep),-90., 0.,
call stats(s2,nx,n,smax)                                    >        'Time-Space evolution (view: -90, 0)',name)
strace(n) = smax                                              call sfc(history,nx,nstep,dt*real(nstep),-30.,20.,
>        'Time-Space evolution (view: -30,+20)',name)

c
c ... Close graphics package and stop.
c
call clsgks

Page 2
stop                                                             real s1(0:nx+1)
end                                                        c
print*,'set boundary conditions here.'
c
return
c                                                                      end
c ============================ ic =====================          c
c IC sets the initial condition                                  c ======================= integrate ===================
c Arguments:                                                     c Integrate forward (advection only) by one time step.
c      s1      real array   holds initial condition              c Arguments:
c      dx      real         grid spacing                         c      s1      real array   values at current step
c      nx      integer      number of grid points                c      s2      real array   values at next step
c                                                                c      c       real         true speed of wave
subroutine ic(s1,dx,nx)                                   c      dx      real         grid spacing
implicit none                                             c      dt      real         time step
integer nx                                                c      nx      integer      number of grid points
real x,dx,s1(0:nx+1)                                      c      linear logical       if true, linear advection;
c                                                                c                           otherwise, nonlinear
integer i                                                 c
real pi,length                                                   subroutine integrate(s1,s2,c,dx,dt,nx,linear)
c                                                                       implicit none
pi = 4.0*atan(1.0)                                               integer i,nx
length = dx*real(nx)                                             real c,dt,dx,s1(0:nx+1),s2(0:nx+1),courant
c                                                                       logical linear
do i = 1,nx                                               c
x = dx*real(i-1)                                              if (linear) then
s1(i) = sin( 2.*pi/length*x )                                    print*,'Put integration code here for linear
c                                                                          do i = 1,nx
return                                                                s2(i) = 0.90*s1(i)
end                                                                 enddo
c                                                                       else
c ============================ bc =====================                    print*,'Put integration code here for nonlinear
c BC sets the boundary conditions                                advection.'
c Arguments:                                                            endif
c      s1      real array   values at current time step
c      nx      integer      main array size, not including             return
c                           extra 'ghost' zones/points                 end
c
subroutine bc(s1,nx)
implicit none
integer nx

Page 3
c                                                             c
c ========================= update ====================       c ========================= stats =====================
c Update: replace old values with new ones                    c Stats computes and prints out the max S values
c Arguments:                                                  c Arguments:
c      s1,s2 real arrays old, new data arrays                 c      s2      real array   latest array values
c      nx     integer             size of arrays              c      nx      integer             size of data array
c                                                             c      n       integer             time step counter
subroutine update(s1,s2,nx)                            c      smax    real         holds max absolute
implicit none                                          c                             value of s2
integer nx                                             c
real s1(0:nx+1),s2(0:nx+1)                                    subroutine stats(s2,nx,n,smax)
integer i                                                     implicit none
integer i,n,nx
s1(1:nx) = s2(1:nx)                                            real s2(0:nx+1),smax
c
return                                                         smax = abs(s2(1))
end                                                            do i = 2,nx
smax = max(smax,abs(s2(i)))
enddo
write(6,10) n,smax
10     format('Step ',i3,', Max = ',f7.5)

return
end

Page 4
/* Arrays and other variables */
Computer Program 1 – Demonstration Code – “C”
float s1[NXDIM],s2[NXDIM],strue[NXDIM];
/*                                                              float strace[MAXSTEP],c,dt,dx,courant,smax;
* ATMS 502 / CS 505 / CSE 566 -- Fall, 2009                    int i,j,n,nstep,nplot,linear;
* Demo for pgm1: Linear and nonlinear advection                char *name = "Your name";
* Brian Jewett, University of Illinois
*                                                        /* Variables for run history */
* Compile with: "ncargcc p1demo.c"
* You'll need the following in your home directory             float history[MAXSTEP][NX], workarray[WORKSIZE];
".soft" file (minus the C comments & spaces here):
*                                                        /* Variables to reverse default black/white colors in NCAR
*     # NCAR graphics                                        Graphics */
*     NCARG_ROOT = "/u/ncsa/bjewett/ncarg"
*     PATH += \$NCARG_ROOT/bin                                  Gcolr_rep rgb1,rgb2;
*/
/* Definitions for routines */
#include <stdio.h>
#include <stdlib.h>                                              void
#include <math.h>                                         ic(),stats(),plot1d(),sfc(),bc(),integrate(),update();

#include <ncarg/ncargC.h>                                 /* Parameters and input .............................. */
#include <ncarg/gks.h>
#define IWTYPE 1                                                printf("NX=%d, BC_WIDTH=%d, I1=%d, I2=%d,
#define WKID   1                                                       NXDIM=%d\n", NX,BC_WIDTH,I1,I2,NXDIM);

main()                                                          c = 1.0;
{                                                               dx = 0.1;

/*                                                              printf("Enter courant number: ");
* Definitions                                                  scanf("%f",&courant);
*/                                                             dt = courant*dx/c;

#define   NX 25                                                 printf("Number of steps for one complete loop =
#define   BC_WIDTH 1                                                  %.0f\n", ( dx*(float)NX / c / dt ));
#define   I1 BC_WIDTH
#define   I2 I1+NX-1                                            printf("Enter number of time steps: ");
#define   NXDIM NX+2*BC_WIDTH                                   scanf("%d",&nstep);
#define   MAXSTEP 200
#define   WORKSIZE (2*NX*MAXSTEP+NX+MAXSTEP)                    printf("Enter plot interval, in steps: ");
scanf("%d",&nplot);

Page 5
printf("Enter 1 for linear, 0 for nonlinear         /*
advection: ");                               * .. Integrate ....................................
scanf("%d",&linear);                                 */
/*                                                               for (n=1; n<=nstep; n++) {
* Open the NCAR Graphics package and set colors.
*/                                                       /*   . . . Set boundary conditions                     */
gopen_gks("stdout",0);                                       bc(s1,I1,I2,NX);
gopen_ws(WKID, NULL, IWTYPE);
gactivate_ws(WKID);                                /*   . . . Compute values at next step                 */
integrate(s1,s2,c,dx,dt,I1,I2,NX,linear);
/* omit following lines to invert black/white
colors */                                        /*   . . . Do array update at end of time step         */
rgb1.rgb.red = 1.;                                            update(s1,s2,I1,I2,NX);
rgb1.rgb.green = 1.;
rgb1.rgb.blue = 1.;                                 /*   . . . Copy latest solution s2() to history array */
rgb2.rgb.red = 0.;                                            for (i=I1,j=0; i<=I2; i++,j++)
rgb2.rgb.green = 0.;                                            history[n-1][j] = s2[i];
rgb2.rgb.blue = 0.;
gset_colr_rep(WKID,0,&rgb1);                        /*   . . . Stats                                       */
gset_colr_rep(WKID,1,&rgb2);                                  stats(s2,I1,I2,NX,n,&smax);
/*                                                                  strace[n-1] = smax;
* X-axis label
*/                                                       /*   . . . Plot fields when needed                     */
c_anotat("I","S",0,0,0,0);                                   if (n == nstep) {
/*                                                                    plot1d(s2,I1,I2,NX,n,smax,1,strue,"Final
* Set default Y plot bounds                                                 solution",name);
*/                                                                 } else if (n%nplot == 0) {
c_agsetf("Y/MINIMUM.",-1.2);                                   plot1d(s2,I1,I2,NX,n,smax,0,strue,
c_agsetf("Y/MAXIMUM.", 1.2);                                          "Solution",name);
}
/*
* Set and plot the initial condition                     /*   . . . Check if problem out of bounds              */
*/                                                                 if (smax > 1.5) {
ic(s1,dx,I1,I2,NX);                                            printf("Stop - solution blowing up at step
stats(s1,I1,I2,NX,0,&smax);                                            %d\n",n);
plot1d(s1,I1,I2,NX,0,smax,0,strue,"Initial                     plot1d(s2,I1,I2,NX,n,smax,0,strue,
condition",name);                                               "Solution blowing up",name);
/*                                                                    nstep = n;
* Copy the initial condition to the "strue" array.                   break;
*/                                                                 }
for (i=I1; i<=I2; i++) strue[i]=s1[i];
}      /* end of time loop n = 1,...,nstep */

Page 6
/*
* Run complete - do final plots                                void ic(s1,dx,i1,i2,nx)
*/                                                             float dx,s1[];
int i1,i2,nx;
/* . . Plot Smax(t)                                    */       {
c_agsetf("Y/MINIMUM.",0.0);                                     int i,j;
c_agsetf("Y/MAXIMUM.",1.5);                                     float x,pi,length;
c_anotat("N","Max value",0,0,0,0);
plot1d(strace,0,nstep-1,nstep,0,smax,0,strue,                  pi = 4.0*atan(1.0);
"Smax vs. time",name);                                  length = dx * (float) nx;

/* . . Plot history surface s(x,t)                    */              for (i=i1,j=1; i<=i2; i++,j++) {
sfc(history,nstep,NX,(dt*(float)nstep),                          x = dx * (float)(j-1);
-90., 0.,                                                      s1[i] = sin( 2.*pi/length*x );
"Time-Space evolution (view: -90, 0)",name);                }
sfc(history,nstep,NX,(dt*(float)nstep),
-30.,20.,                                                    return;
"Time-Space evolution (view: -30,+20)",name);          }
/*
* Close the graphics package and stop.                         /*
*/                                                              * ============================ bc =====================
* BC sets the boundary conditions
gdeactivate_ws(WKID);                                      * Arguments:
gclose_ws(WKID);                                           *     s1     real array   values at current time step
gclose_gks();                                              *     i1,i2 integers      indices bounding array data
*     nx     integer             main array size, not
exit;                                                     including
}                                                                *                         extra 'ghost' zones/points
*/
/*
* ============================ ic =====================        void bc(s1,i1,i2,nx)
* IC sets the initial condition                                float s1[];
* Arguments:                                                   int i1,i2,nx;
*     s1     real array   IC data. Set 1..nx here;             {
*                           [0],[nx+1] = ghost zones                  printf("Set boundary conditions here.\n");
*     dx     real         grid spacing
*     i1,i2 integers      indices bounding array data                return;
*     nx     integer             number of grid points         }
*
*/

Page 7
/*
* ======================= integrate ===================        void update(s1,s2,i1,i2,nx)
* Integrate forward (advection only) by one time step.         float s1[],s2[];
* Arguments:                                                   int i1,i2,nx;
*     s1     real array   values at current step               {
*     s2     real array   values at next step                         int i;
*     c      real         true speed of wave
*     dx     real         grid spacing                               for (i=i1; i<=i2; i++)
*     dt     real         time step                                    s1[i] = s2[i];
*     i1,i2 integers      indices bounding array data
*     nx     integer             number of grid points               return;
*     linear integer             if 1, linear advection;       }
*                         otherwise, nonlinear
*/                                                             /*
void integrate(s1,s2,c,dx,dt,i1,i2,nx,linear)                      * ========================= stats =====================
float s1[],s2[],c,dx,dt;                                           * Stats computes and prints out the max S values
int i1,i2,nx,linear;                                               * Arguments:
{                                                                  *     s2     real array   Latest data. Check 1..nx;
int i;                                                     *                           [0],[nx+1] = ghost zones
float courant;                                             *     i1,i2 integers      indices bounding array data
*     nx     integer             size of data array
if (linear) {                                                *     n      integer             time step counter
printf("Put integration code here for linear               *     smax   real         holds max absolute
for (i=i1; i<=i2; i++)                                     */
s2[i] = 0.90*s1[i];                                    void stats(s2,i1,i2,nx,n,smax)
} else {                                                   int i1,i2,nx,n;
printf("Put integration code here for nonlinear          float s2[],*smax;
}                                                                  int i;
float stmp;
return;
}                                                                      stmp = fabs(s2[i1]);
/*                                                                     for (i=i1+1; i<=i2; i++) {
* ========================= update ====================                if (fabs(s2[i]) > stmp) stmp = fabs(s2[i]);
* Update: replace old values with new ones                           }
* Arguments:
*     s1,s2 real arrays old, new data arrays                         printf("Step %3d, Max = %7.5f\n",n,stmp);
*     i1,i2 integers      indices bounding array data                *smax = stmp;
*     nx     integer             size of arrays                      return;
*                                                              }
*/

Page 8

```
To top