src/expde/evpar/evpar.cc

Go to the documentation of this file.

//    expde: expression templates for partial differential equations.
//    Copyright (C) 2001  Christoph Pflaum
//    This program is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation; either version 2 of the License, or
//    (at your option) any later version.
//
//    This program is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//                 SEE  Notice1.doc made by 
//                 LAWRENCE LIVERMORE NATIONAL LABORATORY
//

// ------------------------------------------------------------
//
// evpar.cc
//
// ------------------------------------------------------------


// Include level 0:

#ifdef COMP_GNUOLD
 #include <iostream.h>
 #include <fstream.h>
 #include <math.h>
 #include <time.h>
#else
 #include <iostream>
 #include <fstream>
 #include <cmath>
 #include <ctime>
#endif 


// Include level 0:
#include "../parser.h"

// Include level 1:
#include "../paramete.h"
#include "../abbrevi.h"
#include "../math_lib/math_lib.h"

// Include level 2:
#include "../basic/basic.h"

// Include level 3:
#include "../domain/domain.h"

// Include level 4:
#include "../formulas/boundy.h"
#include "../formulas/loc_sten.h"


// Include level 5.0
#include "../grid/gpar.h"
#include "../grid/parallel.h"
#include "../grid/mgcoeff.h"
#include "../grid/sto_man.h"
#include "../grid/gridbase.h"
#include "../grid/grid.h"

// Include level 5.1
#include "evpar.h"



Evaluation_Parallelization_object::Evaluation_Parallelization_object(
                                   Grid *g, int max_v) : 
  grid(g) {
  //  max_v is maximal number of variables
  max_var = max_v;
  max_level  = g->Max_level();
  storage_of_update_type = new type_of_update[max_level*(max_var+1)];
  for(int i=max_level*(max_var+1)-1;i>=0;--i) {
    storage_of_update_type[i] = full_update;
  }
  array_vars = new bool[max_var];
  array_num  = new int[max_var];
  array_vars_other = new bool[max_var];
  array_num_other  = new int[max_var];
  // this formula is connected with the formula
  // dir_sons reordering_formula(dir_sons d);

  Parallel_GS_formula[0] = Edir;
  Parallel_GS_formula[1] = Ndir;
  Parallel_GS_formula[2] = Tdir;
  Parallel_GS_formula[3] = Sdir;

  Parallel_GS_formula[4] = Wdir;
  Parallel_GS_formula[5] = Ndir;
  Parallel_GS_formula[6] = Ddir;
  Parallel_GS_formula[7] = Sdir;
};

void Evaluation_Parallelization_object::Start_evaluation_seriel(
                const Run_type_object rob, const int l, const int ss,
                const int num_ops) {
  size_stencil = ss;
  run_object = rob;
  level=l;
  status=0;
  number_ops=num_ops;
}

void Evaluation_Parallelization_object::StartA_varcon_evaluation(int l) {
  level=l;
  for(int i=0;i<max_var;++i) {
    array_vars[i]=false;
  }
};


void Evaluation_Parallelization_object::StartB_varcon_evaluation() {
  type_of_update type_upd;

  anz_num=0;

  for(int i=0;i<max_var;++i) {
    if(array_vars[i]) {
      type_upd = Give_type_of_update(i,level);
      if(type_upd!=full_update) {
        array_num[anz_num]=i;
        ++anz_num;
      }
      Set_type_of_update(i,level,full_update);
    }
  }
  grid->Update_Variable_no_full(array_num,anz_num,level);
};

void Evaluation_Parallelization_object::StartA_evaluation(
                const Run_type_object rob, const int l, const int ss,
                const int num_ops) {
  size_stencil = ss;
  run_object = rob;
  level=l;
  status=0;
  number_ops=num_ops;

  if(size_stencil>1) {
    other_is_fine=other_is_coarse=false;
    for(int i=0;i<max_var;++i) {
      array_vars[i]=array_vars_other[i]=false;
    }
    //    grid->Full_update_Variable();
  }
}


