Skip to content

I
INMOST
Commit 7892a56d authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files
Options

Synchronize

parent 545619ce
Hide whitespace changes
Inline Side-by-side
Showing with 224 additions and 19 deletions
Examples/ADMFD/main.cpp
View file @ 7892a56d
......@@ -17,6 +17,7 @@ using namespace INMOST;
#endif
//shortcuts
typedef Storage::bulk bulk;
typedef Storage::real real;
typedef Storage::integer integer;
typedef Storage::enumerator enumerator;
......@@ -35,11 +36,24 @@ int main(int argc,char ** argv)
#endif
if( argc > 1 )
{
/*
std::cout << "Row: " << sizeof(Sparse::Row) << std::endl;
std::cout << "HessianRow: " << sizeof(Sparse::HessianRow) << std::endl;
std::cout << "Matrix: " << sizeof(Sparse::Matrix) << std::endl;
std::cout << "HessianMatrix: " << sizeof(Sparse::HessianMatrix) << std::endl;
std::cout << "variable: " << sizeof(multivar_expression) << std::endl;
Sparse::Row J;
Sparse::HessianRow H;
unknown a = unknown(0.5,1);
unknown b = unknown(0.1,2);
(a*a*b*(a+b)).GetHessian(1.0,J,1.0,H);
std::cout << "J: "; J.Print();
std::cout << "H: "; H.Print();
J.Clear();
(a*a*b*(a+b)).GetJacobian(1.0,J);
std::cout << "J: "; J.Print();
return 0;
*/
double ttt; // Variable used to measure timing
bool repartition = false; // Is it required to redistribute the mesh?
Mesh * m = new Mesh(); // Create an empty mesh
......@@ -99,7 +113,7 @@ int main(int argc,char ** argv)
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
Tag tag_DMP; // Indicates weather local W matrix satisfy DMP condition
if( m->GetProcessorsNumber() > 1 ) //skip for one processor job
{ // Exchange ghost cells
......@@ -159,10 +173,12 @@ int main(int argc,char ** argv)
}
m->ExchangeData(tag_P,CELL|FACE,0); //Synchronize initial solution with boundary unknowns
tag_W = m->CreateTag("nKGRAD",DATA_REAL,CELL,NONE);
tag_DMP = m->CreateTag("isDMP",DATA_BULK,CELL,NONE);
ttt = Timer();
//Assemble gradient matrix W on cells
int total = 0, dmp = 0;
#if defined(USE_OMP)
//#pragma omp parallel for
#pragma omp parallel for reduction(+:total) reduction(+:dmp)
#endif
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
{
......@@ -340,12 +356,27 @@ int main(int argc,char ** argv)
//std::cout << "W: " << std::endl;
//nKGRAD.Print();
bulk & isDMP = cell->Bulk(tag_DMP);
isDMP = 1;
for(int k = 0; k < NF; ++k)
{
real row_sum = 0;
if( nKGRAD(k,k) < 0.0 ) isDMP = 0;
for(int j = 0; j < NF; ++j)
row_sum += nKGRAD(k,j);
if( row_sum < 0.0 ) isDMP = 0;
for(int j = k+1; j < NF; ++j)
if( nKGRAD(k,j) > 0.0 )
isDMP = 0;
}
++total;
if( isDMP ) ++dmp;
real_array W = cell->RealArrayDV(tag_W); //access data structure for gradient matrix in mesh
W.resize(NF*NF); //resize the structure
std::copy(nKGRAD.data(),nKGRAD.data()+NF*NF,W.data()); //write down the gradient matrix
} //end of loop over cells
std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
std::cout << "Satisfy DMP: " << dmp << " out of " << total << std::endl;
if( m->HaveTag("FORCE") ) //Is there force on the mesh?
{
......@@ -399,8 +430,9 @@ int main(int argc,char ** argv)
{
R.Clear(); //clean up the residual
double tttt = Timer();
int total = 0, dmp = 0;
#if defined(USE_OMP)
#pragma omp parallel for
#pragma omp parallel for reduction(+:total) reduction(+:dmp)
#endif
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) //loop over cells
{
......@@ -408,11 +440,69 @@ int main(int argc,char ** argv)
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
int NF = (int)faces.size();
rMatrix nKGRAD(cell->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
bulk & isDMP = cell->Bulk(tag_DMP);
vMatrix pF(NF,1); //vector of pressure differences on faces
vMatrix FLUX(NF,1); //computed flux on faces
for(int k = 0; k < NF; ++k)
pF(k,0) = (P(faces[k]) - P(cell))*faces[k].Area();
FLUX = nKGRAD*pF; //fluxes on faces
//Change matrix by nonlinear correction for DMP
++total;
if( !isDMP )
//if( false )
{
const real var = 0;//0.25;
vMatrix vnKGRAD = nKGRAD;
vMatrix card(NF,1);
vMatrix dpF(NF,1);
//denominators of corrections
//this is sum of absolute values of differences
for(int k = 0; k < NF; ++k)
{
for(int j = 0; j < NF; ++j)
{
dpF(k,0) += soft_fabs(pF(j,0) - pF(k,0),1.0e-9);
card(k,0) += (soft_fabs(pF(j,0) - pF(k,0),1.0e-9)) / (soft_fabs(pF(j,0) - pF(k,0),1.0e-9)+1.0e-9);
}
}
//actual corrections
variable beta;
for(int k = 0; k < NF; ++k)
{
for(int j = k+1; j < NF; ++j) //if( nKGRAD(k,j) > 0.0 )
{
//beta = soft_fabs(FLUX(k,0),1.0e-9)/(dpF(k,0)) + soft_fabs(FLUX(j,0),1.0e-9)/(dpF(j,0));
beta = soft_max(soft_fabs(FLUX(k,0),1.0e-9)/card(k,0),soft_fabs(FLUX(j,0),1.0e-9)/card(j,0),1.0e-9)/soft_fabs(pF(k,0)-pF(j,0),1.0e-9);
beta = variation(beta,var);
vnKGRAD(k,j) -= (beta);
vnKGRAD(j,k) -= (beta);
vnKGRAD(k,k) += (beta);
vnKGRAD(j,j) += (beta);
}
}
FLUX = vnKGRAD*pF;
bool haveDMP = true;
for(int k = 0; k < NF; ++k)
{
real row_sum = 0;
if( vnKGRAD(k,k) < 0.0 ) haveDMP = false;
for(int j = 0; j < NF; ++j)
row_sum += get_value(vnKGRAD(k,j));
if( row_sum < 0.0 ) haveDMP = false;
for(int j = k+1; j < NF; ++j)
if( vnKGRAD(k,j) > 0.0 )
haveDMP = false;
}
if( !haveDMP )
{
//std::cout << "Failed to correct, matrix:" << std::endl;
//nKGRAD.Print();
}
else ++dmp;
}
else ++dmp;
if( cell.GetStatus() != Element::Ghost )
{
for(int k = 0; k < NF; ++k) //loop over faces of current cell
......@@ -434,6 +524,8 @@ int main(int argc,char ** argv)
}
} //end of loop over cells
std::cout << "Satisfy DMP: " << dmp << " out of " << total << std::endl;
if( tag_F.isValid() )
{
......@@ -452,8 +544,8 @@ int main(int argc,char ** argv)
R.Rescale();
R.GetJacobian().Save("jacobian.mtx",&Text);
R.GetResidual().Save("residual.mtx");
//R.GetJacobian().Save("jacobian.mtx",&Text);
//R.GetResidual().Save("residual.mtx");
std::cout << "Nonlinear residual: " << R.Norm() << "\t\t" << std::endl;
......@@ -464,7 +556,7 @@ int main(int argc,char ** argv)
S.SetMatrix(R.GetJacobian());
S.SetParameterReal("relative_tolerance", 1.0e-14);
S.SetParameterReal("absolute_tolerance", 1.0e-12);
S.SetParameterReal("drop_tolerance", 1.0e-3);
S.SetParameterReal("drop_tolerance", 1.0e-2);
S.SetParameterReal("reuse_tolerance", 1.0e-4);
//std::fill(Update.Begin(),Update.End(),0.0);
if( S.Solve(R.GetResidual(),Update) )
......
Source/Headers/inmost_dense.h
View file @ 7892a56d
......@@ -428,6 +428,13 @@ namespace INMOST
for(enumerator i = 0; i < n*m; ++i)
space[i] = other.space[i];
}
template<typename typeB>
Matrix(const Matrix<typeB> & other) : space(other.Cols()*other.Rows()), n(other.Rows()), m(other.Cols())
{
for(enumerator i = 0; i < n; ++i)
for(enumerator j = 0; j < m; ++j)
(*this)(i,j) = get_value(other(i,j));
}
~Matrix() {}
void Resize(enumerator nrows, enumerator mcols)
{
......@@ -445,6 +452,16 @@ namespace INMOST
m = other.m;
return *this;
}
template<typename typeB>
Matrix & operator =(Matrix<typeB> const & other)
{
if( n*m != other.n*other.m ) space.resize(other.n*other.m);
for(enumerator i = 0; i < other.n*other.m; ++i)
space[i] = get_value(other.space[i]);
n = other.n;
m = other.m;
return *this;
}
// i is in [0,n] - row index
// j is in [0,m] - column index
Var & operator()(enumerator i, enumerator j)
......
Source/Headers/inmost_expression.h
View file @ 7892a56d
......@@ -619,7 +619,8 @@ namespace INMOST
}
};
/// c/x = -c dx / (x*x)
/// (c/x)' = -c dx / (x*x)
/// (c/x)'' = 2 c dx dx / (x*x*x)
template<class A>
class reciprocal_expression : public shell_expression<reciprocal_expression<A> >
{
......@@ -644,7 +645,12 @@ namespace INMOST
}
__INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
{
arg.GetHessian(-multJ*value*reciprocial_val,J,2*multH*value*reciprocial_val*reciprocial_val,H);
Sparse::HessianRow ArgH;
double coefJ = -multJ*value*reciprocial_val;
double signJ = coefJ < 0 ? -1 : 1;
arg.GetHessian(signJ,J,-2*multH*value*reciprocial_val,ArgH);
Sparse::HessianRow::MergeJacobianHessian(2*value*reciprocial_val*reciprocial_val*signJ,J,J,1.0,ArgH,H);
for(INMOST_DATA_ENUM_TYPE k = 0; k < J.Size(); ++k) J.GetValue(k) *= coefJ*signJ;
}
};
......@@ -687,18 +693,21 @@ namespace INMOST
__INLINE INMOST_DATA_REAL_TYPE GetValue() const {return value;}
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
{
arg.GetJacobian( mult * dmult, r);
arg.GetJacobian( (value == 0 ? (mult < 0.0 ? -1 : 1) : 1) * mult * dmult, r);
}
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
{
arg.GetJacobian( mult * dmult, r);
arg.GetJacobian( (value == 0 ? (mult < 0.0 ? -1 : 1) : 1) * mult * dmult, r);
}
__INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ,Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
{
throw NotImplemented;
double a = (value == 0 ? (multJ < 0.0 ? -1 : 1) : 1);
double b = (value == 0 ? (multH < 0.0 ? -1 : 1) : 1);
arg.GetHessian( a * multJ * dmult,J, b*multH * dmult,H);
}
};
// ex
template<class A>
class exp_expression : public shell_expression<exp_expression<A> >
{
......@@ -722,7 +731,12 @@ namespace INMOST
}
__INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
{
arg.GetHessian(multJ*value, J, multH*value, H); //check
Sparse::HessianRow ArgH;
double coefJ = multJ*value;
double signJ = coefJ < 0.0 ? -1 : 1;
arg.GetHessian(signJ, J, multH*value, ArgH); //check
Sparse::HessianRow::MergeJacobianHessian(coefJ*signJ,J,J,1.0,ArgH,H);
for(INMOST_DATA_ENUM_TYPE k = 0; k < J.Size(); ++k) J.GetValue(k) *= coefJ*signJ;
}
};
......@@ -750,7 +764,12 @@ namespace INMOST
}
__INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
{
arg.GetJacobian(multJ*dmult,J,-multH*dmult*dmult,H);
Sparse::HessianRow ArgH;
double coefJ = multJ*dmult;
double signJ = coefJ < 0.0 ? -1 : 1;
arg.GetHessian(signJ, J, 2*multH*dmult, ArgH); //check
Sparse::HessianRow::MergeJacobianHessian(-coefJ*signJ*dmult,J,J,1.0,ArgH,H);
for(INMOST_DATA_ENUM_TYPE k = 0; k < J.Size(); ++k) J.GetValue(k) *= coefJ*signJ;
}
};
......@@ -1002,7 +1021,7 @@ namespace INMOST
Sparse::Row::MergeSortedRows(right.GetValue(),JL,left.GetValue(),JR,J);
//preallocate H to HL.Size+HR.Size+JL.Size*JR.Size
//merge sorted
Sparse::HessianRow::MergeJacobianHessian(1.0,JL,JR,right.GetValue(),HL,left.GetValue(),HR,H);
Sparse::HessianRow::MergeJacobianHessian(2.0,JL,JR,right.GetValue(),HL,left.GetValue(),HR,H);
}
};
......@@ -1046,7 +1065,7 @@ namespace INMOST
Sparse::Row::MergeSortedRows(reciprocal_rval,JL,-value * reciprocal_rval,JR,J);
//preallocate H to HL.Size+HR.Size+JL.Size*JR.Size
//merge sorted
Sparse::HessianRow::MergeJacobianHessian(1.0,JL,JR,reciprocal_rval,HL,2*left.GetValue()*multH,HR,H);
Sparse::HessianRow::MergeJacobianHessian(2.0,JL,JR,reciprocal_rval,HL,2*left.GetValue()*multH,HR,H);
}
};
......
Source/Headers/inmost_sparse.h
View file @ 7892a56d
......@@ -296,6 +296,7 @@ namespace INMOST
/// output = alpha * left + beta *right
static void MergeSortedRows(INMOST_DATA_REAL_TYPE alpha, const HessianRow & left, INMOST_DATA_REAL_TYPE beta, const HessianRow & right, HessianRow & output);
static void MergeJacobianHessian(INMOST_DATA_REAL_TYPE a, const Row & JL, const Row & JR, INMOST_DATA_REAL_TYPE b, const HessianRow & HL, INMOST_DATA_REAL_TYPE c, const HessianRow & HR, HessianRow & output);
static void MergeJacobianHessian(INMOST_DATA_REAL_TYPE a, const Row & JL, const Row & JR, INMOST_DATA_REAL_TYPE b, const HessianRow & H, HessianRow & output);
};
/// This class can be used for shared access to matrix with OpenMP.
......
Source/Solver/sparse.cpp
View file @ 7892a56d
......@@ -1092,7 +1092,7 @@ namespace INMOST
if( j < HR.Size() )
candidate[1] = make_entry(HR.GetIndex(j),c*HR.GetValue(j));
if( k < JL.Size() && l < JR.Size() )
candidate[2] = make_entry(make_index(JL.GetIndex(kk),JR.GetIndex(ll)),2*a*JL.GetValue(kk)*JR.GetValue(ll));
candidate[2] = make_entry(make_index(JL.GetIndex(kk),JR.GetIndex(ll)),a*JL.GetValue(kk)*JR.GetValue(ll));
do
{
//pick smallest
......@@ -1152,7 +1152,7 @@ namespace INMOST
}
}
if( kk < JL.Size() && ll < JR.Size() )
candidate[2] = make_entry(make_index(JL.GetIndex(kk),JR.GetIndex(ll)),2*a*JL.GetValue(kk)*JR.GetValue(ll));
candidate[2] = make_entry(make_index(JL.GetIndex(kk),JR.GetIndex(ll)),a*JL.GetValue(kk)*JR.GetValue(ll));
else
candidate[2] = stub_entry;
}
......@@ -1160,5 +1160,81 @@ namespace INMOST
while(true);
output.Resize(q);
}
void HessianRow::MergeJacobianHessian(INMOST_DATA_REAL_TYPE a, const Row & JL, const Row & JR, INMOST_DATA_REAL_TYPE b, const HessianRow & H, HessianRow & output)
{
// merge three sorted arrays at once
// one of the array is synthesized from JL and JR on the go
static const entry stub_entry = make_entry(make_index(ENUMUNDEF,ENUMUNDEF),0.0);
assert(JL.isSorted());
assert(JR.isSorted());
assert(H.isSorted());
output.Resize(H.Size()+JL.Size()*JR.Size());
INMOST_DATA_ENUM_TYPE i = 0, j = 0, k = 0, l = 0, q = 0, kk = 0, ll = 0, r;
entry candidate[2] = {stub_entry,stub_entry};
if( i < H.Size() )
candidate[0] = make_entry(H.GetIndex(i),b*H.GetValue(i));
if( k < JL.Size() && l < JR.Size() )
candidate[1] = make_entry(make_index(JL.GetIndex(kk),JR.GetIndex(ll)),a*JL.GetValue(kk)*JR.GetValue(ll));
do
{
//pick smallest
r = 0;
if( candidate[1].first < candidate[r].first ) r = 1;
//all candidates are stub - exit
if( candidate[r].first == stub_entry.first ) break;
//record selected entry
if( q > 0 && (output.GetIndex(q-1) == candidate[r].first) )
output.GetValue(q-1) += candidate[r].second;
else
{
output.GetIndex(q) = candidate[r].first;
output.GetValue(q) = candidate[r].second;
++q;
}
if( r == 0 ) //update left hessian index
{
if( ++i < H.Size() )
candidate[0] = make_entry(H.GetIndex(i),b*H.GetValue(i));
else candidate[0] = stub_entry;
}
else //update jacobians indexes
{
if( JR.GetIndex(l) < JL.GetIndex(k) )
{
if( ++kk == JL.Size() )
{
++l;
kk = k;
ll = l;
}
}
else if( JL.GetIndex(k) < JR.GetIndex(l) )
{
if( ++ll == JR.Size() )
{
++k;
kk = k;
ll = l;
}
}
else //values are equal
{
if( ++ll == JR.Size() )
{
++k;
kk = k;
ll = l;
}
}
if( kk < JL.Size() && ll < JR.Size() )
candidate[1] = make_entry(make_index(JL.GetIndex(kk),JR.GetIndex(ll)),a*JL.GetValue(kk)*JR.GetValue(ll));
else
candidate[1] = stub_entry;
}
}
while(true);
output.Resize(q);
}
}
}
\ No newline at end of file
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment