Commit 2014d0bc authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Merge branch 'master' of https://github.com/INMOST-DEV/INMOST

parents fa343ffe fb25619f
......@@ -22,3 +22,6 @@ Source/Solver/solver_k3biilu2/k3d_block.hxx
.DS_Store
/.vs
/out/install/x64-Release
/CMakeSettings.json
......@@ -26,6 +26,19 @@ typedef Storage::real_array real_array;
typedef Storage::var_array var_array;
bool print_niter = false; //save file on nonlinear iterations
bool output_matrix = true;
bool norne_wells = true;
//data for wells
Cell cc[3] = {InvalidCell(),InvalidCell(),InvalidCell()};
real WI[3] = {50000,50000,50000};
real pbhp[3] = {265,105,110};
real ccnt[3][3] =
{
{4.567151e+05, 7.321079e+06, 2.767665e+03},
{4.609346e+05, 7.323503e+06, 2.597767e+03},
{4.595400e+05, 7.326078e+06, 2.803586e+03}
};
//#define OPTIMIZATION
......@@ -101,8 +114,8 @@ int main(int argc,char ** argv)
Tag tag_P; // Pressure
Tag tag_K; // Diffusion tensor
Tag tag_F; // Forcing term
Tag tag_BC; // Boundary conditions
Tag tag_W; // Gradient matrix acting on harmonic points on faces and returning gradient on faces
TagRealArray tag_BC; // Boundary conditions
......@@ -139,6 +152,18 @@ int main(int argc,char ** argv)
m->ExchangeData(tag_K,CELL,0); //Exchange diffusion tensor
}
if( norne_wells )
{
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it) if( it->GetStatus() != Element::Ghost )
{
for(int q = 0; q < 3; ++q) if( it->Inside(ccnt[q]) )
{
cc[q] = it->self();
std::cout << "proc " << m->GetProcessorRank() << " found c" << q << " " << cc[q].LocalID() << std::endl;
}
}
}
if( m->HaveTag("PRESSURE") ) //Is there a pressure on the mesh?
tag_P = m->GetTag("PRESSURE"); //Get the pressure
......@@ -219,9 +244,9 @@ int main(int argc,char ** argv)
//~ W+=(rMatrix::Unit(NF) - R*(R.Transpose()*R).CholeskyInvert()*R.Transpose())*
//~ (2.0/(static_cast<real>(NF)*volume)*(NK*K.CholeskyInvert()*NK.Transpose()).Trace());
double ca = 0, ba = ca;
W = (NK+ca*L*R)*((NK+ca*L*R).Transpose()*R).Invert()*(NK+ca*L*R).Transpose();
W+= L - (1+ba)*(L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
W = (NK)*((NK).Transpose()*R).Invert()*(NK).Transpose();
W+= L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
//~ W = Areas*W*Areas;
//access data structure for gradient matrix in mesh
......@@ -232,7 +257,7 @@ int main(int argc,char ** argv)
std::copy(W.data(),W.data()+NF*NF,store_W.data());
tag_WG[cell].resize(3*NF);
tag_WG(cell,3,NF) = (NK.Transpose()*R+ca*L*R).PseudoInvert(1.0e-12)*(NK+ca*L*R).Transpose();
tag_WG(cell,3,NF) = (NK.Transpose()*R).PseudoInvert(1.0e-12)*(NK).Transpose();
} //end of loop over cells
}
std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
......@@ -254,7 +279,32 @@ int main(int argc,char ** argv)
{ //Main loop for problem solution
Automatizator aut; // declare class to help manage unknowns
Automatizator::MakeCurrent(&aut);
dynamic_variable P(aut,aut.RegisterTag(tag_P,CELL|FACE)); //register pressure as primary unknown
MarkerType unk = m->CreateMarker();
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
m->CellByLocalID(q).SetMarker(unk);
for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
{
Face face = m->FaceByLocalID(q);
face.SetMarker(unk);
if( face.Boundary() )
{
double alpha = 0, beta = 1, gamma = 0;
if( tag_BC.isValid() && face.HaveData(tag_BC) )
{
alpha = tag_BC[face][0];
beta = tag_BC[face][1];
gamma = tag_BC[face][2];
}
if( alpha != 0 && beta == 0 )
{
face.RemMarker(unk);
face.Real(tag_P) = gamma/alpha;
}
}
}
dynamic_variable P(aut,aut.RegisterTag(tag_P,CELL|FACE,unk)); //register pressure as primary unknown
aut.EnumerateEntries(); //enumerate all primary variables
std::cout << "Enumeration done, size " << aut.GetLastIndex() - aut.GetFirstIndex() << std::endl;
......@@ -272,10 +322,21 @@ int main(int argc,char ** argv)
for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
{
Face face = m->FaceByLocalID(q);
if( face.GetStatus() != Element::Ghost )
if( face.GetStatus() != Element::Ghost && face.GetMarker(unk) )
{
if( tag_BC.isValid() && face.HaveData(tag_BC) )
Text.SetAnnotation(P.Index(face),"Pressure guided by boundary condition");
if( face.Boundary() )
{
double alpha = 0, beta = 1, gamma = 0;
if( tag_BC.isValid() && face.HaveData(tag_BC) )
{
alpha = tag_BC[face][0];
beta = tag_BC[face][1];
gamma = tag_BC[face][2];
}
std::stringstream str;
str << "Pressure guided by boundary condition " << alpha << " " << beta << " " << gamma;
Text.SetAnnotation(P.Index(face),str.str());
}
else
Text.SetAnnotation(P.Index(face),"Interface pressure");
}
......@@ -311,7 +372,7 @@ int main(int argc,char ** argv)
FLUX.Resize(NF,1);
for(int k = 0; k < NF; ++k)
pF(k,0) = (P(faces[k]) - P(cell));
pF(k,0) = (P[faces[k]] - P[cell]);
FLUX = W*pF; //fluxes on faces
tag_PG(cell,3,1) = tag_WG(cell,3,NF)*pF;
if( cell.GetStatus() != Element::Ghost )
......@@ -322,6 +383,7 @@ int main(int argc,char ** argv)
for(int k = 0; k < NF; ++k) //loop over faces of current cell
{
if( faces[k].GetStatus() == Element::Ghost ) continue;
if( !faces[k].GetMarker(unk) ) continue;
int index = P.Index(faces[k]);
Locks.Lock(index);
......@@ -363,6 +425,12 @@ int main(int argc,char ** argv)
}
}
if( norne_wells )
{
for(int q = 0; q < 3; ++q) if( cc[q].isValid() )
R[P.Index(cc[q])] += WI[q]*(pbhp[q] - P[cc[q]])*cc[q].Volume();
}
//R[P.Index(m->BeginCell()->self())] = P[m->BeginCell()->self()];
}
std::cout << "assembled in " << Timer() - tttt << "\t\t\t" << std::endl;
......@@ -371,20 +439,53 @@ int main(int argc,char ** argv)
if( R.Norm() < 1.0e-4 ) break;
//~ R.GetJacobian().Save("A.mtx");
//~ R.GetResidual().Save("b.mtx");
//Solver S(Solver::INNER_ILU2);
//~ Solver S(Solver::INNER_MPTILUC);
Solver S(Solver::K3BIILU2);
Solver S(Solver::INNER_MPTILUC);
//~ Solver S(Solver::INNER_MLMPTILUC);
//~ Solver S(Solver::K3BIILU2);
//Solver S("superlu");
S.SetParameter("verbosity","1");
S.SetParameter("verbosity","3");
S.SetParameter("rescale_iterations", "10");
S.SetParameter("relative_tolerance", "1.0e-14");
S.SetParameter("absolute_tolerance", "1.0e-12");
S.SetParameter("drop_tolerance", "1.0e-3");
//~ S.SetParameter("drop_tolerance", "5.0e-2");
//~ S.SetParameter("reuse_tolerance", "2.5e-3");
S.SetParameter("drop_tolerance", "1.0e-2");
S.SetParameter("reuse_tolerance", "1.0e-4");
S.SetParameter("pivot_condition", "20");
double tset = Timer(), titr;
S.SetMatrix(R.GetJacobian());
tset = Timer() - tset;
if( output_matrix )
{
std::cout << "write A.mtx" << std::endl;
R.GetJacobian().Save("A.mtx");//,&Text);
std::cout << "write b.mtx" << std::endl;
R.GetResidual().Save("b.mtx");
std::cout << "write done, solve" << std::endl;
}
titr = Timer();
if( S.Solve(R.GetResidual(),Update) )
bool success = S.Solve(R.GetResidual(),Update);
titr = Timer() - titr;
if( success )
{
std::cout << "Solved system in " << S.Iterations() << " iterations, time for setup " << tset << "s iterations " << titr << "s" << std::endl;
if( output_matrix )
{
std::cout << "write x.mtx" << std::endl;
Update.Save("x.mtx");
//std::cout << "exit" << std::endl;
//exit(-1);
}
#if defined(USE_OMP)
#pragma omp parallel for
#endif
......@@ -392,6 +493,7 @@ int main(int argc,char ** argv)
{
Cell cell = m->CellByLocalID(q);
if( cell->GetStatus() == Element::Ghost ) continue;
if( !cell->GetMarker(unk) ) continue;
cell->Real(tag_P) -= Update[P.Index(cell)];
}
#if defined(USE_OMP)
......@@ -401,6 +503,7 @@ int main(int argc,char ** argv)
{
Face face = m->FaceByLocalID(q);
if( face->GetStatus() == Element::Ghost ) continue;
if( !face->GetMarker(unk) ) continue;
face->Real(tag_P) -= Update[P.Index(face)];
}
m->ExchangeData(tag_P, CELL|FACE, 0);
......
......@@ -8,6 +8,10 @@ bool balance_mesh_refine = false;
bool balance_mesh_coarse = false;
std::string file_format = ".pmf";
#ifndef M_PI
#define M_PI 3.1415926535897932
#endif
int main(int argc, char ** argv)
{
Mesh::Initialize(&argc,&argv);
......
......@@ -1217,6 +1217,107 @@ double display_elem_info(Element e, double top, double left, double interval)
return top-2*interval;
}
double display_elem_info_cout(Element e)
{
printf("%s %d\n",ElementTypeName(e->GetElementType()),e->LocalID());
if( e->GetElementType() == ESET ) printf(" size %d\n",e->getAsSet().Size());
for(Mesh::iteratorTag t = mesh->BeginTag(); t != mesh->EndTag(); ++t) if( t->isDefined(e->GetElementType()) )
{
if( e->HaveData(*t) )
{
char str[131072];
char temp[131072];
str[0] = '\0';
int dsize;
switch(t->GetDataType())
{
case DATA_INTEGER:
{
Storage::integer_array arr = e->IntegerArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{
sprintf(temp,"%s %d",str,arr[k]);
strcpy(str,temp);
}
break;
}
case DATA_REAL:
{
Storage::real_array arr = e->RealArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{
sprintf(temp,"%s %lf",str,arr[k]);
strcpy(str,temp);
}
break;
}
case DATA_BULK:
{
Storage::bulk_array arr = e->BulkArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{
sprintf(temp,"%s %d",str,arr[k]);
strcpy(str,temp);
}
break;
}
case DATA_REFERENCE:
{
Storage::reference_array arr = e->ReferenceArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{
if(arr.at(k) == InvalidHandle()) sprintf(temp,"%s NULL",str);
else sprintf(temp,"%s %s:%d",str,ElementTypeName(arr[k]->GetElementType()),arr[k]->LocalID());
strcpy(str,temp);
}
break;
}
case DATA_REMOTE_REFERENCE:
{
Storage::remote_reference_array arr = e->RemoteReferenceArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{
if(arr.at(k).first == NULL || arr.at(k).second == InvalidHandle()) sprintf(temp,"%s NULL",str);
else sprintf(temp,"%s %p:%s:%d",str,arr[k]->GetMeshLink(),ElementTypeName(arr[k]->GetElementType()),arr[k]->LocalID());
strcpy(str,temp);
}
break;
}
#if defined(USE_AUTODIFF)
case DATA_VARIABLE:
{
Storage::var_array arr = e->VariableArray(*t);
dsize = arr.size();
for(INMOST_DATA_ENUM_TYPE k = 0; k < arr.size(); k++)
{
std::stringstream stream;
stream << arr[k].GetValue() << " {[" << arr[k].GetRow().Size() << "] ";
for(INMOST_DATA_ENUM_TYPE q = 0; q < arr[k].GetRow().Size(); ++q)
{
stream << "(" << arr[k].GetRow().GetValue(q) << "," << arr[k].GetRow().GetIndex(q) << ") ";
}
stream << "}" << std::endl;
sprintf(temp,"%s %s",str,stream.str().c_str());
strcpy(str,temp);
}
break;
}
#endif
}
sprintf(temp,"%s[%d] %s\n %s",t->GetTagName().c_str(),dsize,DataTypeName(t->GetDataType()),str);
strcpy(str,temp);
printf("%s\n",str);
}
}
return 0;
}
void ProcessCommonInput(char inpstr[8192], int inptype)
{
......@@ -1420,6 +1521,7 @@ void ProcessCommonInput(char inpstr[8192], int inptype)
if (disp_e.isValid())
{
display_elem_info_cout(disp_e);
if (disp_e.GetElementType() != ESET)
{
disp_e->Centroid(shift);
......
......@@ -19,6 +19,35 @@ int zoom = 200;
int block_size = 0;
int draw_color = 0;
struct block_t
{
int rows, rowe;
int cols, cole;
block_t(int rows, int rowe, int cols, int cole)
: rows(rows), rowe(rowe), cols(cols), cole(cole){}
void drawgl(int mat_size)
{
glBegin(GL_LINES);
//top line
glVertex2i(cols,mat_size - rows);
glVertex2i(cole,mat_size - rows);
//bottom line
glVertex2i(cols,mat_size - rowe);
glVertex2i(cole,mat_size - rowe);
//left line
glVertex2i(cols,mat_size - rows);
glVertex2i(cols,mat_size - rowe);
//right line
//left line
glVertex2i(cole,mat_size - rows);
glVertex2i(cole,mat_size - rowe);
glEnd();
}
};
std::vector< block_t > blocks;
double amin = 1e20,amax = -1e20;
void printtext(const char * fmt, ...)
......@@ -550,7 +579,14 @@ void draw()
}
glEnd();
}
if( !blocks.empty() )
{
glColor3f(0.5,0,1);
//rows
for(int k = 0; k < blocks.size(); ++k)
blocks[k].drawgl(m->Size());
}
if (CommonInput != NULL)
{
glDisable(GL_DEPTH_TEST);
......@@ -576,7 +612,7 @@ int main(int argc, char ** argv)
{
if( argc < 2 )
{
std::cout << "usage: " << argv[0] << " matrix.mtx [text]" << std::endl;
std::cout << "usage: " << argv[0] << " matrix.mtx [(text [width]|blocks.txt)]" << std::endl;
return -1;
}
m = new Sparse::Matrix();
......@@ -584,15 +620,15 @@ int main(int argc, char ** argv)
if( argc > 2 )
{
int recw = (argc > 3 ? (atoi(argv[3])+4) : 10);
if( recw <= 4 )
{
std::cout << "Bad record width " << argv[3] << " on input, setting 10" << std::endl;
recw = 10;
}
std::string line(recw,'-');
if( std::string(argv[2]) == "text" )
{
int recw = (argc > 4 ? (atoi(argv[3])+4) : 10);
if( recw <= 4 )
{
std::cout << "Bad record width " << argv[3] << " on input, setting 10" << std::endl;
recw = 10;
}
std::string line(recw,'-');
std::cout << std::fixed << std::setprecision(recw-4);
std::cout << " ";
for(int k = 0; k < (int)m->Size(); ++k)
......@@ -633,9 +669,21 @@ int main(int argc, char ** argv)
for(int k = 0; k < (int)m->Size(); ++k)
std::cout << line;
std::cout << "*" << std::endl;
delete m;
return 0;
}
else
{
std::ifstream file(argv[2]);
INMOST_DATA_ENUM_TYPE rows, rowe, cols, cole;
while( !file.eof() )
{
file >> rows >> rowe >> cols >> cole;
blocks.push_back(block_t(rows,rowe,cols,cole));
}
file.close();
}
delete m;
return 0;
}
//ord = new Reorder_ARMS(m,0,m->Size());
......
......@@ -10,9 +10,17 @@ class CubeTransform
Mesh & mesh;
double xyz[8][3];
public:
CubeTransform(Mesh & m) : mesh(m),
xyz({{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}})
{}
CubeTransform(Mesh & m) : mesh(m)
{
xyz[0][0] = 0; xyz[0][1] = 0; xyz[0][2] = 0;
xyz[1][0] = 1; xyz[1][1] = 0; xyz[1][2] = 0;
xyz[2][0] = 0; xyz[2][1] = 1; xyz[2][2] = 0;
xyz[3][0] = 1; xyz[3][1] = 1; xyz[3][2] = 0;
xyz[4][0] = 0; xyz[4][1] = 0; xyz[4][2] = 1;
xyz[5][0] = 1; xyz[5][1] = 0; xyz[5][2] = 1;
xyz[6][0] = 0; xyz[6][1] = 1; xyz[6][2] = 1;
xyz[7][0] = 1; xyz[7][1] = 1; xyz[7][2] = 1;
}
CubeTransform(const CubeTransform & b) : mesh(b.mesh) {}
CubeTransform & operator = (CubeTransform const & b) {mesh = b.mesh; return *this;}
~CubeTransform() {}
......
......@@ -3,6 +3,9 @@
#include "fracture.h"
#include "fix_faults.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
class SliceFault : public Slice
{
......
......@@ -249,7 +249,7 @@ namespace INMOST
}
}
}
std::set<INMOST_DATA_ENUM_TYPE> Pre, Post; //Nonlocal indices
//~ std::set<INMOST_DATA_ENUM_TYPE> Pre, Post; //Nonlocal indices
#if defined(USE_MPI)
int size;
MPI_Comm_size(MPI_COMM_WORLD,&size);
......@@ -283,6 +283,7 @@ namespace INMOST
for(std::map<Mesh *,std::vector<Tag> >::iterator it = exch_tags.begin(); it != exch_tags.end(); ++it)
it->first->ExchangeData(it->second, exch_mask,0);
}
#if 0
//compute out-of-bounds indices
for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] )
{
......@@ -309,6 +310,7 @@ namespace INMOST
}
} //etype
} //it
#endif
// after cycle
}
#endif
......@@ -324,7 +326,8 @@ namespace INMOST
{
merger.resize(MAX_THREADS);
}
merger[OMP_THREAD].Resize(first_num,last_num,std::vector<INMOST_DATA_ENUM_TYPE>(Pre.begin(),Pre.end()),std::vector<INMOST_DATA_ENUM_TYPE>(Post.begin(),Post.end()),false);
//~ merger[OMP_THREAD].Resize(first_num,last_num,std::vector<INMOST_DATA_ENUM_TYPE>(Pre.begin(),Pre.end()),std::vector<INMOST_DATA_ENUM_TYPE>(Post.begin(),Post.end()),false);
merger[OMP_THREAD].Resize(first_num,last_num,false);
}
}
......
......@@ -582,12 +582,14 @@ namespace INMOST
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type, pool_array_t<typename Promote<Var,typeB>::type> >
operator*(typeB coef) const;
#if defined(USE_AUTODIFF)
/// Multiply the matrix by a coefficient.
/// @param coef Coefficient.
/// @return Matrix multiplied by the coefficient.
template<class A>
Matrix<typename Promote<Var,variable>::type, pool_array_t<typename Promote<Var,variable>::type> >
operator*(shell_expression<A> const & coef) const {return operator*(variable(coef));}
#endif //USE_AUTODIFF
/// Multiply the matrix by the coefficient of the same type and store the result.
/// @param coef Coefficient.
/// @return Reference to the current matrix.
......@@ -599,12 +601,14 @@ namespace INMOST
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type, pool_array_t<typename Promote<Var,typeB>::type> >
operator/(typeB coef) const;
#if defined(USE_AUTODIFF)
/// Divide the matrix by a coefficient of a different type.
/// @param coef Coefficient.
/// @return Matrix divided by the coefficient.
template<class A>
Matrix<typename Promote<Var,variable>::type, pool_array_t<typename Promote<Var,variable>::type> >
operator/(shell_expression<A> const & coef) const {return operator/(variable(coef));}
#endif //USE_AUTODIFF
/// Divide the matrix by the coefficient of the same type and store the result.
/// @param coef Coefficient.
/// @return Reference to the current matrix.
......@@ -2342,15 +2346,16 @@ namespace INMOST
assign((*this)(i,j),(*this)(i,j)+other(i,j));
return *this;
}
#if defined(USE_AUTODIFF)
template<>
template<>
__INLINE Matrix<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type> >
__INLINE Matrix<Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<Promote<INMOST_DATA_REAL_TYPE,variable>::type> >
AbstractMatrix<INMOST_DATA_REAL_TYPE>::operator*<variable>(const AbstractMatrix<variable> & other) const
{
//~ std::cout << __FUNCTION__ << std::endl;
assert(Cols() == other.Rows());
Matrix<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<typename Promote<INMOST_DATA_REAL_TYPE,variable>::type> > ret(Rows(),other.Cols()); //check RVO
Matrix<Promote<INMOST_DATA_REAL_TYPE,variable>::type, pool_array_t<Promote<INMOST_DATA_REAL_TYPE,variable>::type> > ret(Rows(),other.Cols()); //check RVO
for(enumerator i = 0; i < Rows(); ++i) //loop rows
{
for(enumerator j = 0; j < other.Cols(); ++j) //loop columns
......@@ -2384,12 +2389,12 @@ namespace INMOST
template<>
template<>
__INLINE Matrix<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type, pool_array_t<typename Promote<variable,INMOST_DATA_REAL_TYPE>::type> >