Distributed-Run-Time-Resource-Manager/RCCE_V2.0/apps/NPB/BT/.svn/text-base/rhs.c.svn-base

//---------------------------------------------------------------------
// 
// 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 "header.h"

void compute_rhs() {

//---------------------------------------------------------------------
//---------------------------------------------------------------------

      int c, i, j, k, m;
      double rho_inv, uijk, up1, um1, vijk, vp1, vm1,
           wijk, wp1, wm1;


//---------------------------------------------------------------------
//     loop over all cells owned by this node                           
//---------------------------------------------------------------------
      for (c = 1; c <= ncells; c++) {

//---------------------------------------------------------------------
//     compute the reciprocal of density, and the kinetic energy, 
//     and the speed of sound.
//---------------------------------------------------------------------
         for (k = -1; k <= cell_size(3,c); k++) {
            for (j = -1; j <= cell_size(2,c); j++) {
               for (i = -1; i <= cell_size(1,c); i++) {
                  rho_inv = 1.0e0/u(1,i,j,k,c);
                  rho_i(i,j,k,c) = rho_inv;
                  us(i,j,k,c) = u(2,i,j,k,c) * rho_inv;
                  vs(i,j,k,c) = u(3,i,j,k,c) * rho_inv;
                  ws(i,j,k,c) = u(4,i,j,k,c) * rho_inv;
                  square(i,j,k,c)     = 0.5e0* (
                       u(2,i,j,k,c)*u(2,i,j,k,c) + 
                       u(3,i,j,k,c)*u(3,i,j,k,c) +
                       u(4,i,j,k,c)*u(4,i,j,k,c) ) * rho_inv;
                  qs(i,j,k,c) = square(i,j,k,c) * rho_inv;
               }
            }
         }

//---------------------------------------------------------------------
// copy the exact forcing term to the right hand side;  because 
// this forcing term is known, we can store it on the whole of every 
// cell,  including the boundary                   
//---------------------------------------------------------------------

         for (k = 0; k <= cell_size(3,c)-1; k++) {
            for (j = 0; j <= cell_size(2,c)-1; j++) {
               for (i = 0; i <= cell_size(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = forcing(m,i,j,k,c);
                  }
               }
            }
         }


//---------------------------------------------------------------------
//     compute xi-direction fluxes 
//---------------------------------------------------------------------
         for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) {
            for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  uijk = us(i,j,k,c);
                  up1  = us(i+1,j,k,c);
                  um1  = us(i-1,j,k,c);

                  rhs(1,i,j,k,c) = rhs(1,i,j,k,c) + dx1tx1 * 
                       (u(1,i+1,j,k,c) - 2.0e0*u(1,i,j,k,c) + 
                       u(1,i-1,j,k,c)) -
                       tx2 * (u(2,i+1,j,k,c) - u(2,i-1,j,k,c));

                  rhs(2,i,j,k,c) = rhs(2,i,j,k,c) + dx2tx1 * 
                       (u(2,i+1,j,k,c) - 2.0e0*u(2,i,j,k,c) + 
                       u(2,i-1,j,k,c)) +
                       xxcon2*con43 * (up1 - 2.0e0*uijk + um1) -
                       tx2 * (u(2,i+1,j,k,c)*up1 - 
                       u(2,i-1,j,k,c)*um1 +
                       (u(5,i+1,j,k,c)- square(i+1,j,k,c)-
                       u(5,i-1,j,k,c)+ square(i-1,j,k,c))*
                       c2);

                  rhs(3,i,j,k,c) = rhs(3,i,j,k,c) + dx3tx1 * 
                       (u(3,i+1,j,k,c) - 2.0e0*u(3,i,j,k,c) +
                       u(3,i-1,j,k,c)) +
                       xxcon2 * (vs(i+1,j,k,c) - 2.0e0*vs(i,j,k,c) +
                       vs(i-1,j,k,c)) -
                       tx2 * (u(3,i+1,j,k,c)*up1 - 
                       u(3,i-1,j,k,c)*um1);

                  rhs(4,i,j,k,c) = rhs(4,i,j,k,c) + dx4tx1 * 
                       (u(4,i+1,j,k,c) - 2.0e0*u(4,i,j,k,c) +
                       u(4,i-1,j,k,c)) +
                       xxcon2 * (ws(i+1,j,k,c) - 2.0e0*ws(i,j,k,c) +
                       ws(i-1,j,k,c)) -
                       tx2 * (u(4,i+1,j,k,c)*up1 - 
                       u(4,i-1,j,k,c)*um1);

                  rhs(5,i,j,k,c) = rhs(5,i,j,k,c) + dx5tx1 * 
                       (u(5,i+1,j,k,c) - 2.0e0*u(5,i,j,k,c) +
                       u(5,i-1,j,k,c)) +
                       xxcon3 * (qs(i+1,j,k,c) - 2.0e0*qs(i,j,k,c) +
                       qs(i-1,j,k,c)) +
                       xxcon4 * (up1*up1 -       2.0e0*uijk*uijk + 
                       um1*um1) +
                       xxcon5 * (u(5,i+1,j,k,c)*rho_i(i+1,j,k,c) - 
                       2.0e0*u(5,i,j,k,c)*rho_i(i,j,k,c) +
                       u(5,i-1,j,k,c)*rho_i(i-1,j,k,c)) -
                       tx2 * ( (c1*u(5,i+1,j,k,c) - 
                       c2*square(i+1,j,k,c))*up1 -
                       (c1*u(5,i-1,j,k,c) - 
                       c2*square(i-1,j,k,c))*um1 );
               }
            }
         }

