Commit 33e6e3dc authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

some updates

parent 5d8217ea
......@@ -110,6 +110,7 @@ ConvectionDiffusion::ConvectionDiffusion(Mesh * _m, Tag _tag_U, Tag _tag_K, Tag
KD(3,3), //difference of tensors
gammaK(1,3), //transversal part of co-normal at cell K
gammaL(1,3), //transversal part of co-normal at cell L
gamma(1,3), //common transversal part for flux
iT(3,3), //heterogeneous interpolation tensor
iC(1,3) //heterogeneous interpolation correction
;
......@@ -121,7 +122,8 @@ ConvectionDiffusion::ConvectionDiffusion(Mesh * _m, Tag _tag_U, Tag _tag_K, Tag
dK, //distance from center to interface at cell K
dL, //distance from center to interface at cell L
lambdaK, //projection of co-normal onto normal at cell K
lambdaL //projection of co-normal onto normal at cell L
lambdaL, //projection of co-normal onto normal at cell L
div //divider
;
const real eps = degenerate_diffusion_regularization;
Cell cK, cL, cU;
......@@ -155,6 +157,7 @@ ConvectionDiffusion::ConvectionDiffusion(Mesh * _m, Tag _tag_U, Tag _tag_K, Tag
KK = rMatrix::FromTensor(cK.RealArray(tag_K).data(),cK.RealArray(tag_K).size());//.Transpose();
KKn = nKL*KK;
lambdaK = nKL.DotProduct(KKn);
if( lambdaK < 0 ) std::cout << __FILE__ << ":" << __LINE__ << " lambdaK=" << lambdaK << std::endl;
gammaK = KKn - lambdaK*nKL;
//Diffusion part
......@@ -169,6 +172,7 @@ ConvectionDiffusion::ConvectionDiffusion(Mesh * _m, Tag _tag_U, Tag _tag_K, Tag
KL = rMatrix::FromTensor(cL.RealArray(tag_K).data(),cL.RealArray(tag_K).size());//.Transpose();
KLn = nKL*KL;
lambdaL = nKL.DotProduct(KLn);
if( lambdaL < 0 ) std::cout << __FILE__ << ":" << __LINE__ << " lambdaL=" << lambdaL << std::endl;
gammaL = KLn - lambdaL*nKL;
......@@ -176,16 +180,36 @@ ConvectionDiffusion::ConvectionDiffusion(Mesh * _m, Tag _tag_U, Tag _tag_K, Tag
R = 0;
if( split_diffusion )
{
uK = A * (lambdaK*lambdaL*(yK-yL) + lambdaL*dK*gammaK + lambdaK*dL*gammaL)/(dK*lambdaL+dL*lambdaK + 1.0e-100);
uL = uK;
T = A * lambdaK*lambdaL/(dK*lambdaL+dL*lambdaK + 1.0e-100);
div = (dK*lambdaL+dL*lambdaK);
if( div > 0 )
{
T = A * lambdaK*lambdaL/div;
gamma = A * (lambdaK*lambdaL*(yK-yL) + lambdaL*dK*gammaK + lambdaK*dL*gammaL)/div;
uK = gamma;
uL = gamma;
}
else if( div < 0 )
{
T = -A * lambdaK*lambdaL/div;
gamma = A * (lambdaK*lambdaL*(yK-yL) + lambdaL*dK*gammaK + lambdaK*dL*gammaL)/div;
uK = 2*A*KKn - gamma;
uL = 2*A*KLn - gamma;
}
else
{
uK = A * KKn;
uL = A * KLn;
T = 0;
}
}
if( !split_diffusion || T < 0 ) //un-splitted diffusion
else //un-splitted diffusion
{
uK = A * KKn;
uL = A * KLn;
T = 0;
}
if( T < 0 ) std::cout << __FILE__ << ":" << __LINE__ << " T=" << T << std::endl;
//internal
......@@ -199,14 +223,31 @@ ConvectionDiffusion::ConvectionDiffusion(Mesh * _m, Tag _tag_U, Tag _tag_K, Tag
//don't forget the area!
if( split_diffusion )
{
uK = KKn - ((xKL - xK) * lambdaK / dK);
T = lambdaK / dK;
if( dK > 0 )
{
T = lambdaK / dK;
gamma = KKn - ((xKL - xK) * lambdaK / dK);
uK = gamma;
}
else if( dK < 0 )
{
T = -lambdaK / dK;
gamma = KKn - ((xKL - xK) * lambdaK / dK);
uK = 2*KKn - gamma;
}
else
{
T = 0;
uK = KKn;
}
}
if( !split_diffusion || T < 0 ) //un-splitted diffusion
else //un-splitted diffusion
{
uK = KKn;
T = 0;
}
if( T < 0 ) std::cout << __FILE__ << ":" << __LINE__ << " T=" << T << std::endl;
real bcconds[3] = {0.0,1.0,0.0}; //pure neumann boundary condition
if( tag_BC.isValid() && fKL.HaveData(tag_BC) ) //are there boundary conditions on face?
......@@ -219,9 +260,9 @@ ConvectionDiffusion::ConvectionDiffusion(Mesh * _m, Tag _tag_U, Tag _tag_K, Tag
}
//account for boundary conditions
R = A*T*bcconds[2]/(bcconds[0] + bcconds[1]*T + 1.0e-100);
uK *=A* bcconds[0]/(bcconds[0] + bcconds[1]*T + 1.0e-100);
T *=A* bcconds[0]/(bcconds[0] + bcconds[1]*T + 1.0e-100);
R = A*T*bcconds[2]/(bcconds[0] + bcconds[1]*T);
uK *=A* bcconds[0]/(bcconds[0] + bcconds[1]*T);
T *=A* bcconds[0]/(bcconds[0] + bcconds[1]*T);
//on BC
if( uK.FrobeniusNorm() > 0 )
......@@ -395,7 +436,7 @@ ConvectionDiffusion::ConvectionDiffusion(Mesh * _m, Tag _tag_U, Tag _tag_K, Tag
bulk flag;
real_array v;
#if defined(USE_OMP)
#pragma omp for
#pragma omp for schedule(dynamic)
#endif
for(integer q = 0; q < m->CellLastLocalID(); ++q ) if( m->isValidCell(q) )
{
......
......@@ -34,18 +34,18 @@ typedef Storage::reference_array ref_array;
const bool use_reference_solution = false;
//SchemeType scheme_type = MPFA;
SchemeType scheme_type = MPFA;
//SchemeType scheme_type = NMPFA_QUADRATIC_CORRECTED;
//SchemeType scheme_type = NMPFA_VAN_ALBADA_CORRECTED;
//SchemeType scheme_type = NMPFA_VAN_ALBADA_CORRECTED;
SchemeType scheme_type = NMPFA_VAN_ALBADA;
//SchemeType scheme_type = NMPFA_VAN_ALBADA;
//SchemeType scheme_type = NMPFA_VAN_LEER;
//SchemeType scheme_type = NTPFA_PICARD;
//SchemeType scheme_type = NTPFA;
OutputMeshFormat output_format = PMF; //current choice of mesh format
//scheme behavior
bool split_diffusion = false; //diffusion is represented by two-point part plus correction
bool split_diffusion = true; //diffusion is represented by two-point part plus correction
bool joint_advection_diffusion = true; //nonlinear weighting is performed for advection and diffusion simultaneously
//time stepping constants
//#define FRACTIONAL //tvd-limited fractional time-stepping scheme, bacwards euler if commented
......@@ -60,11 +60,11 @@ bool check_div = false; //check that the velocity is divergence free.
double nonlinear_iterations = 100; //maximal number of iterations
double nonlinear_abs_tolerance = 1.0e-6; //absolute tolerance
double nonlinear_rel_tolerance = 1.0e-5; //relative tolerance
double regularization = 1.0e-9; //eps regularization in |x| = sqrt(x*x+eps)
double degenerate_diffusion_regularization = 1.0e-8; //eps regularization when diffusion goes to zero
double regularization = 1.0e-32; //eps regularization in |x| = sqrt(x*x+eps)
double degenerate_diffusion_regularization = 1.0e-32; //eps regularization when diffusion goes to zero
//constants for stencils
int max_layers = 3; //maximum number of layers to be considered in interpolation
int max_layers = 5; //maximum number of layers to be considered in interpolation
ElementType bridge_layers = FACE; //used as bridge elements for layers
......@@ -634,10 +634,10 @@ int main(int argc,char ** argv)
Solver S(linear_solver_type);
S.SetParameterReal("relative_tolerance", 1.0e-14);
S.SetParameterReal("absolute_tolerance", 1.0e-12);
S.SetParameterReal("relative_tolerance", 1.0e-18);
S.SetParameterReal("absolute_tolerance", 1.0e-15);
S.SetParameterReal("drop_tolerance", 1.0e-2);
S.SetParameterReal("reuse_tolerance", 1.0e-4);
S.SetParameterReal("reuse_tolerance", 1.0e-3);
S.SetMatrix(R.GetJacobian());
std::fill(Update.Begin(),Update.End(),0.0);
if( S.Solve(R.GetResidual(),Update) )
......
......@@ -2,7 +2,7 @@
using namespace INMOST;
#define HARMONIC_VERSION
//#define HARMONIC_VERSION
//shortcuts
typedef Storage::real real;
......@@ -43,6 +43,7 @@ bool find_stencils(Cell cK,
//remember condition number
std::vector<real> mincond(compute.size());
std::vector<real> minsum(compute.size());
std::vector<integer> minwgt(compute.size());
//fill data
......@@ -50,6 +51,7 @@ bool find_stencils(Cell cK,
{
minwgt[k] = max_layers*3;
mincond[k] = 1.0e+100;
minsum[k] = 1.0e+100;
compute[k].computed = false;
}
......@@ -387,10 +389,11 @@ bool find_stencils(Cell cK,
(rMatrix(compute[l].v,1,3)*V*Sinv*U.Transpose()).Print();
}
}
//check factorization
//if( (I - rMatrix::Unit(3)).FrobeniusNorm() < 1.0e-8 )
{
for(integer l = 0; l < (integer)compute.size(); ++l) if( cond < mincond[l] && combination_wgt < minwgt[l] )
for(integer l = 0; l < (integer)compute.size(); ++l) //if( cond < mincond[l] )//&& combination_wgt < minwgt[l] )
{
//current vector
v = rMatrix(compute[l].v,1,3);
......@@ -404,19 +407,33 @@ bool find_stencils(Cell cK,
coef_sum = 0.0;
for(integer q = 0; q < 3; ++q)
coef_sum += fabs(coef(0,q));
for(integer q = 0; q < 3; ++q)
{
if( coef(0,q) < 0.0 )
{
if( coef(0,q)/coef_sum < -1.0e-5 )
nonnegative = false;
else //truncate small values
coef(0,q) = 0.0;
//else //truncate small values
// coef(0,q) = 0.0;
}
}
//it's positive and condition number is improved
if( nonnegative )
if( nonnegative && coef_sum < minsum[l] )
{
/*
if( (I - rMatrix::Unit(3)).FrobeniusNorm() > 1.0e-8 )
{
#pragma omp critical
{
std::cout << "deviation from unit: " << std::endl;
(I - rMatrix::Unit(3)).Print();
std::cout << "deviation from conormal: " << std::endl;
(v-coef*A).Print();
}
}
*/
//remember stencil
if( !compute[l].computed )
{
......@@ -437,6 +454,8 @@ bool find_stencils(Cell cK,
compute[l].nonnegative = nonnegative;
*compute[l].rhs = compute[l].sign*(bndrhs[C[0]]*coef(0,0) + bndrhs[C[1]]*coef(0,1) + bndrhs[C[2]]*coef(0,2));
mincond[l] = cond;
minsum[l] = coef_sum;
minwgt[l] = combination_wgt;
} //if nonnegative
} //check coefficients
} //loop over stencils
......
......@@ -2,7 +2,7 @@
using namespace INMOST;
bool output_file = false;
bool output_file = true;
bool balance_mesh = true;
bool balance_mesh_refine = false;
bool balance_mesh_coarse = false;
......
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