Distributed-Run-Time-Resource-Manager/RCCE_V2.0/apps/NPB/LU/erhs.c

// 
// Copyright 2010 Intel Corporation
// 
//    Licensed under the Apache License, Version 2.0 (the "License");
//    you may not use this file except in compliance with the License.
//    You may obtain a copy of the License at
// 
//        http://www.apache.org/licenses/LICENSE-2.0
// 
//    Unless required by applicable law or agreed to in writing, software
//    distributed under the License is distributed on an "AS IS" BASIS,
//    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//    See the License for the specific language governing permissions and
//    limitations under the License.
// 
#include "applu_share.h"
#include "applu_macros.h"

void erhs(){

//   compute the right hand side based on exact solution

//c---------------------------------------------------------------------
//c  local variables
//c---------------------------------------------------------------------
      int i, j, k, m;
      int iglob, jglob;
      int iex;
      int L1, L2;
      int ist1, iend1;
      int jst1, jend1;
      double  dsspm;
      double  xi, eta, zeta;
      double  q;
      double  u21, u31, u41;
      double  tmp;
      double  u21i, u31i, u41i, u51i;
      double  u21j, u31j, u41j, u51j;
      double  u21k, u31k, u41k, u51k;
      double  u21im1, u31im1, u41im1, u51im1;
      double  u21jm1, u31jm1, u41jm1, u51jm1;
      double  u21km1, u31km1, u41km1, u51km1;

      dsspm = dssp;

      for (k=1; k<=nz; k++)
         for (j=1; j<=ny; j++)
            for (i=1; i<=nx; i++)
               for (m=1; m<=5; m++)
                  frct( m, i, j, k ) = 0.0;

      for (k=1; k<=nz; k++) {
         zeta = ( (double)(k-1) ) / ( nz - 1 );
         for (j=1; j<=ny; j++) {
            jglob = jpt + j;
            eta = ( (double)(jglob-1) ) / ( ny0 - 1 );
            for (i=1; i<=nx; i++) {
               iglob = ipt + i;
               xi = ( (double)(iglob-1) ) / ( nx0 - 1 );
               for (m=1; m<=5; m++) {
                  rsd(m,i,j,k) =  ce(m,1)
                       + ce(m,2) * xi
                       + ce(m,3) * eta
                       + ce(m,4) * zeta
                       + ce(m,5) * xi * xi
                       + ce(m,6) * eta * eta
                       + ce(m,7) * zeta * zeta
                       + ce(m,8) * xi * xi * xi
                       + ce(m,9) * eta * eta * eta
                       + ce(m,10) * zeta * zeta * zeta
                       + ce(m,11) * xi * xi * xi * xi
                       + ce(m,12) * eta * eta * eta * eta
                       + ce(m,13) * zeta * zeta * zeta * zeta;
               }
            }
         }
      }

//c---------------------------------------------------------------------
//c   xi-direction flux differences
//c---------------------------------------------------------------------
//c
//c   iex = flag : iex = 0  north/south communication
//c              : iex = 1  east/west communication
//c
//c---------------------------------------------------------------------
      iex   = 0;

//c---------------------------------------------------------------------
//c   communicate and receive/send two rows of data
//c---------------------------------------------------------------------
//
      exchange_3(rsd,iex);

      L1 = 0;
      if (north == -1) L1 = 1;
      L2 = nx + 1;
      if (south == -1) L2 = nx;

      ist1 = 1;
      iend1 = nx;
      if (north == -1) ist1 = 4;
      if (south == -1) iend1 = nx - 3;

      for (k=2; k<=nz-1; k++)
         for (j=jst; j<=jend; j++)
            for (i=L1; i<=L2; i++) {
               flux(1,i,j,k) = rsd(2,i,j,k);
               u21 = rsd(2,i,j,k) / rsd(1,i,j,k);
               q = 0.50 * (  rsd(2,i,j,k) * rsd(2,i,j,k)
                               + rsd(3,i,j,k) * rsd(3,i,j,k)
                               + rsd(4,i,j,k) * rsd(4,i,j,k) )
                            / rsd(1,i,j,k);
               flux(2,i,j,k) = rsd(2,i,j,k) * u21 + c2 * 
                               ( rsd(5,i,j,k) - q );
               flux(3,i,j,k) = rsd(3,i,j,k) * u21;
               flux(4,i,j,k) = rsd(4,i,j,k) * u21;
               flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u21;
      }

      for (k=2; k<=nz-1; k++)
         for (j=jst; j<=jend; j++) {
            for (i=ist; i<=iend; i++)
               for (m=1; m<=5; m++)
                  frct(m,i,j,k) =  frct(m,i,j,k)
                         - tx2 * ( flux(m,i+1,j,k) - flux(m,i-1,j,k) );
            for (i=ist; i<=L2; i++) {
               tmp = 1.0 / rsd(1,i,j,k);

               u21i = tmp * rsd(2,i,j,k);
               u31i = tmp * rsd(3,i,j,k);
               u41i = tmp * rsd(4,i,j,k);
               u51i = tmp * rsd(5,i,j,k);

               tmp = 1.0 / rsd(1,i-1,j,k);

               u21im1 = tmp * rsd(2,i-1,j,k);
               u31im1 = tmp * rsd(3,i-1,j,k);
               u41im1 = tmp * rsd(4,i-1,j,k);
               u51im1 = tmp * rsd(5,i-1,j,k);

               flux(2,i,j,k) = (4.0/3.0) * tx3 * 
                              ( u21i - u21im1 );
               flux(3,i,j,k) = tx3 * ( u31i - u31im1 );
               flux(4,i,j,k) = tx3 * ( u41i - u41im1 );
               flux(5,i,j,k) = 0.50 * ( 1.0 - c1*c5 )
                    * tx3 * ( ( u21i*u21i + u31i*u31i + u41i*u41i )
                            - ( u21im1*u21im1 + u31im1*u31im1 + u41im1*u41im1 ) )
                    + (1.0/6.0)
                    * tx3 * ( u21i*u21i - u21im1*u21im1 )
                    + c1 * c5 * tx3 * ( u51i - u51im1 );
            }

            for (i=ist; i<=iend; i++) {
               frct(1,i,j,k) = frct(1,i,j,k)
                    + dx1 * tx1 * (            rsd(1,i-1,j,k)
                                   - 2.0 * rsd(1,i,j,k)
                                   +           rsd(1,i+1,j,k) );
               frct(2,i,j,k) = frct(2,i,j,k)
                 + tx3 * c3 * c4 * ( flux(2,i+1,j,k) - flux(2,i,j,k) )
                    + dx2 * tx1 * (            rsd(2,i-1,j,k)
                                   - 2.0 * rsd(2,i,j,k)
                                   +           rsd(2,i+1,j,k) );
               frct(3,i,j,k) = frct(3,i,j,k)
                 + tx3 * c3 * c4 * ( flux(3,i+1,j,k) - flux(3,i,j,k) )
                    + dx3 * tx1 * (            rsd(3,i-1,j,k)
                                   - 2.0 * rsd(3,i,j,k)
                                   +           rsd(3,i+1,j,k) );
               frct(4,i,j,k) = frct(4,i,j,k)
                  + tx3 * c3 * c4 * ( flux(4,i+1,j,k) - flux(4,i,j,k) )
                    + dx4 * tx1 * (            rsd(4,i-1,j,k)
                                   - 2.0 * rsd(4,i,j,k)
                                   +           rsd(4,i+1,j,k) );
               frct(5,i,j,k) = frct(5,i,j,k)
                 + tx3 * c3 * c4 * ( flux(5,i+1,j,k) - flux(5,i,j,k) )
                    + dx5 * tx1 * (            rsd(5,i-1,j,k)
                                   - 2.0 * rsd(5,i,j,k)
                                   +           rsd(5,i+1,j,k) );
            }

//c---------------------------------------------------------------------
//c   Fourth-order dissipation
//c---------------------------------------------------------------------
            if (north == -1) {
             for (m=1; m<=5; m++) {
               frct(m,2,j,k) = frct(m,2,j,k)
                 - dsspm * ( + 5.0 * rsd(m,2,j,k)
                             - 4.0 * rsd(m,3,j,k)
                             +           rsd(m,4,j,k) );
               frct(m,3,j,k) = frct(m,3,j,k)
                 - dsspm * ( - 4.0 * rsd(m,2,j,k)
                             + 6.0 * rsd(m,3,j,k)
                             - 4.0 * rsd(m,4,j,k)
                             +           rsd(m,5,j,k) );
             }
            }

            for (i=ist1; i<=iend1; i++)
               for (m=1; m<=5; m++)
                  frct(m,i,j,k) = frct(m,i,j,k) 
                    - dsspm * (            rsd(m,i-2,j,k)
                               - 4.0 * rsd(m,i-1,j,k)
                               + 6.0 * rsd(m,i,j,k)
                               - 4.0 * rsd(m,i+1,j,k)
                               +           rsd(m,i+2,j,k) );

            if (south == -1) {
             for (m=1; m<=5; m++) {
               frct(m,nx-2,j,k) = frct(m,nx-2,j,k)
                 - dsspm * (             rsd(m,nx-4,j,k)
                             - 4.0 * rsd(m,nx-3,j,k)
                             + 6.0 * rsd(m,nx-2,j,k)
                             - 4.0 * rsd(m,nx-1,j,k)  );
               frct(m,nx-1,j,k) = frct(m,nx-1,j,k)
                 - dsspm * (             rsd(m,nx-3,j,k)
                             - 4.0 * rsd(m,nx-2,j,k)
                             + 5.0 * rsd(m,nx-1,j,k) );
             }
            }

         }

//c---------------------------------------------------------------------
//c   eta-direction flux differences
//c---------------------------------------------------------------------
//c
//c   iex = flag : iex = 0  north/south communication
//c              : iex = 1  east/west communication
//c
//c---------------------------------------------------------------------
      iex   = 1;

//c---------------------------------------------------------------------
//c   communicate and receive/send two rows of data
//c---------------------------------------------------------------------

      exchange_3(rsd,iex);

      L1 = 0;
      if (west == -1) L1 = 1;
      L2 = ny + 1;
      if (east == -1) L2 = ny;

      jst1 = 1;
      jend1 = ny;
      if (west == -1) jst1 = 4;
      if (east == -1) jend1 = ny - 3;

      for (k=2; k<=nz-1; k++)
         for (j=L1; j<=L2; j++)
            for (i=ist; i<=iend; i++) {
               flux(1,i,j,k) = rsd(3,i,j,k);
               u31 = rsd(3,i,j,k) / rsd(1,i,j,k);
               q = 0.50 * (  rsd(2,i,j,k) * rsd(2,i,j,k)
                               + rsd(3,i,j,k) * rsd(3,i,j,k)
                               + rsd(4,i,j,k) * rsd(4,i,j,k) )
                            / rsd(1,i,j,k);
               flux(2,i,j,k) = rsd(2,i,j,k) * u31 ;
               flux(3,i,j,k) = rsd(3,i,j,k) * u31 + c2 * 
                             ( rsd(5,i,j,k) - q );
               flux(4,i,j,k) = rsd(4,i,j,k) * u31;
               flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u31;
      }

      for (k=2; k<=nz-1; k++) {
         for (i=ist; i<=iend; i++) 
            for (j=jst; j<=jend; j++) 
               for (m=1; m<=5; m++)
                  frct(m,i,j,k) =  frct(m,i,j,k)
                       - ty2 * ( flux(m,i,j+1,k) - flux(m,i,j-1,k) );

         for (j=jst; j<=L2; j++)
            for (i=ist; i<=iend; i++) {
               tmp = 1.0 / rsd(1,i,j,k);

               u21j = tmp * rsd(2,i,j,k);
               u31j = tmp * rsd(3,i,j,k);
               u41j = tmp * rsd(4,i,j,k);
               u51j = tmp * rsd(5,i,j,k);

               tmp = 1.0 / rsd(1,i,j-1,k);

               u21jm1 = tmp * rsd(2,i,j-1,k);
               u31jm1 = tmp * rsd(3,i,j-1,k);
               u41jm1 = tmp * rsd(4,i,j-1,k);
               u51jm1 = tmp * rsd(5,i,j-1,k);

               flux(2,i,j,k) = ty3 * ( u21j - u21jm1 );
               flux(3,i,j,k) = (4.0/3.0) * ty3 * 
                             ( u31j - u31jm1 );
               flux(4,i,j,k) = ty3 * ( u41j - u41jm1 );
               flux(5,i,j,k) = 0.50 * ( 1.0 - c1*c5 )
                    * ty3 * ( (   u21j*u21j +   u31j*u31j +   u41j*u41j )
                            - ( u21jm1*u21jm1 + u31jm1*u31jm1 + u41jm1*u41jm1 ) )
                    + (1.0/6.0)
                    * ty3 * ( u31j*u31j - u31jm1*u31jm1 )
                    + c1 * c5 * ty3 * ( u51j - u51jm1 );
         }

         for (j=jst; j<=jend; j++) 
            for (i=ist; i<=iend; i++) {
               frct(1,i,j,k) = frct(1,i,j,k)
                    + dy1 * ty1 * (            rsd(1,i,j-1,k)
                                   - 2.0 * rsd(1,i,j,k)
                                   +           rsd(1,i,j+1,k) );
               frct(2,i,j,k) = frct(2,i,j,k)
                + ty3 * c3 * c4 * ( flux(2,i,j+1,k) - flux(2,i,j,k) )
                    + dy2 * ty1 * (            rsd(2,i,j-1,k)
                                   - 2.0 * rsd(2,i,j,k)
                                   +           rsd(2,i,j+1,k) );
               frct(3,i,j,k) = frct(3,i,j,k)
                + ty3 * c3 * c4 * ( flux(3,i,j+1,k) - flux(3,i,j,k) )
                    + dy3 * ty1 * (            rsd(3,i,j-1,k)
                                   - 2.0 * rsd(3,i,j,k)
                                   +           rsd(3,i,j+1,k) );
               frct(4,i,j,k) = frct(4,i,j,k)
                + ty3 * c3 * c4 * ( flux(4,i,j+1,k) - flux(4,i,j,k) )
                    + dy4 * ty1 * (            rsd(4,i,j-1,k)
                                   - 2.0 * rsd(4,i,j,k)
                                   +           rsd(4,i,j+1,k) );
               frct(5,i,j,k) = frct(5,i,j,k)
                + ty3 * c3 * c4 * ( flux(5,i,j+1,k) - flux(5,i,j,k) )
                    + dy5 * ty1 * (            rsd(5,i,j-1,k)
                                   - 2.0 * rsd(5,i,j,k)
                                   +           rsd(5,i,j+1,k) );
         }

//c---------------------------------------------------------------------
//c   fourth-order dissipation
//c---------------------------------------------------------------------
         if (west == -1) {
            for (i=ist; i<=iend; i++) 
             for (m=1; m<=5; m++) {
               frct(m,i,2,k) = frct(m,i,2,k)
                 - dsspm * ( + 5.0 * rsd(m,i,2,k)
                             - 4.0 * rsd(m,i,3,k)
                             +           rsd(m,i,4,k) );
               frct(m,i,3,k) = frct(m,i,3,k)
                 - dsspm * ( - 4.0 * rsd(m,i,2,k)
                             + 6.0 * rsd(m,i,3,k)
                             - 4.0 * rsd(m,i,4,k)
                             +           rsd(m,i,5,k) );
            }
         }

         for (j=jst1; j<=jend1; j++)
            for (i=ist; i<=iend; i++) 
               for (m=1; m<=5; m++) 
                  frct(m,i,j,k) = frct(m,i,j,k)
                    - dsspm * (            rsd(m,i,j-2,k)
                              - 4.0 * rsd(m,i,j-1,k)
                              + 6.0 * rsd(m,i,j,k)
                              - 4.0 * rsd(m,i,j+1,k)
                              +           rsd(m,i,j+2,k) );

         if (east == -1) {
            for (i=ist; i<=iend; i++) 
             for (m=1; m<=5; m++) {
               frct(m,i,ny-2,k) = frct(m,i,ny-2,k)
                 - dsspm * (             rsd(m,i,ny-4,k)
                             - 4.0 * rsd(m,i,ny-3,k)
                             + 6.0 * rsd(m,i,ny-2,k)
                             - 4.0 * rsd(m,i,ny-1,k)  );
               frct(m,i,ny-1,k) = frct(m,i,ny-1,k)
                 - dsspm * (             rsd(m,i,ny-3,k)
                             - 4.0 * rsd(m,i,ny-2,k)
                             + 5.0 * rsd(m,i,ny-1,k)  );
             }
         }

      }

//c---------------------------------------------------------------------
//c   zeta-direction flux differences
//c---------------------------------------------------------------------
      for (k=1; k<=nz; k++) 
         for (j=jst; j<=jend; j++) 
            for (i=ist; i<=iend; i++) {
               flux(1,i,j,k) = rsd(4,i,j,k);
               u41 = rsd(4,i,j,k) / rsd(1,i,j,k);
               q = 0.50 * (  rsd(2,i,j,k) * rsd(2,i,j,k)
                               + rsd(3,i,j,k) * rsd(3,i,j,k)
                               + rsd(4,i,j,k) * rsd(4,i,j,k) )
                            / rsd(1,i,j,k);
               flux(2,i,j,k) = rsd(2,i,j,k) * u41 ;
               flux(3,i,j,k) = rsd(3,i,j,k) * u41 ;
               flux(4,i,j,k) = rsd(4,i,j,k) * u41 + c2 * 
                               ( rsd(5,i,j,k) - q );
               flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u41;
      }

      for (k=2; k<=nz-1; k++)
         for (j=jst; j<=jend; j++) 
            for (i=ist; i<=iend; i++) 
               for (m=1; m<=5; m++)
                  frct(m,i,j,k) =  frct(m,i,j,k)
                        - tz2 * ( flux(m,i,j,k+1) - flux(m,i,j,k-1) );

      for (k=2; k<=nz; k++) 
         for (j=jst; j<=jend; j++) 
            for (i=ist; i<=iend; i++) {
               tmp = 1.0 / rsd(1,i,j,k);

               u21k = tmp * rsd(2,i,j,k);
               u31k = tmp * rsd(3,i,j,k);
               u41k = tmp * rsd(4,i,j,k);
               u51k = tmp * rsd(5,i,j,k);

               tmp = 1.0 / rsd(1,i,j,k-1);

               u21km1 = tmp * rsd(2,i,j,k-1);
               u31km1 = tmp * rsd(3,i,j,k-1);
               u41km1 = tmp * rsd(4,i,j,k-1);
               u51km1 = tmp * rsd(5,i,j,k-1);

               flux(2,i,j,k) = tz3 * ( u21k - u21km1 );
               flux(3,i,j,k) = tz3 * ( u31k - u31km1 );
               flux(4,i,j,k) = (4.0/3.0) * tz3 * ( u41k 
                             - u41km1 );
               flux(5,i,j,k) = 0.50 * ( 1.0 - c1*c5 )
                    * tz3 * ( (   u21k*u21k +   u31k*u31k +   u41k*u41k )
                            - ( u21km1*u21km1 + u31km1*u31km1 + u41km1*u41km1 ) )
                    + (1.0/6.0)
                    * tz3 * ( u41k*u41k - u41km1*u41km1 )
                    + c1 * c5 * tz3 * ( u51k - u51km1 );
      }

      for (k=2; k<=nz-1; k++)
         for (j=jst; j<=jend; j++) 
            for (i=ist; i<=iend; i++) {
               frct(1,i,j,k) = frct(1,i,j,k)
                    + dz1 * tz1 * (            rsd(1,i,j,k+1)
                                   - 2.0 * rsd(1,i,j,k)
                                   +           rsd(1,i,j,k-1) );
               frct(2,i,j,k) = frct(2,i,j,k)
                + tz3 * c3 * c4 * ( flux(2,i,j,k+1) - flux(2,i,j,k) )
                    + dz2 * tz1 * (            rsd(2,i,j,k+1)
                                   - 2.0 * rsd(2,i,j,k)
                                   +           rsd(2,i,j,k-1) );
               frct(3,i,j,k) = frct(3,i,j,k)
                + tz3 * c3 * c4 * ( flux(3,i,j,k+1) - flux(3,i,j,k) )
                    + dz3 * tz1 * (            rsd(3,i,j,k+1)
                                   - 2.0 * rsd(3,i,j,k)
                                   +           rsd(3,i,j,k-1) );
               frct(4,i,j,k) = frct(4,i,j,k)
                + tz3 * c3 * c4 * ( flux(4,i,j,k+1) - flux(4,i,j,k) )
                    + dz4 * tz1 * (            rsd(4,i,j,k+1)
                                   - 2.0 * rsd(4,i,j,k)
                                   +           rsd(4,i,j,k-1) );
               frct(5,i,j,k) = frct(5,i,j,k)
                + tz3 * c3 * c4 * ( flux(5,i,j,k+1) - flux(5,i,j,k) )
                    + dz5 * tz1 * (            rsd(5,i,j,k+1)
                                   - 2.0 * rsd(5,i,j,k)
                                   +           rsd(5,i,j,k-1) );
      }

//c---------------------------------------------------------------------
//c   fourth-order dissipation
//c---------------------------------------------------------------------
      for (j=jst; j<=jend; j++) 
         for (i=ist; i<=iend; i++) 
            for (m=1; m<=5; m++) {
               frct(m,i,j,2) = frct(m,i,j,2)
                 - dsspm * ( + 5.0 * rsd(m,i,j,2)
                             - 4.0 * rsd(m,i,j,3)
                             +           rsd(m,i,j,4) );
               frct(m,i,j,3) = frct(m,i,j,3)
                 - dsspm * (- 4.0 * rsd(m,i,j,2)
                            + 6.0 * rsd(m,i,j,3)
                            - 4.0 * rsd(m,i,j,4)
                            +           rsd(m,i,j,5) );
      }

      for (k=4; k<=nz-3; k++)
         for (j=jst; j<=jend; j++) 
            for (i=ist; i<=iend; i++) 
               for (m=1; m<=5; m++) 
                  frct(m,i,j,k) = frct(m,i,j,k)
                    - dsspm * (           rsd(m,i,j,k-2)
                              - 4.0 * rsd(m,i,j,k-1)
                              + 6.0 * rsd(m,i,j,k)
                              - 4.0 * rsd(m,i,j,k+1)
                              +           rsd(m,i,j,k+2) );

      for (j=jst; j<=jend; j++) 
         for (i=ist; i<=iend; i++) 
            for (m=1; m<=5; m++) {
               frct(m,i,j,nz-2) = frct(m,i,j,nz-2)
                 - dsspm * (            rsd(m,i,j,nz-4)
                            - 4.0 * rsd(m,i,j,nz-3)
                            + 6.0 * rsd(m,i,j,nz-2)
                            - 4.0 * rsd(m,i,j,nz-1)  );
               frct(m,i,j,nz-1) = frct(m,i,j,nz-1)
                 - dsspm * (             rsd(m,i,j,nz-3)
                             - 4.0 * rsd(m,i,j,nz-2)
                             + 5.0 * rsd(m,i,j,nz-1)  );
      }

      return;
}