Commit 3a09573c authored by Kirill Terekhov's avatar Kirill Terekhov

Some changes

parent cf1303c0
......@@ -95,11 +95,11 @@ int main(int argc,char ** argv)
Tag tag_K; // Diffusion tensor
Tag tag_F; // Forcing term
Tag tag_BC; // Boundary conditions
Tag tag_M; // Stiffness matrix
Tag tag_I; // Temporary local indexation
Tag tag_H; // Harmonic points
Tag tag_W; // Gradient matrix acting on harmonic points on faces and returning gradient on faces
Tag tag_W; // Gradient matrix
Tag tag_D; // Entries for scaling matrix D
Tag tag_T; // Transmissibility
Tag tag_Y; // Transverse component of tensor
if( m->GetProcessorsNumber() > 1 ) //skip for one processor job
{ // Exchange ghost cells
......@@ -148,31 +148,32 @@ int main(int argc,char ** argv)
{
tag_BC = m->GetTag("BOUNDARY_CONDITION");
//initialize unknowns at boundary
}
//initialize unknowns at boundary
for(Mesh::iteratorFace face = m->BeginFace(); face != m->EndFace(); ++face) if( face->Boundary() ) //if( face->HaveData(tag_BC) )
face->Real(tag_P) = 0.0; //create zero entry
m->ExchangeData(tag_P,CELL|FACE,0); //Synchronize initial solution with boundary unknowns
//run in a loop to identify boundary pressures so that they enter as primary variables
//3 entries are coordinates and last entry is the coefficient for FrontCell
tag_H = m->CreateTag("HARMONIC_POINT",DATA_REAL,FACE,NONE,4);
tag_H = m->CreateTag("HARMONIC_POINT_",DATA_REAL,FACE,NONE,4);
tag_Y = m->CreateTag("GAMMA",DATA_REAL,FACE,NONE,3);
tag_T = m->CreateTag("TRANSMISSIBILITY",DATA_REAL,FACE,NONE,1);
// Assemble coefficients for harmonic averaging matrix H
for(Mesh::iteratorFace fKL = m->BeginFace(); fKL != m->EndFace(); ++fKL) //go over faces
{
real_array H = fKL->RealArrayDF(tag_H); //access data structure
real & T = fKL->RealDF(tag_T);
real_array Y = fKL->RealArrayDF(tag_Y);
Cell cK = fKL->BackCell(); //get cell K from the back of the normal
Cell cL = fKL->FrontCell(); //get cell L to the front of the normal
if( !cL.isValid() ) //this is boundary face
{
fKL->Centroid(H.data()); //set center of face as harmonic point
H[3] = 0.0; //set coefficient
continue;
}
real dK, dL, lK, lL, D;
real coefK, coefL, coefQ, coefDiv;
real coefK, coefL, coefQ, coefDiv,aF;
rMatrix y(3,1); //harmonic point
rMatrix yK(3,1); //projection of the center of cell K onto surface
......@@ -181,271 +182,131 @@ int main(int argc,char ** argv)
rMatrix xL(3,1); //center of cell L
rMatrix xKL(3,1); //point on face
rMatrix nKL(3,1); //normal to face
rMatrix lKs(3,1); //reminder of the co-normal vector Kn in cell K when projection to normal is substracted, i.e. (Kn - n n.Kn)
rMatrix lLs(3,1); //reminder of the co-normal vector Kn in cell L when projection to normal is substracted, i.e. (Kn - n n.Kn)
rMatrix Ks(3,1); //reminder of the co-normal vector Kn in cell K when projection to normal is substracted, i.e. (Kn - n n.Kn)
rMatrix Ls(3,1); //reminder of the co-normal vector Kn in cell L when projection to normal is substracted, i.e. (Kn - n n.Kn)
rMatrix gamma(3,1); //transverse part of co-normal not accounted with two-point transmissbility
rMatrix KK = rMatrix::FromTensor(cK->RealArray(tag_K).data(),cK->RealArray(tag_K).size()).Transpose(); //diffusion tensor in cell K
rMatrix KL = rMatrix::FromTensor(cL->RealArray(tag_K).data(),cL->RealArray(tag_K).size()).Transpose(); //diffusion tensor in cell L
cK->Centroid(xK.data()); //retrive center of cell K
cL->Centroid(xL.data()); //retrive center of cell L
fKL->Centroid(xKL.data()); //retrive center of face
fKL->OrientedUnitNormal(cK,nKL.data()); //retrive unit normal to face
aF = fKL->Area();
lKs = KK*nKL; //get conormal in cell K
lLs = KL*nKL; //get conormal in cell L
lK = nKL.DotProduct(lKs); //find projection of conormal onto normal in cell K
lL = nKL.DotProduct(lLs); //find projection of conormal onto normal in cell L
lKs -= nKL*lK; //obtain reminder in cell K
lLs -= nKL*lL; //obtain reminder in cell L
D = -nKL.DotProduct(xKL); // Compute constant in equation for plane (nx,ny,nz,D)
dK = nKL.DotProduct(xK)+D; //compute distance to the center of cell K
dL = nKL.DotProduct(xL)+D; //compute distance to the center of cell L
D = nKL.DotProduct(xKL); // Compute constant in equation for plane (nx,ny,nz,D)
dK = D-nKL.DotProduct(xK); //compute distance to the center of cell K
if( dK < 0 ) std::cout << "incorrect sign for dK " << dK << std::endl;
if( dK*dL > 0 ) //check consistency of geometry
Ks = KK*nKL; //get conormal in cell K
lK = nKL.DotProduct(Ks); //find projection of conormal onto normal in cell K
yK = xK + nKL*dK; //compute projection of the center of cell K onto face
if( !cL.isValid() ) //this is boundary face
{
std::cout << "Cell centers are on the same side from the face" << std::endl;
T = lK / dK * aF;
//here Ks is full co-normal
gamma = Ks - (xKL - xK) * lK / dK;
y = xKL;
coefL = 0;
}
else
{
rMatrix KL = rMatrix::FromTensor(cL->RealArray(tag_K).data(),cL->RealArray(tag_K).size()).Transpose(); //diffusion tensor in cell L
dK = fabs(dK); //signs should be encorporated automatically
dL = fabs(dL);
cL->Centroid(xL.data()); //retrive center of cell L
if( dK < 1.0e-9 || dL < 1.0e-9 ) std::cout << "Cell center is located on the face" << std::endl;
Ls = KL*nKL; //get conormal in cell L
lL = nKL.DotProduct(Ls); //find projection of conormal onto normal in cell L
Ks -= nKL*lK; //obtain reminder in cell K
Ls -= nKL*lL; //obtain reminder in cell L
yK = xK + nKL*dK; //compute projection of the center of cell K onto face
yL = xL - nKL*dL; //compute projection of the center of cell L onto face
dL = D-nKL.DotProduct(xL); //compute distance to the center of cell L
if( dL > 0 ) std::cout << "incorrect sign for dK " << dK << std::endl;
//compute coefficients for harmonic point
coefDiv = (lL*dK+lK*dL);
coefL = lL*dK/coefDiv;
coefK = lK*dL/coefDiv;
coefQ = dK*dL/coefDiv;
if( dK*dL > 0 ) //check consistency of geometry
std::cout << "Cell centers are on the same side from the face" << std::endl;
dL = -dL;
y = yK*coefK + yL*coefL + (lKs - lLs)*coefQ; //calculate position of harmonic point
if( dK < 1.0e-9 || dL < 1.0e-9 ) std::cout << "Cell center is located on the face, dK " << dK << " dL " << dL << std::endl;
std::copy(y.data(),y.data()+3,H.data()); //store position of harmonic point
H[3] = coefL; //store multiplier for pressure in FrontCell, the other is coefK = 1 - coefL
} //end of loop over faces
tag_M = m->CreateTag("INNER_PRODUCT",DATA_REAL,CELL,NONE);
//assemble inner product matrix M acting on faces for each cell
for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell)
{
real xP[3]; //center of the cell
real nF[3]; //normal to the face
real aF, vP = cell->Volume(); //area of the face
cell->Centroid(xP); //obtain cell center
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
enumerator NF = (enumerator)faces.size(), k; //number of faces;
yL = xL - nKL*dL; //compute projection of the center of cell L onto face
rMatrix IP(NF,NF), N(NF,3), R(NF,3); //matrices for inner product
//compute coefficients for harmonic point
coefDiv = (lL*dK+lK*dL);
coefL = lL*dK/coefDiv;
coefK = lK*dL/coefDiv;
coefQ = dK*dL/coefDiv;
//assemble matrices R and N
k = 0;
for(ElementArray<Face>::iterator face = faces.begin(); face != faces.end(); ++face)
{
aF = face->Area();
real_array yF = face->RealArrayDF(tag_H); //point on face is harmonic point
face->OrientedUnitNormal(cell->self(),nF);
// assemble R matrix, follows chapter 3.4.3 in the book "Mimetic Finite Difference Method for Elliptic Problems" by Lourenco et al
R(k,0) = (yF[0]-xP[0])*vP;
R(k,1) = (yF[1]-xP[1])*vP;
R(k,2) = (yF[2]-xP[2])*vP;
// assemble N matrix, follow chapter 3.4. in the book and projection operator definition from
// chapter 2.3.1 in "Mimetic scalar products of differential forms" by Brezzi et al
N(k,0) = nF[0]*aF;
N(k,1) = nF[1]*aF;
N(k,2) = nF[2]*aF;
k++;
} //end of loop over faces
T = lL*lK/coefDiv*aF;
//inner product definition from Corollary 4.2, "Mimetic scalar products of differential forms" by Brezzi et al
IP = R*(R.Transpose()*N).Invert().first*R.Transpose(); // Consistency part
IP += (rMatrix::Unit(NF) - N*(N.Transpose()*N).Invert().first*N.Transpose())*(2.0/static_cast<real>(NF)*IP.Trace()); //Stability part
//assert(IP.isSymmetric()); //test positive definitiness as well!
/*
if( !IP.isSymmetric() )
{
std::cout << "unsymmetric" << std::endl;
IP.Print();
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
} //end of loop over cells
gamma = (yK-yL)*T + coefL*Ks + coefK*Ls;
y = yK*coefK + yL*coefL + (Ks - Ls)*coefQ; //calculate position of harmonic point
}
std::copy(gamma.data(),gamma.data()+3,Y.data()); //store co-normal reminder
std::copy(y.data(),y.data()+3,H.data()); //store position of harmonic point
H[3] = coefL; //store multiplier for pressure in FrontCell, the other is coefK = 1 - coefL
} //end of loop over faces
tag_W = m->CreateTag("GRAD",DATA_REAL,CELL,NONE);
tag_D = m->CreateTag("DIAG",DATA_VARIABLE,CELL,NONE);
//tag_D = m->CreateTag("DIAG",DATA_VARIABLE,CELL,NONE);
//Assemble all-positive gradient matrix W on cells
for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell)
{
//real_array D = cell->VariableArrayDV(tag_D); //resize scaling matrix D for the future use
real_array W = cell->RealArrayDV(tag_W); //access data structure for gradient matrix in mesh
real xP[3]; //center of the cell
real nF[3]; //normal to the face
real vP = cell->Volume();
cell->Centroid(xP);
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
enumerator NF = (enumerator)faces.size(), k; //number of faces;
dynarray<real,64> NG(NF,0), G(NF,0); // count how many times gradient was calculated for face
int NF = (int)faces.size(); //number of faces;
rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),cell->RealArrayDF(tag_K).size()); //get permeability for the cell
rMatrix GRAD(NF,NF), NK(NF,3), R(NF,3); //big gradient matrix, co-normals, directions
rMatrix GRADF(NF,NF), NKF(NF,3), RF(NF,3); //corresponding matrices with masked negative values when multiplied by individual co-normals
rMatrix GRADFstab(NF,NF);
rMatrix H(NF,NF);//,T(NF,NF);
GRAD.Zero();
k = 0;
for(ElementArray<Face>::iterator face = faces.begin(); face != faces.end(); ++face) //loop over faces
H.Zero();
//T.Zero();
for(int k = 0; k < NF; ++k) //loop over faces
{
real_array yF = face->RealArrayDF(tag_H); //point on face is harmonic point
face->OrientedUnitNormal(cell->self(),nF);
real aF = faces[k].Area();
real_array yF = faces[k]->RealArrayDF(tag_H); //point on face is harmonic point
real_array gammaF = faces[k]->RealArrayDF(tag_Y);
faces[k]->OrientedUnitNormal(cell->self(),nF);
// assemble matrix of directions
R(k,0) = (yF[0]-xP[0]);
R(k,1) = (yF[1]-xP[1]);
R(k,2) = (yF[2]-xP[2]);
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) = gammaF[0];
//NK(k,1) = gammaF[1];
//NK(k,2) = gammaF[2];
NK(k,0) = nK(0,0);
NK(k,1) = nK(0,1);
NK(k,2) = nK(0,2);
k++;
H(k,k) = aF * (faces[k].BackCell() == cell->self() ? 1.0-yF[3] : yF[3]);
//T(k,k) = faces[k]->RealDF(tag_T);
} //end of loop over faces
GRAD = NK*(NK.Transpose()*R).Invert().first*NK.Transpose(); //stability part
//std::cout << "consistency" << std::endl;
//GRADF.Print();
GRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert().first*R.Transpose())*(2.0/NF*GRADF.Trace());
/*
std::cout << "stability" << std::endl;
GRADFstab.Print();
GRADF += GRADFstab;
std::cout << "consistency+stability" << std::endl;
GRADF.Print();
k = 0;
for(ElementArray<Face>::iterator face = faces.begin(); face != faces.end(); ++face) //loop over faces
{
real nK[3] = {NK(k,0),NK(k,1),NK(k,2)};
real N = 0.0;
for(enumerator l = 0; l < NF; ++l) //run through faces
{
if( nK[0]*NK(l,0) + nK[1]*NK(l,1) + nK[2]*NK(l,2) > 0.0 && nK[0]*R(l,0) + nK[1]*R(l,1) + nK[2]*R(l,2) > 0.0 ) //all entries are positive
{
NKF(l,0) = NK(l,0);
NKF(l,1) = NK(l,1);
NKF(l,2) = NK(l,2);
RF(l,0) = R(l,0);
RF(l,1) = R(l,1);
RF(l,2) = R(l,2);
G[l] = 1; //mask active rows
NG[l]++; //count number of equations
N++; //count number of active rows
}
else
{
NKF(l,0) = NKF(l,1) = NKF(l,2) = 0.0;
RF(l,0) = RF(l,1) = RF(l,2) = 0.0;
G[l] = 0; //mask active rows
}
}
//NKF.Print();
//RF.Print();
GRADF.Zero();
GRADF = NKF*(NKF.Transpose()*RF).Invert()*NKF.Transpose(); //stability part
std::cout << "consistency" << std::endl;
GRADF.Print();
GRADFstab.Zero();
GRADFstab = (FromDiagonal(G.data(),NF) - RF*(RF.Transpose()*RF).Invert()*RF.Transpose())*(2.0/N*GRADF.Trace());
real alpha = 1.0;
for(enumerator i = 0; i < NF; ++i)
{
for(enumerator j = 0; j < NF; ++j)
if( GRADF(i,j) + alpha*GRADFstab(i,j) < 0.0 )
{
alpha = std::min(alpha,-GRADF(i,j)/GRADFstab(i,j));
}
}
std::cout << "stability" << std::endl;
GRADFstab.Print();
std::cout << "alpha " << alpha << std::endl;
GRADF += GRADFstab*alpha; //consistency part
std::cout << "consistency+stability" << std::endl;
GRADF.Print();
GRAD += GRADF;
//GRAD.Print();
k++;
} //end of loop over faces
GRAD = GRAD*FromDiagonalInverse(NG.data(),NF); //normalize gradients if we computed them multiple times
//std::cout << "Gradient matrix: " << std::endl;
//GRAD.Print();
*/
real_array W = cell->RealArrayDV(tag_W); //access data structure for gradient matrix in mesh
GRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert().first*R.Transpose())*(2.0/(vP*NF)*GRAD.Trace());
GRAD = H*GRAD*H;// + T;
W.resize(NF*NF); //resize the structure
std::copy(GRAD.data(),GRAD.data()+NF*NF,W.data()); //write down the gradient matrix
cell->VariableArrayDV(tag_D).resize(NF); //resize scaling matrix D for the future use
} //end of loop over cells
if( m->HaveTag("FORCE") ) //Is there force on the mesh?
{
Tag tag_F0 = m->GetTag("FORCE"); //initial force
assert(tag_F0.isDefined(CELL)); //assuming it was defined on cells
tag_F = m->CreateTag("RECOMPUTED_FORCE",DATA_REAL,CELL,NONE,1); //recomputed force on face
// compute expression H^T M H f, where H is harmonic averaging and M is inner product over faces
for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell)
{
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
enumerator NF = (enumerator)faces.size(), k = 0;
Cell cK = cell->self();
real gK = cK->Real(tag_F0), gL; //force on current cell and on adjacent cell
rMatrix M(cK->RealArrayDF(tag_M).data(),NF,NF); //inner product matrix
rMatrix gF(NF,1); //vector of forces on faces
rMatrix MgF(NF,1); //vector of forces multiplied by inner product matrix
//first average force onto faces with harmonic mean
k = 0;
for(ElementArray<Face>::iterator face = faces.begin(); face != faces.end(); ++face) //go over faces
{
real_array H = face->RealArrayDF(tag_H); //access structure for harmonic averaging
if( !face->Boundary() ) //internal face
{
Cell cL = cK->Neighbour(face->self());
real aL = (cK == face->BackCell()) ? H[3] : 1-H[3]; //harmonic coefficient for neighbour cell
real aK = 1 - aL; //harmonic coefficient for current cell
gL = cL->Real(tag_F0); //retrive force on adjacent cell
gF(k,0) = gK*aK + gL*aL; //harmonic averaging
}
else gF(k,0) = gK;
k++;
} //end of loop over faces
MgF = M*gF; //multiply by inner product matrix
//now apply H^T for current cell only
real & rgK = cell->Real(tag_F); //access force in current cell
rgK = 0.0;
k = 0;
for(ElementArray<Face>::iterator face = faces.begin(); face != faces.end(); ++face) //go over faces
{
real_array H = face->RealArrayDF(tag_H); //access structure for harmonic averaging
if( !face->Boundary() ) //internal face
{
Cell cL = cell->Neighbour(face->self());
real aL = (cK == face->BackCell()) ? H[3] : 1-H[3]; //harmonic coefficient for neighbour cell
real aK = 1 - aL; //harmonic coefficient for current cell
real & rgL = face->FrontCell()->RealDF(tag_F); // access force in adjacent cell
//harmonic distribution
rgK += MgF(k,0)*aK; // application of H^T with inner product in current cell
rgL += MgF(k,0)*aL; // application of H^T with inner product in current cell
}
else rgK += MgF(k,0);
k++;
} //end of loop over faces
} //end of loop over cells
m->ExchangeData(tag_F,CELL,0); //synchronize force
} // end of force
} //end of initialize data
......@@ -455,180 +316,121 @@ 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
Residual R("",aut.GetFirstIndex(),aut.GetLastIndex());
Sparse::Matrix Jacobian("",aut.GetFirstIndex(),aut.GetLastIndex()); //matrix for jacobian
Sparse::Vector Update("",aut.GetFirstIndex(),aut.GetLastIndex()); //vector for update
Sparse::Vector Residual("",aut.GetFirstIndex(),aut.GetLastIndex()); //vector for residual
real Residual_norm = 0.0;
int nit = 0;
do
{
std::fill(Residual.Begin(),Residual.End(),0.0); //clean up the residual
R.Clear(); //clean up the residual
//First we need to evaluate the gradient at each cell for scaling matrix D
for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell) //loop over cells
{
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
enumerator NF = (enumerator)faces.size(), k = 0;
Cell cK = cell->self();
variable pK = P(cK), pL; //pressure for back and front cells
rMatrix GRAD(cK->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
int NF = (int)faces.size();
rMatrix GRAD(cell->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
vMatrix pF(NF,1); //vector of pressures on faces
vMatrix FLUX(NF,1); //computed flux on faces
for(ElementArray<Face>::iterator face = faces.begin(); face != faces.end(); ++face)
vMatrix FLUX(NF,1); //computed fluxes
for(int k = 0; k < NF; ++k)
{
real_array H = face->RealArrayDF(tag_H); //access structure for harmonic averaging
if( !face->Boundary() ) //internal face
{
Cell cL = cK->Neighbour(face->self()); //
real aL = (cK == face->BackCell()) ? H[3] : 1-H[3]; //harmonic coefficient for neighbour cell
real aK = 1 - aL; //harmonic coefficient for current cell
pL = P(face->FrontCell()); //retrive force on adjacent cell
pF(k,0) = pK*aK + pL*aL;
}
else if( face->HaveData(tag_P) ) //boundary condition
pF(k,0) = P(face->self()); //get pressure on boundary
if( faces[k]->Boundary() )
pF(k,0) = P(faces[k]) - P(cell->self());
else
pF(k,0) = pK; //get pressure in current cell
pF(k,0) -= pK; //substract pressure in the center to get differences
k++;
pF(k,0) = P(faces[k]->FrontCell()) - P(cell->self());
}
FLUX = GRAD*pF; //fluxes on faces
var_array D = cell->VariableArrayDV(tag_D); //access data structure on mesh for variations
for(enumerator l = 0; l < NF; ++l) D[l] = FLUX(l,0); //copy the computed flux with variations into mesh
} //end of loop over cells
//Now we need to assemble and transpose nonlinear gradient matrix
for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell) //loop over cells
{
ElementArray<Face> facesK = cell->getFaces(), facesL; //obtain faces of the cell
enumerator NF = (enumerator)facesK.size(), k, l;
Cell cK = cell->self();
variable pK = P(cK), pL; //pressure for back and front cells
rMatrix GRAD(cK->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
vMatrix D(NF,NF); //Matrix for diagonal
vMatrix pF(NF,1); //vector of pressures on faces
D.Zero();
//assemble diagonal matrix, loop through faces and access corresponding entries
k = 0;
for(ElementArray<Face>::iterator faceK = facesK.begin(); faceK != facesK.end(); ++faceK) //loop over faces of current cell
for(int k = 0; k < NF; ++k)
{
if( !faceK->Boundary() ) //face is on boundary
R[P.Index(cell->self())] += FLUX(k,0);
if( faces[k]->Boundary() )
{
variable DK(cK->VariableArrayDV(tag_D)[k]);
Cell cL = cK->Neighbour(faceK->self()); //retrive neighbour
facesL = cL->getFaces(); //retrive faces of other cells
l = 0; //count face number
for(ElementArray<Face>::iterator faceL = facesL.begin(); faceL != facesL.end(); ++faceL) //loop over faces of other cell
if( faces[k]->HaveData(tag_BC) )
{
if( faceK->self() == faceL->self() ) //faces match - get flux entry
{
variable DL(cL->VariableArrayDV(tag_D)[l]); //retrive flux entry as variable
variable absDK = sqrt(DK*DK+reg_abs) + reg_div;
variable absDL = sqrt(DL*DL+reg_abs) + reg_div;
D(k,k) = absDL/(absDK+absDL); // compute absolute value and write to D matrix
break;
}
l++; //advance the face number
real_array bc = faces[k]->RealArray(tag_BC);
R[P.Index(faces[k])] -= bc[0]*P(faces[k]) + bc[1]*FLUX(k,0) - bc[2];
}
else R[P.Index(faces[k])] -= FLUX(k,0);
}
else
{
D(k,k) = 1; //should somehow use the gradient defined on boundary
}
real_array H = faceK->RealArrayDF(tag_H); //access structure for harmonic averaging
if( !faceK->Boundary() ) //internal face
{
Cell cL = cK->Neighbour(faceK->self()); //
real aL = (cK == faceK->BackCell()) ? H[3] : 1-H[3]; //harmonic coefficient for neighbour cell
real aK = 1 - aL; //harmonic coefficient for current cell
pL = P(faceK->FrontCell()); //retrive force on adjacent cell
pF(k,0) = pK*aK + pL*aL;
}
else if( faceK->HaveData(tag_P) ) //boundary condition
pF(k,0) = P(faceK->self()); //get pressure on boundary
else
pF(k,0) = pK; //get pressure in current cell
pF(k,0) -= pK; //substract pressure in the center to get differences
k++;
R[P.Index(faces[k]->FrontCell())] -= FLUX(k,0);
}
vMatrix DGRAD = D*GRAD; //multiply the matrix with diagonal
rMatrix M(cK->RealArrayDV(tag_M).data(),NF,NF); //inner product matrix
vMatrix DIVGRADHp = DGRAD.Transpose()*M*DGRAD*pF; //cell-wise divgradp
//redistribute divgradp onto cells by applying H^T
enumerator indK = P.Index(cK->self()),indL;
Storage::real & ResidK = Residual[indK];
Sparse::Row & JacK = Jacobian[indK]; //access force in current cell
k = 0;
for(ElementArray<Face>::iterator face = facesK.begin(); face != facesK.end(); ++face) //go over faces
{
real_array H = face->RealArrayDF(tag_H); //access structure for harmonic averaging
if( !face->Boundary() ) //internal face
{
Cell cL = cell->Neighbour(face->self());
real aL = (cK == face->BackCell()) ? H[3] : 1-H[3]; //harmonic coefficient for neighbour cell
real aK = 1 - aL; //harmonic coefficient for current cell
indL = P.Index(face->FrontCell());
Storage::real & ResidL = Residual[indL];
Sparse::Row & JacL = Jacobian[indL]; // access force in adjacent cell
//harmonic distribution
merger.Merge(JacK,1.0,JacK,-aK,DIVGRADHp(k,0).GetRow());// application of H^T with inner product in current cell
ResidK += DIVGRADHp(k,0).GetValue()*aK;
merger.Merge(JacL,1.0,JacL,-aL,DIVGRADHp(k,0).GetRow());// application of H^T with inner product in current cell
ResidL += DIVGRADHp(k,0).GetValue()*aL;
}
else
{
indL = P.Index(face->self());
Storage::real & ResidL = Residual[indL];
Sparse::Row & JacL = Jacobian[indL]; // access force in adjacent cell
merger.Merge(JacK,1.0,JacK,-1.0,DIVGRADHp(k,0).GetRow());
ResidK += DIVGRADHp(k,0).GetValue();
merger.Merge(JacL,1.0,JacL,-1.0,DIVGRADHp(k,0).GetRow());
ResidL += DIVGRADHp(k,0).GetValue();
}
k++;
} //end of loop over faces
}
} //end of loop over cells
if( tag_F.isValid() )
{
for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell)
Residual[P.Index(cell->self())] += cell->RealDF(tag_F);
R[P.Index(cell->self())] += cell->RealDF(tag_F)*cell->Volume();
}
R.Rescale();
Residual_norm = 0.0;
for(Sparse::Vector::iterator it = Residual.Begin(); it != Residual.End(); ++it) Residual_norm += (*it)*(*it); //compute residual norm on current iteration
std::cout << "Nonlinear residual: " << R.Norm() << std::endl;
std::cout << "Nonlinear residual: " << Residual_norm << std::endl;
if( R.Norm() < 1.0e-4 ) break;
Solver S(Solver::INNER_MPTILUC);
S.SetMatrix(Jacobian);
S.SetParameterReal("drop_tolerance",1.0e-3);
S.SetParameterReal("reuse_tolerance",1.0e-6);
S.SetParameterReal("drop_tolerance",1.0e-2);
S.SetParameterReal("reuse_tolerance",1.0e-4);
S.SetParameterEnum("gmres_substeps",4);
if( S.Solve(Residual,Update) )
R.GetJacobian().Save("jacobian.mtx");
R.GetResidual().Save("residual.mtx");
S.SetMatrix(R.GetJacobian());
if( S.Solve(R.GetResidual(),Update) )
{
for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell)
cell->Real(tag_P) += Update[P.Index(cell->self())];
cell->Real(tag_P) -= Update[P.Index(cell->self())];
for(Mesh::iteratorFace face = m->BeginFace(); face != m->EndFace(); ++face) if( face->HaveData(tag_P) )
face->Real(tag_P) -= Update[P.Index(face->self())];
}
else
{
std::cout << "Unable to solve: " << S.GetReason() << std::endl;
break;
}
} while( Residual_norm > 1.0e-4 ); //check the residual norm
} while( R.Norm() > 1.0e-4 && nit < 10 ); //check the residual norm
std::cout << "Solved problem in " << Timer() - ttt << " seconds with " << nit << " iterations " << std::endl;
if( m->HaveTag("REFERENCE_SOLUTION") )
{
Tag tag_E = m->CreateTag("ERRROR",DATA_REAL,CELL,NONE,1);
Tag tag_R = m->GetTag("REFERENCE_SOLUTION");
real C, L2, volume;
C = L2 = volume = 0.0;
for( int q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
{
Cell cell = m->CellByLocalID(q);
real err = cell->Real(tag_P) - cell->Real(tag_R);
real vol = cell->Volume();
if( C < fabs(err) ) C = fabs(err);
L2 += err*err*vol;
volume += vol;
cell->Real(tag_E) = err;
}
L2 = sqrt(L2/volume);
std::cout << "Error on cells, C-norm " << C << " L2-norm " << L2 << std::endl;
}
m->Save("out.pmf");
if( m->GetProcessorsNumber() == 1 )
m->Save("out.vtk");
else
m->Save("out.pvtk");
}
delete m; //clean up the mesh
......
......@@ -40,97 +40,97 @@ namespace INMOST
///Definition for data type of topology error or event. This type may be extended later to 64 bits
/// to accomodate more topological errors or events.
typedef INMOST_DATA_ENUM_TYPE TopologyCheck;
///Throw TopologyError exception on error.
static const TopologyCheck THROW_EXCEPTION = 0x00000001;