Computer Program 1 - Demonstration Code - Fortran by rrk61112

VIEWS: 10 PAGES: 8

									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     NCARG_ROOT = "/u/ncsa/bjewett/ncarg"                                 read*,nstep
c     PATH += $NCARG_ROOT/bin
c                                                                          print*,'Enter plot interval, in steps:'
       program pgm1                                                        read*,nplot
       implicit none
c                                                                          print*,'Enter L for linear, N for nonlinear
c ... Definitions ...............................                                 advection:'
c                                                                          read*,reply
       integer nx,maxstep                                                  if (reply.eq.'L'.or.reply.eq.'l') then
       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*1 reply                                               gscr below.
       character*15 name                                         c
       logical linear                                                   call opngks
       parameter (name="Your name")                                     call gscr(1,0,1.,1.,1.)
                                                                        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
       enddo                                                     advection.'
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
               advection.\n");                                     *                           value of s2
        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;
               advection.\n");                                   {
      }                                                                  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