//---------------------------------------------------------------------
//     add fourth order xi-direction dissipation               
//---------------------------------------------------------------------
         if (start(1,c) > 0) {
            for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) {
               for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
                  i = 1;
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c)- dssp * 
                          ( 5.0e0*u(m,i,j,k,c) - 4.0e0*u(m,i+1,j,k,c) +
                          u(m,i+2,j,k,c));
                  }

                  i = 2;
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp * 
                          (-4.0e0*u(m,i-1,j,k,c) + 6.0e0*u(m,i,j,k,c) -
                          4.0e0*u(m,i+1,j,k,c) + u(m,i+2,j,k,c));
                  }
               }
            }
         }

         for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) {
            for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
               for (i = 3*start(1,c); i <= cell_size(1,c)-3*end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp * 
                          (  u(m,i-2,j,k,c) - 4.0e0*u(m,i-1,j,k,c) + 
                          6.0*u(m,i,j,k,c) - 4.0e0*u(m,i+1,j,k,c) + 
                          u(m,i+2,j,k,c) );
                  }
               }
            }
         }
         

         if (end(1,c) > 0) {
            for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) {
               for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
                  i = cell_size(1,c)-3;
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
                          ( u(m,i-2,j,k,c) - 4.0e0*u(m,i-1,j,k,c) + 
                          6.0e0*u(m,i,j,k,c) - 4.0e0*u(m,i+1,j,k,c) );
                  }

                  i = cell_size(1,c)-2;
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
                          ( u(m,i-2,j,k,c) - 4.e0*u(m,i-1,j,k,c) +
                          5.e0*u(m,i,j,k,c) );
                  }
               }
            }
         }

