From 2fbc6d4c5d0b72ef9a31bfb67c98c56d3063a0bb Mon Sep 17 00:00:00 2001
From: iwia06 <christoph.pflaum@fau.de>
Date: Fri, 22 Jan 2021 12:12:17 +0100
Subject: [PATCH] performance Optimierung
---
program/main.cc | 17 ++-
program/program.pro | 2 +-
program/source/extemp/variable.cc | 18 +++
program/source/extemp/variable.h | 131 ++++++++++++++++-
program/source/extemp/variable_cc.h | 161 +++++++++++++++++++++
program/source/interpol/interpolCellVar.cc | 45 ++++++
program/source/interpol/interpolCellVar.h | 31 +++-
7 files changed, 399 insertions(+), 6 deletions(-)
diff --git a/program/main.cc b/program/main.cc
index 6e43d53..2ac9d1f 100644
--- a/program/main.cc
+++ b/program/main.cc
@@ -98,12 +98,25 @@ int main(int argc, char** argv)
u = Cos(X*M_PI)*Cos(Y*M_PI)*Cos(Z*M_PI);
+ /*
Interpolate_from_cell_to_point<Stencil_constant_z> interpolGrid(iteration);
+ interpolGrid.init(&blockgrid);
+ interpolGrid.doInterpolation(&u_cell,&interpol);
+*/
- interpolGrid.init(&blockgrid);
+ FastInterpolate_from_cell_to_point interpolGrid;
+ interpolGrid.init(&blockgrid);
+ interpolGrid.doInterpolation(&u_cell,&interpol);
- interpolGrid.doInterpolation(&u_cell,&interpol);
+
+ std::ofstream DATEIA;
+ DATEIA.open("interpol.vtk");
+ interpol.Print_VTK(DATEIA);
+ DATEIA.close();
+
+// cout << " Maximum: " << Maximum(interpol) << " Minimum: " << Minimum(interpol) << endl;
+
cout << " Maximum : " << L_infty(u) << endl;
cout << " Error : " << L_infty(interpol-u) << endl;
diff --git a/program/program.pro b/program/program.pro
index ca587ad..dd6241a 100644
--- a/program/program.pro
+++ b/program/program.pro
@@ -17,7 +17,7 @@ unix: MOC_DIR = MOC_FILES
# getting pathes:
-include(../config.pri)
+include(../configMy.pri)
INCLUDEPATH += $${UGBLOCKS_PATH}
diff --git a/program/source/extemp/variable.cc b/program/source/extemp/variable.cc
index 22395c2..1af0d3b 100644
--- a/program/source/extemp/variable.cc
+++ b/program/source/extemp/variable.cc
@@ -60,6 +60,24 @@ void Variable< std::complex<double> >::Update_back(int id_hex) const {
my_rank, comm);
}
+template <>
+void Variable<double>::AddUpFromHex(int id_hex) const {
+ AddUp_back_fromHex<double>(id_hex, blockgrid,
+ data_hexahedra, data_quadrangles,
+ data_edges, data_points,
+ my_rank, comm);
+}
+
+template <>
+void Variable< std::complex<double> >::AddUpFromHex(int id_hex) const {
+ AddUp_back_fromHex< std::complex<double> >(id_hex, blockgrid,
+ data_hexahedra, data_quadrangles,
+ data_edges, data_points,
+ my_rank, comm);
+}
+
+
+
template <>
template <>
void Variable<double>::Update<hexahedronEl>(int id_hex) const{
diff --git a/program/source/extemp/variable.h b/program/source/extemp/variable.h
index 1c58fe0..fea62d1 100644
--- a/program/source/extemp/variable.h
+++ b/program/source/extemp/variable.h
@@ -306,13 +306,21 @@ class Variable : public Expr<Variable< DTyp > >, public Object_based_on_ug {
inline DTyp Give_data ( params_in, D3vector &posVector ) const ;
void UpdateHexahedra(); ///> Updated alle Randdaten in allen Hexahedra
- void UpdateHexahedraBack();
+ void UpdateHexahedraBack();
/// todo update nur falls notwendig
template <elementTyp TYP_EL>
void Update ( int id ) const;
void Update_back(int id) const; ///> fuer interpolation slice 2D -> 3D
+
+ void AddUpFromHex(); ///> fuer Interpolation cell -> point
+ void AddUpFromHex(int id) const; ///> fuer Interpolation cell -> point
+ void SetAllZero(); ///> fuer Interpolation cell -> point
+ void AddFromOne(); ///> fuer Interpolation cell -> point
+ template <class Zell>
+ void AddFromZell(Zell& ZellVar); ///> fuer Interpolation cell -> point
+
inline DTyp Give_data_quadrangle_middle ( int id, int i, int j, int Nx ) const;
inline DTyp Give_data_quadrangle_interior ( int id, int i, int j, int Nx ) const;
@@ -430,6 +438,15 @@ void Variable<DTyp>::UpdateHexahedraBack() {
}
}
+template <class DTyp>
+void Variable<DTyp>::AddUpFromHex() {
+ for(int id_hex=0;id_hex<ug->Give_number_hexahedra();++id_hex) {
+ AddUpFromHex(id_hex);
+ }
+}
+
+
+
template <class DTyp>
int Variable<DTyp>::checkThreshold ( double threshold, int kz ) {
int id_hex = 1;
@@ -1494,11 +1511,123 @@ void Variable<DTyp>::operator= ( const Expr<A>& a ) {
}
+template <class DTyp>
+void Variable<DTyp>::SetAllZero() {
+ int Nx, Ny, Nz;
+
+ bool useOpenMP_here = UGBlocks::useOpenMP;
+
+ // hexahedra
+ for(int id = 0;id < ug->Give_number_hexahedra();++id ) {
+ if ( ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
+ Nx = blockgrid->Give_Nx_hexahedron ( id );
+ Ny = blockgrid->Give_Ny_hexahedron ( id );
+ Nz = blockgrid->Give_Nz_hexahedron ( id );
+
+ #pragma omp parallel for num_threads(UGBlocks::numThreadsToTake) if(useOpenMP_here)
+ for(int k = 0;k <= Nz;++k ) {
+ for(int j = 0;j <= Ny;++j )
+ for(int i = 0;i <= Nx;++i )
+ data_hexahedra[id][i+ ( Nx+1 ) * ( j+ ( Ny+1 ) *k ) ] = 0.0;
+ }
+ }
+ }
+
+ // quadrangles
+ for ( int id = 0;id < ug->Give_number_quadrangles();++id ) {
+
+ if ( ug->Give_quadrangle ( id )->my_object ( my_rank ) ) {
+ Nx = blockgrid->Give_Nx_quadrangle ( id );
+ Ny = blockgrid->Give_Ny_quadrangle ( id );
+
+ // #pragma omp parallel for private(i) num_threads(UGBlocks::numThreadsToTake) if(useOpenMP_here)
+ for(int j = 0;j <= Ny;++j )
+ for (int i = 0;i <= Nx;++i ) {
+ data_quadrangles[id][i+ ( Nx+1 ) *j] = 0.0;
+ }
+ }
+ }
+
+ // edges
+ for ( int id = 0;id < ug->Give_number_edges();++id ) {
+ if ( ug->Give_edge ( id )->my_object ( my_rank ) ) {
+ Nx = blockgrid->Give_Nx_edge ( id );
+
+ for(int i = 0;i <= Nx;++i )
+ data_edges[id][i] = 0.0;
+ }
+ }
+
+ // points
+ for(int id = 0;id < ug->Give_number_points();++id ) {
+ if(ug->Give_point ( id )->my_object ( my_rank ) ) {
+ data_points[id][0] = 0.0;
+ }
+ }
+}
+template <class DTyp>
+template <class Zell>
+void Variable<DTyp>::AddFromZell(Zell& ZellVar) {
+ int Nx, Ny, Nz;
+
+ // hexahedra
+ for(int id = 0;id < ug->Give_number_hexahedra();++id ) {
+ if(ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
+ Nx = blockgrid->Give_Nx_hexahedron ( id );
+ Ny = blockgrid->Give_Ny_hexahedron ( id );
+ Nz = blockgrid->Give_Nz_hexahedron ( id );
+
+ for(int k = 0;k < Nz;++k ) {
+ for(int j = 0;j < Ny;++j )
+ for(int i = 0;i < Nx;++i ) {
+ DTyp value = ZellVar.Give_cell_hexahedra(id,Ind_loc_matrix_hexahedra(i,j,k),i,j,k,Nx,Ny);
+ data_hexahedra[id][i + (Nx+1) * (j + (Ny+1)*k ) ] += value;
+ data_hexahedra[id][i+1 + (Nx+1) * (j + (Ny+1)*k ) ] += value;
+ data_hexahedra[id][i + (Nx+1) * (j+1 + (Ny+1)*k ) ] += value;
+ data_hexahedra[id][i+1 + (Nx+1) * (j+1 + (Ny+1)*k ) ] += value;
+ data_hexahedra[id][i + (Nx+1) * (j + (Ny+1)*(k+1) ) ] += value;
+ data_hexahedra[id][i+1 + (Nx+1) * (j + (Ny+1)*(k+1) ) ] += value;
+ data_hexahedra[id][i + (Nx+1) * (j+1 + (Ny+1)*(k+1) ) ] += value;
+ data_hexahedra[id][i+1 + (Nx+1) * (j+1 + (Ny+1)*(k+1) ) ] += value;
+ }
+ }
+ }
+ }
+}
+template <class DTyp>
+void Variable<DTyp>::AddFromOne() {
+ int Nx, Ny, Nz;
+
+ // hexahedra
+ for(int id = 0;id < ug->Give_number_hexahedra();++id ) {
+ if(ug->Give_hexahedron ( id )->my_object ( my_rank ) ) {
+ Nx = blockgrid->Give_Nx_hexahedron ( id );
+ Ny = blockgrid->Give_Ny_hexahedron ( id );
+ Nz = blockgrid->Give_Nz_hexahedron ( id );
+
+
+ for(int k = 0;k < Nz;++k ) {
+ for(int j = 0;j < Ny;++j )
+ for(int i = 0;i < Nx;++i ) {
+ DTyp value = 1.0;
+ data_hexahedra[id][i + (Nx+1) * (j + (Ny+1)*k ) ] += value;
+ data_hexahedra[id][i+1 + (Nx+1) * (j + (Ny+1)*k ) ] += value;
+ data_hexahedra[id][i + (Nx+1) * (j+1 + (Ny+1)*k ) ] += value;
+ data_hexahedra[id][i+1 + (Nx+1) * (j+1 + (Ny+1)*k ) ] += value;
+ data_hexahedra[id][i + (Nx+1) * (j + (Ny+1)*(k+1) ) ] += value;
+ data_hexahedra[id][i+1 + (Nx+1) * (j + (Ny+1)*(k+1) ) ] += value;
+ data_hexahedra[id][i + (Nx+1) * (j+1 + (Ny+1)*(k+1) ) ] += value;
+ data_hexahedra[id][i+1 + (Nx+1) * (j+1 + (Ny+1)*(k+1) ) ] += value;
+ }
+ }
+ }
+ }
+}
/////////////////////////////////////////////////////////////////////////
// 5. Nconst , ... funtions for copy of data
diff --git a/program/source/extemp/variable_cc.h b/program/source/extemp/variable_cc.h
index 87c68dc..7b30420 100644
--- a/program/source/extemp/variable_cc.h
+++ b/program/source/extemp/variable_cc.h
@@ -1083,6 +1083,167 @@ assert( rank_rec == 0);
+template <class DTyp>
+void AddUp_back_fromHex(int id_hex, Blockgrid* blockgrid, DTyp** data_hexahedra,
+ DTyp** data_quadrangles, DTyp** data_edges, DTyp** data_points,
+ int my_rank, MPI_Comm comm ) {
+ int i,j;
+ int id_quad, id_edge, id_point;
+ int Ni,Nj,Nconst;
+ int Nxedge;
+ int Nxquad,Nyquad;
+ int Nxhex,Nyhex,Nzhex;
+ int id_SW, id_SE, id_NW, idsp;
+ int id_W, id_E;
+ dir3D_sons c_SW ( ( dir3D_sons ) 0 ), c_SE ( ( dir3D_sons ) 0 ), c_NW ( ( dir3D_sons ) 0 );
+ dir3D_sons c_W ( ( dir3D_sons ) 0 ), c_E ( ( dir3D_sons ) 0 );
+ dir3D c_x, c_y;
+ Unstructured_grid* ug = blockgrid->Give_unstructured_grid();
+
+ int rank_send, rank_rec; // for parallel
+
+ Nxhex = blockgrid->Give_Nx_hexahedron ( id_hex );
+ Nyhex = blockgrid->Give_Ny_hexahedron ( id_hex );
+ Nzhex = blockgrid->Give_Nz_hexahedron ( id_hex );
+
+ rank_rec = ug->Give_hexahedron ( id_hex )->Give_rank(); // for parallel
+
+// keine MPI Paralleliserung!!!!!!!!!!!!!!
+assert( rank_rec == 0);
+
+ // quadrangles
+ for(int fc=0; fc<6; ++fc ) {
+ id_quad = ug->Give_hexahedron ( id_hex )->Give_id_quadrangle ( ( dir3D ) fc );
+ rank_send = ug->Give_quadrangle ( id_quad )->Give_rank(); // for parallel
+ if(rank_send==my_rank || rank_rec==my_rank ) { // for parallel
+ Nxquad = blockgrid->Give_Nx_quadrangle ( id_quad );
+ Nyquad = blockgrid->Give_Ny_quadrangle ( id_quad );
+
+ // a) coordinate points of quadrangle
+ id_SW = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SWdir2D );
+ id_SE = ug->Give_quadrangle ( id_quad )->Give_id_corner ( SEdir2D );
+ id_NW = ug->Give_quadrangle ( id_quad )->Give_id_corner ( NWdir2D );
+
+ // b) Find corresponding points in hexahedron
+ for(int c=0; c<8; ++c ) {
+ idsp = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
+ if ( idsp==id_SW ) c_SW = ( dir3D_sons ) c;
+ else if ( idsp==id_SE ) c_SE = ( dir3D_sons ) c;
+ else if ( idsp==id_NW ) c_NW = ( dir3D_sons ) c;
+ }
+
+ c_x = direction_from_to ( c_SW,c_SE );
+ c_y = direction_from_to ( c_SW,c_NW );
+
+ // x - direction
+ Nconst = Give_Nconst ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
+ Ni = Give_Ni ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
+
+ // y - direction
+ Nconst = Give_Nconst ( ( dir3D ) c_y, Nxhex, Nyhex, Nzhex ) + Nconst;
+ Nj = Give_Ni ( ( dir3D ) c_y, Nxhex, Nyhex, Nzhex );
+
+ // empty - direction
+ Nconst = Give_Nconst ( opposite3D ( ( dir3D ) fc ),
+ Nxhex, Nyhex, Nzhex ) + Nconst;
+ if(rank_send == rank_rec ) { // for parallel
+ for(j=1; j<Nyquad; ++j )
+ for(i=1;i<Nxquad; ++i ) {
+ // vorwaerts:
+ // data_hexahedra[id_hex][i*Ni+j*Nj+Nconst] =
+ // data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j];
+ data_quadrangles[id_quad][i+ ( Nxquad+1 ) *j] +=
+ data_hexahedra[id_hex][i*Ni+j*Nj+Nconst];
+ }
+ }
+ else {
+ assert(false);
+ }
+ }
+ }
+
+ // edges
+ for ( int ed=0; ed<12; ++ed ) {
+ id_edge = ug->Give_hexahedron ( id_hex )->Give_id_edge ( ( Edges_cell ) ed );
+ rank_send = ug->Give_edge ( id_edge )->Give_rank(); // for parallel
+ if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
+ Nxedge = blockgrid->Give_Nx_edge ( id_edge );
+
+ // a) coordinate points of edge
+ id_W = ug->Give_edge ( id_edge )->Give_id_corner_W();
+ id_E = ug->Give_edge ( id_edge )->Give_id_corner_E();
+
+ // b) Find corresponding points in hexahedron
+ for ( int c=0; c<8; ++c ) {
+ idsp = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
+ if ( idsp==id_W ) c_W = ( dir3D_sons ) c;
+ else if ( idsp==id_E ) c_E = ( dir3D_sons ) c;
+ }
+
+ c_x = direction_from_to ( c_W,c_E );
+
+ // x - direction
+ Nconst = Give_Nconst ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
+ Ni = Give_Ni ( ( dir3D ) c_x, Nxhex, Nyhex, Nzhex );
+
+ // A. empty - direction
+ Nconst = Give_Nconst ( opposite3D ( OrthoA ( ( Edges_cell ) ed ) ),
+ Nxhex, Nyhex, Nzhex ) + Nconst;
+
+ // B. empty - direction
+ Nconst = Give_Nconst ( opposite3D ( OrthoB ( ( Edges_cell ) ed ) ),
+ Nxhex, Nyhex, Nzhex ) + Nconst;
+
+ if(rank_send == rank_rec ) { // for parallel
+ for(i=1; i<Nxedge; ++i ) {
+ // vorwaerts
+ //data_hexahedra[id_hex][i*Ni+Nconst] = data_edges[id_edge][i];
+ data_edges[id_edge][i] += data_hexahedra[id_hex][i*Ni+Nconst];
+ }
+ }
+ else {
+ assert(false);
+ }
+ }
+ }
+
+ // points
+ for ( int c=0; c<8; ++c ) {
+ id_point = ug->Give_hexahedron ( id_hex )->Give_id_corner ( ( dir3D_sons ) c );
+ rank_send = ug->Give_point ( id_point )->Give_rank(); // for parallel
+ if ( rank_send==my_rank || rank_rec==my_rank ) { // for parallel
+
+ // x- direction
+ Nconst = Give_Nconst ( opposite3D ( xCoord ( ( dir3D_sons ) c ) ),
+ Nxhex, Nyhex, Nzhex );
+
+ // y- direction
+ Nconst = Give_Nconst ( opposite3D ( yCoord ( ( dir3D_sons ) c ) ),
+ Nxhex, Nyhex, Nzhex ) + Nconst;
+
+ // z- direction
+ Nconst = Give_Nconst ( opposite3D ( zCoord ( ( dir3D_sons ) c ) ),
+ Nxhex, Nyhex, Nzhex ) + Nconst;
+
+ if(rank_send == rank_rec ) // for parallel
+ // vorwaerts
+ //data_hexahedra[id_hex][Nconst] = data_points[id_point][0];
+ {
+ data_points[id_point][0] += data_hexahedra[id_hex][Nconst];
+ }
+ else {
+ assert(false);
+ }
+ }
+ }
+}
+
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////
+
diff --git a/program/source/interpol/interpolCellVar.cc b/program/source/interpol/interpolCellVar.cc
index 702cd42..fa703e5 100644
--- a/program/source/interpol/interpolCellVar.cc
+++ b/program/source/interpol/interpolCellVar.cc
@@ -26,6 +26,50 @@ using namespace ::_COLSAMM_;
#include "assert.h"
+FastInterpolate_from_cell_to_point::FastInterpolate_from_cell_to_point() {
+ weight = NULL;
+}
+
+void FastInterpolate_from_cell_to_point::deleteData() {
+ if(weight!=NULL) delete weight;
+ weight = NULL;
+}
+
+void FastInterpolate_from_cell_to_point::init(Blockgrid* blockgrid) {
+ deleteData();
+ weight = new Variable<double>(*blockgrid);
+ weight->SetAllZero();
+ weight->AddFromOne();
+ weight->AddUpFromHex();
+ weight->UpdateHexahedra();
+
+/*
+ std::ofstream DATEIA;
+
+ DATEIA.open("weight.vtk");
+ weight->Print_VTK(DATEIA);
+ DATEIA.close();
+
+ cout << " Maximum: " << Maximum(*weight) << " Minimum: " << Minimum(*weight) << endl;
+ */
+}
+
+void FastInterpolate_from_cell_to_point::doInterpolation(Cell_variable<double>* varCell, Variable<double>* varPoint) {
+ assert(weight!=NULL);
+ varPoint->SetAllZero();
+ varPoint->AddFromZell(*varCell);
+ varPoint->AddUpFromHex();
+
+ (*varPoint) = (*varPoint) / (*weight);
+ varPoint->UpdateHexahedra();
+
+
+
+}
+
+
+/////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////
template<>
@@ -68,3 +112,4 @@ void Interpolate_from_cell_to_point<Stencil_constant_z>::init(Blockgrid* blockgr
}
}
+
diff --git a/program/source/interpol/interpolCellVar.h b/program/source/interpol/interpolCellVar.h
index 7348969..9682ddd 100644
--- a/program/source/interpol/interpolCellVar.h
+++ b/program/source/interpol/interpolCellVar.h
@@ -48,6 +48,27 @@
#define INTERPOLCELLVAR_H
+/** \addtogroup InterpolationOperators **/
+/* @{ */
+/**
+ * Interpolation from Cell_variable to Variable <br>
+ * FAST VERSION <br>
+ * Sometimes this version is a little bit less accurate than Interpolate_from_cell_to_point in case of smooth functions <br>
+ * But in case of a non smooth function, this version is more appropriate
+**/
+class FastInterpolate_from_cell_to_point {
+ public:
+ FastInterpolate_from_cell_to_point();
+
+ void init(Blockgrid* blockgrid_);
+ void deleteData();
+ void doInterpolation(Cell_variable<double>* varCell, Variable<double>* varPoint);
+
+ ~FastInterpolate_from_cell_to_point() { deleteData(); }
+private:
+ Variable<double>* weight;
+};
+
@@ -56,8 +77,11 @@
/**
- * interpolation by Gauss-Seidel from Cell_variable to Variable
- * template might be Stencil_constant_z
+ * Interpolation by Gauss-Seidel from Cell_variable to Variable <br>
+ * template might be Stencil_constant_z <br>
+ * SLOW VERSION <br>
+ * Sometimes this version is more accurate than FastInterpolate_from_cell_to_point <br>
+ * An example showed an increase of factor 3 of accuracy in case of a smooth function
**/
template<class StenTyp>
class Interpolate_from_cell_to_point {
@@ -81,6 +105,8 @@ private:
int numberGSinterations;
};
+/* @} */
+
template<class StenTyp>
void Interpolate_from_cell_to_point<StenTyp>::doInterpolation(Cell_variable<double>* varCell, Variable<double>* varPoint) {
@@ -125,4 +151,5 @@ void Interpolate_from_cell_to_point<StenTyp>::deleteData() {
}
}
+
#endif // INTERPOLCELLVAR_H
--
GitLab