void Evaluation_Parallelization_object::StartB_evaluation(
                                        const int number_variable,
                                        const bool GS_expr) {
  anz_num=anz_num_other=0;

  /*
  //test
  GS_expression = false;
  if(size_stencil>1) {
    // for left hand side:
    for(int i=0;i<max_var;++i) {
      if(array_vars[i]) {
        if(i==number_variable)
          GS_expression = true;
      }
    }
  }
  if(GS_expression!=GS_expr) {
    if(grid->Give_my_rank()==0)
      cout << "\n error! " << GS_expr << endl;
  }
  */


  GS_expression = GS_expr;

  type_of_update type_upd;

  if(size_stencil>1) {
    // coarse and fine part
    update_other_bla = false;
    for(int i=0;i<max_var;++i) {
      if(array_vars_other[i]) {
        if(other_is_fine) {
          type_upd = Give_type_of_update(i,level+1);
          if(type_upd!=full_update) {
            array_num_other[anz_num_other]=i;
            ++anz_num_other;
            if(type_upd==no_update) {
              update_other_bla = true;
            }
            Set_type_of_update(i,level+1,full_update);
          }
        }
        if(other_is_coarse) {
          type_upd = Give_type_of_update(i,level-1);
          if(type_upd!=full_update) {
            array_num_other[anz_num_other]=i;
            ++anz_num_other;
            if(type_upd==no_update) {
              update_other_bla = true;
            }
            Set_type_of_update(i,level-1,full_update);
          }
        }
      }
    }
    if(GS_expr==true) {
      // Gauss Seidel
      update_ENT = false;
      update_WSD = false;
      // for left hand side:
      for(int i=0;i<max_var;++i) {
        if(array_vars[i]) {
          type_upd = Give_type_of_update(i,level);
          array_num[anz_num]=i;
          ++anz_num;
          if(type_upd==no_update) {
            update_ENT = true;
          }
          //old
          //      Set_type_of_update(i,level,partial_update);
          //new
          Set_type_of_update(i,level,full_update);
        }
      }
      // for left hand side:
      Set_type_of_update(number_variable,level,partial_update);
      left_level  = level;
      left_update = partial_update;
    }
    else {
      // Jacobi
      update_ENT = false;
      update_WSD = true;
      // for right hand side:
      for(int i=0;i<max_var;++i) {
        if(array_vars[i]) {
          type_upd = Give_type_of_update(i,level);
          if(type_upd!=full_update) {
            array_num[anz_num]=i;
            ++anz_num;
            if(type_upd==no_update) {
              update_ENT = true;
            }
          }
          //old
          //      Set_type_of_update(i,level,no_update);
          //new
          Set_type_of_update(i,level,full_update);
        }
      }
      // for left hand side:
      Set_type_of_update(number_variable,level,no_update);
      left_level  = level;
      left_update = no_update;
    }
    
    if(other_is_fine && other_is_coarse) 
      cout << " Restriction and Prolongation in one expression "
           << " is not allowed!! " << endl;
  }
  else {
    // for left hand side:
    Set_type_of_update(number_variable,level,no_update);
    left_level  = level;
    left_update = no_update;
  }
    /*
    grid->Update_Variable(Wdir,array_num,anz_num,level);
    grid->Update_Variable(Edir,array_num,anz_num,level);
    grid->Update_Variable(Ndir,array_num,anz_num,level);
    grid->Update_Variable(Sdir,array_num,anz_num,level);
    grid->Update_Variable(Tdir,array_num,anz_num,level);
    grid->Update_Variable(Ddir,array_num,anz_num,level);
    */
}


