Skip to content

GitLab

Commit 7d30608c authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Synchronize 

Adapted ADFVdiscr for Residual class.

Matrix<Var>::DotProduct will promote type depending on the argument.
parent 64f740bd
Hide whitespace changes
Inline Side-by-side
Examples/ADFVDiscr/main.cpp
View file @ 7d30608c
......@@ -133,8 +133,8 @@ int main(int argc,char ** argv)
ttt = Timer();
Solver S(Solver::INNER_ILU2); // Specify the linear solver to ASM+ILU2+BiCGStab one
S.SetParameterReal("absolute_tolerance",1e-8);
Sparse::Matrix A; // Declare the matrix of the linear system to be solved
Sparse::Vector x,b; // Declare the solution and the right-hand side vectors
Residual R; // Residual vector
Sparse::Vector Update; // Declare the solution and the right-hand side vectors
Mesh::GeomParam table;
......@@ -154,10 +154,9 @@ int main(int argc,char ** argv)
aut.EnumerateDynamicTags();
// Set the indeces intervals for the matrix and vectors
A.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
x.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
b.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
//~ std::cout << m->GetProcessorRank() << " A,x,b interval " << idmin << ":" << idmax << " size " << idmax-idmin << std::endl;
R.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
Update.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
//~ std::cout << m->GetProcessorRank() << " A,x,b interval " << idmin << ":" << idmax << " size " << idmax-idmin << std::endl;
dynamic_variable Phi(aut,iphi);
// Solve \nabla \cdot \nabla phi = f equation
//for( Mesh::iteratorFace face = m->BeginFace(); face != m->EndFace(); ++face )
......@@ -195,14 +194,9 @@ int main(int argc,char ** argv)
bnd_pnt[2] = r1_cnt[2] + dist * f_nrm[2];
T = r1->Real(tensor_K) * f_area / dist;
//flux = T * (func(bnd_pnt,0) - variable(aut,r1,iphi));
flux = T * (func(bnd_pnt,0) - Phi(r1));
A[i1].Lock();
r.PushRow(1.0,A[i1]);
r.AddRow(-1.0,flux.GetRow());
r.RetriveRow(A[i1]);
r.Clear();
b[i1] -= flux.GetValue();
A[i1].Unlock();
R[i1].Lock();
R[i1] -= T * (func(bnd_pnt,0) - Phi(r1));
R[i1].Unlock();
}
else
{
......@@ -216,26 +210,17 @@ int main(int argc,char ** argv)
flux = T * (Phi(r2) - Phi(r1));
if( s1 != Element::Ghost )
{
A[i1].Lock();
r.PushRow(1.0,A[i1]);
r.AddRow(-1.0,flux.GetRow());
r.RetriveRow(A[i1]);
r.Clear();
b[i1] -= flux.GetValue();
A[i1].Unlock();
R[i1].Lock();
R[i1] -= flux;
R[i1].Unlock();
}
if( s2 != Element::Ghost )
{
A[i2].Lock();
r.PushRow(1.0,A[i2]);
r.AddRow(1.0,flux.GetRow());
r.RetriveRow(A[i2]);
r.Clear();
b[i2] += flux.GetValue();
A[i2].Unlock();
R[i2].Lock();
R[i2] += flux;
R[i2].Unlock();
}
}
}
}
#if defined(USE_OMP)
......@@ -245,7 +230,7 @@ int main(int argc,char ** argv)
{
Cell cell = Cell(m,ComposeCellHandle(icell));
if( cell->GetStatus() != Element::Ghost )
b[cell->Integer(id)] += cell->Mean(func_rhs, cell->Real(tensor_K)) * cell->Volume();
R[cell->Integer(id)] += cell->Mean(func_rhs, cell->Real(tensor_K)) * cell->Volume();
}
BARRIER
if( m->GetProcessorRank() == 0 ) std::cout << "Matrix assemble: " << Timer()-ttt << std::endl;
......@@ -255,16 +240,16 @@ int main(int argc,char ** argv)
if( argc > 3 ) // Save the matrix and RHS if required
{
ttt = Timer();
A.Save(std::string(argv[2])); // "A.mtx"
b.Save(std::string(argv[3])); // "b.rhs"
R.GetJacobian().Save(std::string(argv[2])); // "A.mtx"
R.GetResidual().Save(std::string(argv[3])); // "b.rhs"
BARRIER
if( m->GetProcessorRank() == 0 ) std::cout << "Save matrix \"" << argv[2] << "\" and RHS \"" << argv[3] << "\": " << Timer()-ttt << std::endl;
}
ttt = Timer();
S.SetMatrix(A); // Compute the preconditioner for the original matrix
S.Solve(b,x); // Solve the linear system with the previously computted preconditioner
S.SetMatrix(R.GetJacobian()); // Compute the preconditioner for the original matrix
S.Solve(R.GetResidual(),Update); // Solve the linear system with the previously computted preconditioner
BARRIER
if( m->GetProcessorRank() == 0 )
......@@ -293,7 +278,7 @@ int main(int argc,char ** argv)
{
Storage::real old = cell->Real(phi);
Storage::real exact = cell->Mean(func, 0); // Compute the mean value of the function over the cell
Storage::real res = x[aut.GetDynamicIndex(cell->self(),iphi)];
Storage::real res = Update[aut.GetDynamicIndex(cell->self(),iphi)];
Storage::real sol = old-res;
Storage::real err = fabs (sol - exact);
if (err > local_err_C) local_err_C = err;
......
Examples/ADMFD/main.cpp
View file @ 7d30608c
......@@ -99,6 +99,9 @@ int main(int argc,char ** argv)
Tag tag_W; // Gradient matrix acting on harmonic points on faces and returning gradient on faces
Tag tag_D; // Entries for scaling matrix D
Tag tag_LF; // Coefficients of linearly computed fluxes on faces, 2 of them per internal face
Tag tag_MF; // Fluxes after inner product
tag_MF = m->CreateTag("IPFLUX",DATA_VARIABLE,FACE,NONE,1);
if( m->GetProcessorsNumber() > 1 ) //skip for one processor job
{ // Exchange ghost cells
......@@ -139,14 +142,14 @@ int main(int argc,char ** argv)
srand(1); // Randomization
tag_P = m->CreateTag("PRESSURE",DATA_REAL,CELL|FACE,NONE,1); // Create a new tag for the pressure
for(Mesh::iteratorElement e = m->BeginElement(CELL|FACE); e != m->EndElement(); ++e) //Loop over mesh cells
e->Real(tag_P) = (rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
e->Real(tag_P) = 0;//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
}
if( !tag_P.isDefined(FACE) )
{
tag_P = m->CreateTag("PRESSURE",DATA_REAL,FACE,NONE,1);
for(Mesh::iteratorElement e = m->BeginElement(FACE); e != m->EndElement(); ++e) //Loop over mesh cells
e->Real(tag_P) = (rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
e->Real(tag_P) = 0;//(rand()*1.0)/(RAND_MAX*1.0); // Prescribe random value in [0,1]
}
......@@ -197,11 +200,11 @@ int main(int argc,char ** argv)
N(k,2) = nK(0,2);
} //end of loop over faces
// formula (5.31)
IP = R*(R.Transpose()*N).Invert().first*R.Transpose(); // Consistency part
IP = R*(R.Transpose()*N).Invert(true).first*R.Transpose(); // Consistency part
// formula (5.32)
IP += (rMatrix::Unit(NF) - N*(N.Transpose()*N).Invert().first*N.Transpose())*(1.0/(static_cast<real>(NF)*vP)*(R*K.Invert().first*R.Transpose()).Trace()); //Stability part
IP += (rMatrix::Unit(NF) - N*(N.Transpose()*N).Invert(true).first*N.Transpose())*(2.0/(static_cast<real>(NF)*vP)*(R*K.Invert(true).first*R.Transpose()).Trace()); //Stability part
//assert(IP.isSymmetric()); //test positive definitiness as well!
/*
if( !IP.isSymmetric() )
{
std::cout << "unsymmetric" << std::endl;
......@@ -209,7 +212,7 @@ int main(int argc,char ** argv)
std::cout << "check" << std::endl;
(IP-IP.Transpose()).Print();
}
*/
real_array M = cell->RealArrayDV(tag_M); //access data structure for inner product matrix in mesh
M.resize(NF*NF); //resize variable array
std::copy(IP.data(),IP.data()+NF*NF,M.data()); //write down the inner product matrix
......@@ -235,25 +238,24 @@ int main(int argc,char ** argv)
//for(int k = 0; k < 3; ++k) S(k,k) = sqrt(S(k,k));
//rMatrix K = U*S*V;
rMatrix GRAD(NF,NF), NK(NF,3), R(NF,3); //big gradient matrix, co-normals, directions
GRAD.Zero();
for(int k = 0; k < NF; ++k) //loop over faces
{
aF = faces[k]->Area();
faces[k]->Centroid(yF);
faces[k]->OrientedUnitNormal(cell->self(),nF);
// assemble matrix of directions
R(k,0) = (yF[0]-xP[0])*aF;
R(k,1) = (yF[1]-xP[1])*aF;
R(k,2) = (yF[2]-xP[2])*aF;
R(k,0) = (yF[0]-xP[0]);//*aF;
R(k,1) = (yF[1]-xP[1]);//*aF;
R(k,2) = (yF[2]-xP[2]);//*aF;
// assemble matrix of co-normals
rMatrix nK = rMatrix::FromVector(nF,3).Transpose()*K;
NK(k,0) = nK(0,0);
NK(k,1) = nK(0,1);
NK(k,2) = nK(0,2);
} //end of loop over faces
GRAD = NK*(NK.Transpose()*R).Invert().first*NK.Transpose(); //stability part
//GRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert().first*R.Transpose())*(1.0/(static_cast<real>(NF)*vP)*(NK*K.Invert().first*NK.Transpose()).Trace());
GRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert().first*R.Transpose())*(2.0/(static_cast<real>(NF))*GRAD.Trace());
GRAD = NK*(NK.Transpose()*R).Invert(true).first*NK.Transpose(); //stability part
GRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert(true).first*R.Transpose())*(2.0/(static_cast<real>(NF)*vP)*(NK*K.Invert(true).first*NK.Transpose()).Trace());
//GRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert().first*R.Transpose())*(2.0/(static_cast<real>(NF))*GRAD.Trace());
real_array W = cell->RealArrayDV(tag_W); //access data structure for gradient matrix in mesh
W.resize(NF*NF); //resize the structure
std::copy(GRAD.data(),GRAD.data()+NF*NF,W.data()); //write down the gradient matrix
......@@ -276,7 +278,6 @@ int main(int argc,char ** argv)
{ //Main loop for problem solution
Automatizator aut(m); // declare class to help manage unknowns
Automatizator::MakeCurrent(&aut);
Sparse::RowMerger & merger = aut.GetMerger(); //get structure that helps matrix assembly
dynamic_variable P(aut,aut.RegisterDynamicTag(tag_P,CELL|FACE)); //register pressure as primary unknown
variable calc; //declare variable that helps calculating the value with variations
aut.EnumerateDynamicTags(); //enumerate all primary variables
......@@ -304,6 +305,8 @@ int main(int argc,char ** argv)
do
{
R.Clear(); //clean up the residual
for( int q = 0; q < m->FaceLastLocalID(); ++q ) if( m->isValidFace(q) )
m->FaceByLocalID(q)->Variable(tag_MF) = 0.0;
//First we need to evaluate the gradient at each cell for scaling matrix D
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) //loop over cells
{
......@@ -314,17 +317,22 @@ int main(int argc,char ** argv)
rMatrix GRAD(cK->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
vMatrix pF(NF,1); //vector of pressure differences on faces
vMatrix FLUX(NF,1); //computed flux on faces
vMatrix MFLUX(NF,1);
rMatrix M(cell->RealArrayDV(tag_M).data(),NF,NF); //inner product matrix
for(int k = 0; k < NF; ++k)
pF(k,0) = P(faces[k]) - P(cK);
FLUX = GRAD*pF; //fluxes on faces
for(int k = 0; k < NF; ++k) //copy the computed flux value with variations into mesh
faces[k]->VariableArray(tag_LF)[(faces[k]->BackCell() == cell)? 0 : 1] = FLUX(k,0);
MFLUX = M*FLUX;
for(int k = 0; k < NF; ++k)
faces[k]->Variable(tag_MF) += MFLUX(k,0)*(faces[k].FaceOrientedOutside(cell) ? 1 : -1);
} //end of loop over cells
//Now we need to assemble and transpose nonlinear gradient matrix
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) ) //loop over cells
{
const real eps1 = 1.0e-7;
const real eps1 = 1.0e-3;
const real eps2 = 1.0e-9;
Cell cell = m->CellByLocalID(q);
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
......@@ -348,11 +356,27 @@ int main(int argc,char ** argv)
{
variable & v = LF[(faces[k]->BackCell() == cell) ? 1 : 0]; //neighbour flux value
//single flux definition
FLUX(k,0) = u;//(soft_abs(u,eps1)+eps2)*(soft_abs(v,eps1)+eps2)/(soft_abs(u,eps1)+soft_abs(v,eps1)+2*eps2)*(soft_sign(u,eps1) + soft_sign(v,eps1)); //restore matrix of fluxes
D(k,k) = 1.0;//FLUX(k,0) / u;
/*
FLUX(k,0) = (soft_abs(u,eps1)+eps2)*(soft_abs(v,eps1)+eps2)/(soft_abs(u,eps1)+soft_abs(v,eps1)+2*eps2)*(soft_sign(u,eps1) + soft_sign(v,eps1)); //restore matrix of fluxes
if( u*v > 0 )
D(k,k) = 2.0*(soft_abs(v,eps1)+eps2)/(soft_abs(u,eps1)+soft_abs(v,eps1)+2*eps2);
else
D(k,k) = 0;
*/
//D(k,k) = FLUX(k,0) / u;
//dual flux definition
//FLUX(k,0) = u*D(k,k);
// FLUX(k,0) = ((u*u+eps2)*v+(v*v+eps2)*u)/(u*u+v*v + 2*eps2);
//D(k,k) = 1;//(v*v+eps2)/(u*u+v*v+2*eps2);
if( u*v < 0 )
D(k,k) = (2*(soft_abs(v,eps1)+eps2)/(soft_abs(u,eps1)+soft_abs(v,eps1)+2*eps2));
else
D(k,k) = eps2;
D(k,k) = 1;
FLUX(k,0) = u*D(k,k);
}
//FLUX(k,0) = faces[k]->Variable(tag_MF)*(faces[k].FaceOrientedOutside(cell) ? 1 : -1);
}
vMatrix DIV = -(D*GRAD).Transpose()*M; //cell-wise div
......@@ -379,10 +403,10 @@ int main(int argc,char ** argv)
if( tag_BC.isValid() && faces[k].HaveData(tag_BC) )
{
real_array BC = faces[k].RealArray(tag_BC);
R[index] += BC[0]*P(faces[k]) + BC[1]*DIVKGRAD(k,0) - BC[2];
R[index] -= BC[0]*P(faces[k]) + BC[1]*DIVKGRAD(k,0) - BC[2];
}
else
R[index] += DIVKGRAD(k,0);
R[index] -= DIVKGRAD(k,0);
}
}
......@@ -391,7 +415,7 @@ int main(int argc,char ** argv)
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
{
Cell cell = m->CellByLocalID(q);
if( cell->HaveData(tag_F) ) R[P.Index(cell)] -= cell->Real(tag_F)*cell->Volume();
if( cell->HaveData(tag_F) ) R[P.Index(cell)] += cell->Real(tag_F)*cell->Volume();
}
}
......@@ -403,6 +427,8 @@ int main(int argc,char ** argv)
Solver S(Solver::INNER_MPTILUC);
S.SetMatrix(R.GetJacobian());
S.SetParameterReal("drop_tolerance", 1.0e-4);
S.SetParameterReal("reuse_tolerance", 1.0e-6);
//std::fill(Update.Begin(),Update.End(),0.0);
if( S.Solve(R.GetResidual(),Update) )
{
......@@ -421,6 +447,7 @@ int main(int argc,char ** argv)
else
{
std::cout << "Unable to solve: " << S.GetReason() << std::endl;
break;
}
} while( R.Norm() > 1.0e-4 ); //check the residual norm
......@@ -437,7 +464,7 @@ int main(int argc,char ** argv)
Cell cell = m->CellByLocalID(q);
real err = cell->Real(tag_P) - cell->Real(tag_R);
real vol = cell->Volume();
if( C < err ) C = err;
if( C < fabs(err) ) C = fabs(err);
L2 += err*err*vol;
volume += vol;
cell->Real(tag_E) = err;
......@@ -453,7 +480,7 @@ int main(int argc,char ** argv)
Face face = m->FaceByLocalID(q);
real err = face->Real(tag_P) - face->Real(tag_R);
real vol = (face->BackCell()->Volume() + (face->FrontCell().isValid() ? face->FrontCell()->Volume() : 0))*0.5;
if( C < err ) C = err;
if( C < fabs(err) ) C = fabs(err);
L2 += err*err*vol;
volume += vol;
face->Real(tag_E) = err;
......
Source/Headers/inmost_dense.h
View file @ 7d30608c
......@@ -527,7 +527,7 @@ namespace INMOST
}
return ret;
}
std::pair<Matrix,bool> Invert() const
std::pair<Matrix,bool> Invert(bool print_fail = false) const
{
std::pair<Matrix,bool> ret = std::make_pair(Matrix(m,n),true);
Matrix At = Transpose(); //m by n matrix
......@@ -601,6 +601,7 @@ namespace INMOST
if( ok ) AtA(i,i) = AtA(i,i) < 0.0 ? - 1.0e-12 : 1.0e-12;
else
{
if( print_fail ) std::cout << "Failed to invert matrix" << std::endl;
ret.second = false;
return ret;
}
......@@ -679,12 +680,15 @@ namespace INMOST
}
return true;
}
Var DotProduct(const Matrix & other) const
template<typename typeB>
typename Promote<Var,typeB>::type DotProduct(const Matrix<typeB> & other) const
{
assert(n == other.n);
assert(m == other.m);
Var ret = 0.0;
for(enumerator i = 0; i < n*m; ++i) ret += space[i]*other.space[i];
typename Promote<Var,typeB>::type ret = 0.0;
for(enumerator i = 0; i < n; ++i)
for(enumerator j = 0; j < m; ++j)
ret += ((*this)(i,j))*other(i,j);
return ret;
}
Var FrobeniusNorm()
......
Source/Headers/inmost_expression.h
View file @ 7d30608c
......@@ -308,6 +308,8 @@ namespace INMOST
INMOST_DATA_REAL_TYPE & value;
Sparse::Row & entries;
public:
void Lock() {entries.Lock();}
void Unlock() {entries.Unlock();}
/// Constructor, set links to the provided value and entries
multivar_expression_reference(INMOST_DATA_REAL_TYPE & _value, Sparse::Row & _entries)
: value(_value), entries(_entries) {}
......@@ -1319,14 +1321,6 @@ namespace INMOST
for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) ret += (*it)*(*it);
return sqrt(ret);
}
void Lock(INMOST_DATA_ENUM_TYPE row)
{
jacobian[row].Lock();
}
void Unlock(INMOST_DATA_ENUM_TYPE row)
{
jacobian[row].Unlock();
}
};
}
......
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