//---------------------------------------------------------------------
//     compute eta-direction fluxes 
//---------------------------------------------------------------------
         for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) {
            for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  vijk = vs(i,j,k,c);
                  vp1  = vs(i,j+1,k,c);
                  vm1  = vs(i,j-1,k,c);
                  rhs(1,i,j,k,c) = rhs(1,i,j,k,c) + dy1ty1 * 
                       (u(1,i,j+1,k,c) - 2.0e0*u(1,i,j,k,c) + 
                       u(1,i,j-1,k,c)) -
                       ty2 * (u(3,i,j+1,k,c) - u(3,i,j-1,k,c));
                  rhs(2,i,j,k,c) = rhs(2,i,j,k,c) + dy2ty1 * 
                       (u(2,i,j+1,k,c) - 2.0e0*u(2,i,j,k,c) + 
                       u(2,i,j-1,k,c)) +
                       yycon2 * (us(i,j+1,k,c) - 2.0e0*us(i,j,k,c) + 
                       us(i,j-1,k,c)) -
                       ty2 * (u(2,i,j+1,k,c)*vp1 - 
                       u(2,i,j-1,k,c)*vm1);
                  rhs(3,i,j,k,c) = rhs(3,i,j,k,c) + dy3ty1 * 
                       (u(3,i,j+1,k,c) - 2.0e0*u(3,i,j,k,c) + 
                       u(3,i,j-1,k,c)) +
                       yycon2*con43 * (vp1 - 2.0e0*vijk + vm1) -
                       ty2 * (u(3,i,j+1,k,c)*vp1 - 
                       u(3,i,j-1,k,c)*vm1 +
                       (u(5,i,j+1,k,c) - square(i,j+1,k,c) - 
                       u(5,i,j-1,k,c) + square(i,j-1,k,c))
                       *c2);
                  rhs(4,i,j,k,c) = rhs(4,i,j,k,c) + dy4ty1 * 
                       (u(4,i,j+1,k,c) - 2.0e0*u(4,i,j,k,c) + 
                       u(4,i,j-1,k,c)) +
                       yycon2 * (ws(i,j+1,k,c) - 2.0e0*ws(i,j,k,c) + 
                       ws(i,j-1,k,c)) -
                       ty2 * (u(4,i,j+1,k,c)*vp1 - 
                       u(4,i,j-1,k,c)*vm1);
                  rhs(5,i,j,k,c) = rhs(5,i,j,k,c) + dy5ty1 * 
                       (u(5,i,j+1,k,c) - 2.0e0*u(5,i,j,k,c) + 
                       u(5,i,j-1,k,c)) +
                       yycon3 * (qs(i,j+1,k,c) - 2.0e0*qs(i,j,k,c) + 
                       qs(i,j-1,k,c)) +
                       yycon4 * (vp1*vp1       - 2.0e0*vijk*vijk + 
                       vm1*vm1) +
                       yycon5 * (u(5,i,j+1,k,c)*rho_i(i,j+1,k,c) - 
                       2.0e0*u(5,i,j,k,c)*rho_i(i,j,k,c) +
                       u(5,i,j-1,k,c)*rho_i(i,j-1,k,c)) -
                       ty2 * ((c1*u(5,i,j+1,k,c) - 
                       c2*square(i,j+1,k,c)) * vp1 -
                       (c1*u(5,i,j-1,k,c) - 
                       c2*square(i,j-1,k,c)) * vm1);
               }
            }
         }

//---------------------------------------------------------------------
//     add fourth order eta-direction dissipation         
//---------------------------------------------------------------------
         if (start(2,c) > 0) {
            for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) {
               j = 1;
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c)- dssp * 
                          ( 5.0e0*u(m,i,j,k,c) - 4.0e0*u(m,i,j+1,k,c) +
                          u(m,i,j+2,k,c));
                  }
               }

               j = 2;
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp * 
                          (-4.0e0*u(m,i,j-1,k,c) + 6.0e0*u(m,i,j,k,c) -
                          4.0e0*u(m,i,j+1,k,c) + u(m,i,j+2,k,c));
                  }
               }
            }
         }

         for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) {
            for (j = 3*start(2,c); j <= cell_size(2,c)-3*end(2,c)-1; j++) {
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp * 
                          (  u(m,i,j-2,k,c) - 4.0e0*u(m,i,j-1,k,c) + 
                          6.0*u(m,i,j,k,c) - 4.0e0*u(m,i,j+1,k,c) + 
                          u(m,i,j+2,k,c) );
                  }
               }
            }
         }
         
         if (end(2,c) > 0) {
            for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) {
               j = cell_size(2,c)-3;
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
                          ( u(m,i,j-2,k,c) - 4.0e0*u(m,i,j-1,k,c) + 
                          6.0e0*u(m,i,j,k,c) - 4.0e0*u(m,i,j+1,k,c) );
                  }
               }

               j = cell_size(2,c)-2;
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
                          ( u(m,i,j-2,k,c) - 4.e0*u(m,i,j-1,k,c) +
                          5.e0*u(m,i,j,k,c) );
                  }
               }
            }
         }