evaluation_typ Evaluation_Parallelization_object::Give_next_evaluation_typ() {
  // -> parallel
  if(parallel_version) {
    if(status==0 && size_stencil>1) {
      if(other_is_fine) {
        if(update_other_bla==true) {
          // update from no -> full
          grid->Update_Variable_no_full(array_num_other,anz_num_other,level+1);
        }
        else {
          // test
            grid->Update_Variable_no_full(array_num,anz_num,level);

          // update from partial -> full
          grid->Update_Variable(Edir,array_num_other,anz_num_other,level+1);
          grid->Update_Variable(Ndir,array_num_other,anz_num_other,level+1);
          grid->Update_Variable(Tdir,array_num_other,anz_num_other,level+1);
        }
      }
      if(other_is_coarse) {
       if((level-1)>=(grid->Give_n_parallel()-1)) {
         if(update_other_bla==true) {      
           // update from no -> partial
 grid->Update_Variable_no_full(array_num_other,anz_num_other,level-1);
           /*
           grid->Update_Variable(Wdir,array_num_other,anz_num_other,level-1);
           grid->Update_Variable(Sdir,array_num_other,anz_num_other,level-1);
           grid->Update_Variable(Ddir,array_num_other,anz_num_other,level-1);
           */
         }
        }
        else {
          grid->Update_Variable_for_prolongation(array_num_other,
                                                 anz_num_other,level-1);

          /*
          grid->Update_Variable_for_prolongation(Wdir,array_num_other,
                                                 anz_num_other,level-1);
          grid->Update_Variable_for_prolongation(Edir,array_num_other,
                                                 anz_num_other,level-1);
          grid->Update_Variable_for_prolongation(Ndir,array_num_other,
                                                 anz_num_other,level-1);
          grid->Update_Variable_for_prolongation(Sdir,array_num_other,
                                                 anz_num_other,level-1);
          grid->Update_Variable_for_prolongation(Tdir,array_num_other,
                                                 anz_num_other,level-1);
          grid->Update_Variable_for_prolongation(Ddir,array_num_other,
                                                 anz_num_other,level-1);
          */
        }
      }
    }
    if(GS_expression==false) {
      if(status==0) {// FOR JACOBI
        if(size_stencil>1 && anz_num>0) {
          if(update_ENT &&  update_WSD) {
            // update from no -> full
            grid->Update_Variable_no_full(array_num,anz_num,level);
          }
          else {
            // update from partial -> full

            // test
            grid->Update_Variable_no_full(array_num,anz_num,level);

            grid->Update_Variable(Edir,array_num,anz_num,level);
            grid->Update_Variable(Ndir,array_num,anz_num,level);
            grid->Update_Variable(Tdir,array_num,anz_num,level);
          }
        }
        ++status;
        if(run_object.r_nearb==true) {
          cellpoi_pointer = grid->Start_P_cellpoi(level);
          
          if(size_stencil>1)  return cellf_large_evaluation; 
          else                return cellf_1_evaluation; 
        }
      }
      int color, kind;
      //      while(status<25) {
      while(status<33) {
        // color:   0,1, ..., 7
        // kind:    0,1,2
        // beachte: 3*8  = 24 
        //          24+1 = 25
        // beachte: 4*8  = 32 
        //          32+1 = 33
        color = (status-1)/4;
        kind  = (status-1)%4;

        if(kind==0) {
          ++status;
          if(run_object.r_interior==true) {
            interior_pointer = grid->Start_P_interior(level,color);
            if(max_level==level) {
              if(size_stencil==25) return interior_25_evaluation_fine;
              if(size_stencil==17) return interior_17_evaluation_fine;
              if(size_stencil==15) return interior_15_evaluation_fine;
              else                 return interior_1_evaluation_fine; 
            }
            else {
              if(size_stencil==25) return interior_25_evaluation_coarse;
              if(size_stencil==17) return interior_17_evaluation_coarse;
              if(size_stencil==15) return interior_15_evaluation_coarse;
              else                 return interior_1_evaluation_coarse;  
            }
          }
          kind=1;
        }
        if(kind==1) {
          ++status;
          if(run_object.r_nearb==true) {
            nearb_pointer = grid->Start_P_nearb(level,color);
            
            if(size_stencil>=17) return nearb_17_evaluation;
            if(size_stencil==15) return nearb_15_evaluation;
            else                 return nearb_1_evaluation;
          }
          kind=2;
        }
        if(kind==2) {
          ++status;
          Bo2p_pointer = grid->Start_P_Bo2p(level,color);

          if(run_object.r_boundary==1) {
            return boundary_evaluation_all; 
          }
          if(run_object.r_boundary<0) {
            return boundary_evaluation_subset; 
          }
        }
        if(kind==3) {
          ++status;
          nearb_pointer = grid->Start_P_exteri(level,color);
          if(run_object.r_boundary==1) {
            if(size_stencil>=17) return nearb_17_evaluation;
            if(size_stencil==15) return nearb_15_evaluation;
            else                 return nearb_1_evaluation;
          }
          if(run_object.r_boundary<0) {
            if(size_stencil>=17) return exteri_17_evaluation_subset;
            if(size_stencil==15) return exteri_15_evaluation_subset;
            else                 return exteri_1_evaluation_subset;
          }
        }
      }
      return stop_evaluation; 
    }
    else { // FOR GAUSS-SEIDEL
      if(status==0) {
        ++status;
        if(size_stencil>1 && anz_num>0) {
          if(update_ENT) {
          // test
            grid->Update_Variable_no_full(array_num,anz_num,level);

            // update from no -> partial
            grid->Update_Variable(Wdir,array_num,anz_num,level);
            grid->Update_Variable(Sdir,array_num,anz_num,level);
            grid->Update_Variable(Ddir,array_num,anz_num,level);
          }
        }
        if(run_object.r_nearb==true) {
          cellpoi_pointer = grid->Start_P_cellpoi(level);
          
          if(size_stencil>1)  return cellf_large_evaluation; 
          else                return cellf_1_evaluation; 
        }
      }
      
      int color, kind;
      //      while(status<25) {
      while(status<33) {
        // color:   0,1, ..., 7
        // kind:    0,1,2
        // beachte: 3*8  = 24 
        //          24+1 = 25
        // beachte: 4*8  = 32 
        //          32+1 = 33
        color = (status-1)/4;
        kind  = (status-1)%4;
        
        if(size_stencil>1 && color!=0 && kind==0 && anz_num>0) {

          // test
            grid->Update_Variable_no_full(array_num,anz_num,level);

          /*
          if(color>=10) {
            grid->Update_Variable(Wdir,array_num,anz_num,level);
            grid->Update_Variable(Edir,array_num,anz_num,level);
            grid->Update_Variable(Ndir,array_num,anz_num,level);
            grid->Update_Variable(Sdir,array_num,anz_num,level);
            grid->Update_Variable(Tdir,array_num,anz_num,level);
            grid->Update_Variable(Ddir,array_num,anz_num,level);
          }
          */
          // Eigentlich sollte man folgendes nehmen:
          grid->Update_Variable(Parallel_GS_formula[color-1],array_num,
                                anz_num,level);
          
        }
        
        if(kind==0) {
          ++status;
          if(run_object.r_interior==true) {
            interior_pointer = grid->Start_P_interior(level,color);
            
            if(max_level==level) {
              if(size_stencil==25) return interior_25_evaluation_fine;
              if(size_stencil==17) return interior_17_evaluation_fine;
              if(size_stencil==15) return interior_15_evaluation_fine;
              else                 return interior_1_evaluation_fine; 
            }
            else {
              if(size_stencil==25) return interior_25_evaluation_coarse;
              if(size_stencil==17) return interior_17_evaluation_coarse;
              if(size_stencil==15) return interior_15_evaluation_coarse;
              else                 return interior_1_evaluation_coarse;  
            }
          }
          kind=1;
        }
        if(kind==1) {
          ++status;
          if(run_object.r_nearb==true) {
            nearb_pointer = grid->Start_P_nearb(level,color);
            
            if(size_stencil>=17) return nearb_17_evaluation;
            if(size_stencil==15) return nearb_15_evaluation;
            else                 return nearb_1_evaluation;
          }
          kind=2;
        }
        if(kind==2) {
          ++status;
          Bo2p_pointer = grid->Start_P_Bo2p(level,color);
          
          if(run_object.r_boundary==1) {
            return boundary_evaluation_all; 
          }
          if(run_object.r_boundary<0) {
            return boundary_evaluation_subset; 
          }
        }
        if(kind==3) {
          ++status;
          nearb_pointer = grid->Start_P_exteri(level,color);
          if(run_object.r_boundary==1) {
            if(size_stencil>=17) return nearb_17_evaluation;
            if(size_stencil==15) return nearb_15_evaluation;
            else                 return nearb_1_evaluation;
          }
          if(run_object.r_boundary<0) {
            if(size_stencil>=17) return exteri_17_evaluation_subset;
            if(size_stencil==15) return exteri_15_evaluation_subset;
            else                 return exteri_1_evaluation_subset;
          }
        }
      }
      if(size_stencil>1) {
        grid->Update_Variable(Parallel_GS_formula[7],array_num,anz_num,level);
      }
      return stop_evaluation; 
    }
  }
  // serial
  else {
    if(status==0) {
      ++status;
      if(run_object.r_nearb==true) {
        cellpoi_pointer = grid->Start_P_cellpoi(level);
        
        if(size_stencil>1)  return cellf_large_evaluation; 
        else                return cellf_1_evaluation; 
      }
    }
    int color, kind;
    //      while(status<25) {
    while(status<33) {
      // color:   0,1, ..., 7
      // kind:    0,1,2
      // beachte: 3*8  = 24 
      //          24+1 = 25
      // beachte: 4*8  = 32 
      //          32+1 = 33
      color = (status-1)/4;
      kind  = (status-1)%4;
      
      if(kind==0) {
        ++status;
        if(run_object.r_interior==true) {
          interior_pointer = grid->Start_P_interior(level,color);
          if(max_level==level) {
            if(size_stencil==25) return interior_25_evaluation_fine;
            if(size_stencil==17) return interior_17_evaluation_fine;
            if(size_stencil==15) return interior_15_evaluation_fine;
            else                 return interior_1_evaluation_fine; 
          }
          else {
            if(size_stencil==25) return interior_25_evaluation_coarse;
            if(size_stencil==17) return interior_17_evaluation_coarse;
            if(size_stencil==15) return interior_15_evaluation_coarse;
            else                 return interior_1_evaluation_coarse;  
          }
        }
        kind=1;
      }
      if(kind==1) {
        ++status;
        if(run_object.r_nearb==true) {
          nearb_pointer = grid->Start_P_nearb(level,color);
          
          if(size_stencil>=17) return nearb_17_evaluation;
          if(size_stencil==15) return nearb_15_evaluation;
          else                 return nearb_1_evaluation;
        }
        kind=2;
      }
      if(kind==2) {
        ++status;
        Bo2p_pointer = grid->Start_P_Bo2p(level,color);
        
        if(run_object.r_boundary==1) {
          return boundary_evaluation_all; 
        }
        if(run_object.r_boundary<0) {
          return boundary_evaluation_subset; 
        }
      }
      if(kind==3) {
        ++status;
        nearb_pointer = grid->Start_P_exteri(level,color);
        if(run_object.r_boundary==1) {
          if(size_stencil>=17) return nearb_17_evaluation;
          if(size_stencil==15) return nearb_15_evaluation;
          else                 return nearb_1_evaluation;
        }
        if(run_object.r_boundary<0) {
          if(size_stencil>=17) return exteri_17_evaluation_subset;
          if(size_stencil==15) return exteri_15_evaluation_subset;
          else                 return exteri_1_evaluation_subset;
        }
      }
    }
    return stop_evaluation; 
  }
};

---