//---------------------------------------------------------------------
//     compute zeta-direction fluxes 
//---------------------------------------------------------------------
         for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) {
            for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  wijk = ws(i,j,k,c);
                  wp1  = ws(i,j,k+1,c);
                  wm1  = ws(i,j,k-1,c);

                  rhs(1,i,j,k,c) = rhs(1,i,j,k,c) + dz1tz1 * 
                       (u(1,i,j,k+1,c) - 2.0e0*u(1,i,j,k,c) + 
                       u(1,i,j,k-1,c)) -
                       tz2 * (u(4,i,j,k+1,c) - u(4,i,j,k-1,c));
                  rhs(2,i,j,k,c) = rhs(2,i,j,k,c) + dz2tz1 * 
                       (u(2,i,j,k+1,c) - 2.0e0*u(2,i,j,k,c) + 
                       u(2,i,j,k-1,c)) +
                       zzcon2 * (us(i,j,k+1,c) - 2.0e0*us(i,j,k,c) + 
                       us(i,j,k-1,c)) -
                       tz2 * (u(2,i,j,k+1,c)*wp1 - 
                       u(2,i,j,k-1,c)*wm1);
                  rhs(3,i,j,k,c) = rhs(3,i,j,k,c) + dz3tz1 * 
                       (u(3,i,j,k+1,c) - 2.0e0*u(3,i,j,k,c) + 
                       u(3,i,j,k-1,c)) +
                       zzcon2 * (vs(i,j,k+1,c) - 2.0e0*vs(i,j,k,c) + 
                       vs(i,j,k-1,c)) -
                       tz2 * (u(3,i,j,k+1,c)*wp1 - 
                       u(3,i,j,k-1,c)*wm1);
                  rhs(4,i,j,k,c) = rhs(4,i,j,k,c) + dz4tz1 * 
                       (u(4,i,j,k+1,c) - 2.0e0*u(4,i,j,k,c) + 
                       u(4,i,j,k-1,c)) +
                       zzcon2*con43 * (wp1 - 2.0e0*wijk + wm1) -
                       tz2 * (u(4,i,j,k+1,c)*wp1 - 
                       u(4,i,j,k-1,c)*wm1 +
                       (u(5,i,j,k+1,c) - square(i,j,k+1,c) - 
                       u(5,i,j,k-1,c) + square(i,j,k-1,c))
                       *c2);
                  rhs(5,i,j,k,c) = rhs(5,i,j,k,c) + dz5tz1 * 
                       (u(5,i,j,k+1,c) - 2.0e0*u(5,i,j,k,c) + 
                       u(5,i,j,k-1,c)) +
                       zzcon3 * (qs(i,j,k+1,c) - 2.0e0*qs(i,j,k,c) + 
                       qs(i,j,k-1,c)) +
                       zzcon4 * (wp1*wp1 - 2.0e0*wijk*wijk + 
                       wm1*wm1) +
                       zzcon5 * (u(5,i,j,k+1,c)*rho_i(i,j,k+1,c) - 
                       2.0e0*u(5,i,j,k,c)*rho_i(i,j,k,c) +
                       u(5,i,j,k-1,c)*rho_i(i,j,k-1,c)) -
                       tz2 * ( (c1*u(5,i,j,k+1,c) - 
                       c2*square(i,j,k+1,c))*wp1 -
                       (c1*u(5,i,j,k-1,c) - 
                       c2*square(i,j,k-1,c))*wm1);
               }
            }
         }

//---------------------------------------------------------------------
//     add fourth order zeta-direction dissipation                
//---------------------------------------------------------------------
         if (start(3,c) > 0) {
            k = 1;
            for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c)- dssp * 
                          ( 5.0e0*u(m,i,j,k,c) - 4.0e0*u(m,i,j,k+1,c) +
                          u(m,i,j,k+2,c));
                  }
               }
            }

            k = 2;
            for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp * 
                          (-4.0e0*u(m,i,j,k-1,c) + 6.0e0*u(m,i,j,k,c) -
                          4.0e0*u(m,i,j,k+1,c) + u(m,i,j,k+2,c));
                  }
               }
            }
         }

         for (k = 3*start(3,c); k <= cell_size(3,c)-3*end(3,c)-1; k++) {
            for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp * 
                          (  u(m,i,j,k-2,c) - 4.0e0*u(m,i,j,k-1,c) + 
                          6.0*u(m,i,j,k,c) - 4.0e0*u(m,i,j,k+1,c) + 
                          u(m,i,j,k+2,c) );
                  }
               }
            }
         }
         
         if (end(3,c) > 0) {
            k = cell_size(3,c)-3;
            for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
                          ( u(m,i,j,k-2,c) - 4.0e0*u(m,i,j,k-1,c) + 
                          6.0e0*u(m,i,j,k,c) - 4.0e0*u(m,i,j,k+1,c) );
                  }
               }
            }

            k = cell_size(3,c)-2;
            for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) - dssp *
                          ( u(m,i,j,k-2,c) - 4.e0*u(m,i,j,k-1,c) +
                          5.e0*u(m,i,j,k,c) );
                  }
               }
            }
         }

         for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) {
            for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) {
               for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) {
                  for (m = 1; m <= 5; m++) {
                     rhs(m,i,j,k,c) = rhs(m,i,j,k,c) * dt;
                  }
               }
            }
         }

      }
      
      return;
}