/*
$Id: adaptnl.cc,v 1.2 1996/10/04 15:06:22 roitzsch Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "adaptnl.h"
#include "utils.h"
#include "numerics.h"
#include "problemnl.h"
#include "nonlin.h"
#include "dirichletA.h" // !!!
#include "elements.h"
#include "triang.h"
#include "int.h"
#include "linsystem.h"
#include "precondnl.h"
#include "sysmatml.h"
#include "family.h"
#include "cmdpars.h"
extern CmdPars Cmd;
//-------------------------------------------------------------------------
RK1:: RK1(Element* elementDLY0, Element* elementLag0, Interface* interface0)
: DLY(elementDLY0, interface0), lagrangeElementDLY(elementLag0),
doRefine(1)
{ }
//-------------------------------------------------------------------------
RK1:: ~RK1() { delete lagrangeElementDLY; }
//-------------------------------------------------------------------------
Bool RK1:: estimateError(Problem& problem)
{
if (infoErrorEstimator) cout << "\n\tError Est.: RK1";
distributeError(problem);
enforcedPMediaRefinement(problem);
return True;
}
//-------------------------------------------------------------------------
void RK1:: distributeError(Problem& problem)
{
int i;
Vector<Num> error(1), uQ(1), ADiag(1);
solveQuadDefectProblem(problem, error, uQ, ADiag);
interfaceDLY->solToHB(error);
EnergyError = 0.0;
for (i=problem.u.high()+1; i<=error.high(); ++i)
{
error[i] = error[i] * error[i] * ADiag[i]; // error on edge
EnergyError += Abs(error[i]); // sum up global Error
}
EnergyError *= 0.5;
distributeEdgeErrors(problem, error);
}
//-------------------------------------------------------------------------
void RK1:: solveQuadDefectProblem(Problem& problem, Vector<Num>& error,
Vector<Num>& uQ, Vector<Num>& ADiag)
{
int i, j;
Num sum;
PATCH* patch;
MESH* mesh = problem.Mesh();
Vector<Num>& u = problem.u;
const int edNodeM1 = u.high();
int edNode = edNodeM1;
interfaceDLY->setHighOrderNodes(&edNode);
error.resize(edNode);
ADiag.resize(edNode);
uQ.resize(edNode);
FORALL(error,i) error[i] = ADiag[i] = 0.0;
Vector<Num> uPrevOnNewMesh(1);
TransientProblem* trP = 0;
if ((trP = problem.castToTransientProblem()))
{
uPrevOnNewMesh.resize(edNodeM1+1, edNode);
trP->solutionToNewMesh(uPrevOnNewMesh, interfaceDLY);
}
const int low = problem.element->NoOfNodes()+1;
const int high = elementDLY->NoOfNodes();
Matrix<Num> AElem(high, high);
Vector<Num> bElem(high);
Vector<int> globalNodes(high);
Jacobian Jac;
Vector<Real> LNorm(high);
Vector<Real> lagrangeAreaWeights(error.high());
FORALL(lagrangeAreaWeights,i) lagrangeAreaWeights[i] = 0.0;
Matrix<int> FullDiag(high,high), LowDiag(high,high);
FORALL_ROWS(FullDiag,i)
FORALL_COLUMNS(FullDiag,j) FullDiag(i,j) = LowDiag(i,j) = False;
for (i=1; i<=high; ++i) FullDiag(i,i) = True;
for (i=1; i< low; ++i) LowDiag(i,i) = True;
mesh->resetElemIter();
while ((patch = mesh->elemIterAll()))
{
patch->compJinv(Jac);
elementDLY->initAb(AElem,bElem);
problem.assemble(*elementDLY, *patch, Jac, AElem, bElem,
elementDLY->DLYPattern, True);
interfaceDLY->getGlobalNodes(patch, globalNodes);
for (i=low; i<=high; ++i) // sum up diagonal AQQ
{ // and right hand side
if (globalNodes[i] == 0) continue;
ADiag[globalNodes[i]] += AElem(i,i);
error[globalNodes[i]] += bElem[i];
sum = 0.0;
for (j=1; j<low; ++j) // subtract 'coupling terms' from rhs
{
if (globalNodes[j]) sum += AElem(i,j) * u[globalNodes[j]];
}
error[globalNodes[i]] -= sum;
}
problem.assemble(*lagrangeElementDLY, *patch, Jac, AElem, bElem,
&LowDiag, True);
for (i=1; i<low; ++i) // sum up order 1 diagonal terms of Lagrange-Elem.
{
if (globalNodes[i]) ADiag[globalNodes[i]] += AElem(i,i);
}
FORALL(LNorm,i) LNorm[i] = 0.0;
lagrangeElementDLY->assembleLNorm(*patch, LNorm, Jac, &FullDiag);
for (i=1; i<=high; ++i)
if (globalNodes[i]) lagrangeAreaWeights[globalNodes[i]] += LNorm[i];
}
if (trP) // transient problem: contrib. of previous step to rhs
{
const NonLinearity* nonLin = problem.getNonLinearity();
FORALL(uPrevOnNewMesh,i)
error[i] += lagrangeAreaWeights[i] *
nonLin->HValue(real(uPrevOnNewMesh[i]));
}
problem.Ab->Residual(error,u);
interfaceDLY->rhsToNB(error); // residual from hier. to nodal Basis
FORALL(uQ,i) uQ[i] = 0.0;
FORALL(u,i) uQ[i] = u[i];
interfaceDLY->solToNB(uQ); // interpolate solution from lin. to
// quad nodal Basis
nonLinCorrection(error, uQ, ADiag, lagrangeAreaWeights, problem);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void RK1:: nonLinCorrection(Vector<Num>& error, const Vector<Num>& uQ,
const Vector<Num>& ADiag,
const Vector<Real>& lagrangeAreaWeight,
Problem& problem)
{
int i;
Num res, a;
Real uppDefO, lowDefO; // dummies
Bool critical;
const int dim = error.high();
Vector<Real> extUppO(dim), extLowO(dim);
const NonLinearity* nonLin = problem.getNonLinearity();
nonLin->updateObstacles(extUppO, extLowO, *interfaceDLY);
// -- one non-linear Jacobi-step:
for (i=1; i<=dim; ++i)
{
res = error[i];
nonLin->defectCorrection(error[i], res, uQ[i], ADiag[i], a, uppDefO,
lowDefO, critical, extUppO[i], extLowO[i],
lagrangeAreaWeight[i]);
}
}
//-------------------------------------------------------------------------
void RK1:: enforcedPMediaRefinement(Problem& problem)
{
int i;
PATCH* patch;
MESH* mesh = problem.Mesh();
doRefine.resize(mesh->noOfElements());
FORALL(doRefine,i) doRefine[i] = False;
NonLinearity* NonLin = problem.getNonLinearity();
if (NonLin->castToPorousMedia())
{
// -- force elements at the 'frontal boundary' to be refined:
Vector<Num>& u = problem.u;
Vector<int> nodes(problem.element->NoOfNodes());
Num uMin = u[1];
Num uMax = u[1];
FORALL(u,i)
{
if (u[i] < uMin) uMin = u[i];
if (u[i] > uMax) uMax = u[i];
}
Num upperBound = 1e-4; // upperBound=2.525e-5
upperBound *= (uMax - uMin);
int no = 0;
mesh->resetElemIter();
while ((patch = mesh->elemIterAll()))
{
++no;
problem.interface->getGlobalNodes(patch, nodes);
uMin = 2.0*upperBound;
Bool allNodesCritical = True;
FORALL(nodes,i)
if (nodes[i])
{
if (!NonLin->isCritical(nodes[i])) allNodesCritical = False;
if (u[nodes[i]] < uMin) uMin = u[nodes[i]];
}
if (uMin < upperBound) doRefine[no] = True;
if (allNodesCritical) doRefine[no] = False;
}
// -- extend forced refinement pattern:
int noOfNeighbPatches=0;
no = 0;
mesh->resetElemIter();
while ((patch = mesh->elemIterAll()))
{
++no;
if (doRefine[no]) patch->setMark();
else patch->unMark();
noOfNeighbPatches = patch->noOfNeighbours();
}
Vector<PATCH*> neighb(noOfNeighbPatches);
no = 0;
mesh->resetElemIter();
while ((patch = mesh->elemIterAll()))
{
++no;
if (!doRefine[no])
{
patch->getNeighbours(neighb);
FORALL(neighb,i)
{
if (neighb[i])
if (neighb[i]->marked())
{
doRefine[no] = True;
break;
}
}
}
}
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
RKA:: RKA(Element* elementDLY0, Element* elementLag0, Interface* interface0)
: RK1(elementDLY0, elementLag0, interface0)
{ }
//-------------------------------------------------------------------------
Bool RKA:: estimateError(Problem& problem)
{
if (infoErrorEstimator) cout << "\n\tError Est.: RKA";
distributeError(problem);
enforcedPMediaRefinement(problem);
return True;
}
//-------------------------------------------------------------------------
void RKA:: distributeError(Problem& problem)
{
int i;
Vector<Num> error(1), uQ(1);
MLSparseMatrix AQ(problem.Ab->symmetry, problem.SpaceDim());
solveQuadProblem(problem, error, uQ, AQ);
AQ.Mult(uQ,error);
EnergyError = 0.5 * Abs(cdot(uQ,error));
interfaceDLY->solToHB(error);
for (i=problem.u.high()+1; i<=error.high(); ++i)
{
error[i] = error[i] * error[i] * AQ.Diag(i); // error on edge
}
distributeEdgeErrors(problem, error);
//cout << "\n-- EnergyErr " << EnergyError << " EdgeErr " << EdgeError
// << " d=" << 100*Abs(EdgeError-EnergyError)/EnergyError << " %\n";
}
//-------------------------------------------------------------------------
void RKA:: solveQuadProblem(Problem& problem, Vector<Num>& error,
Vector<Num>& uQ, MLSparseMatrix& AQ)
{
int i;
PATCH* patch;
MESH* mesh = problem.Mesh();
const Vector<Num>& u = problem.u;
const int edNodeM1 = u.high();
int edNode = edNodeM1;
interfaceDLY->setHighOrderNodes(&edNode);
uQ.resize(edNode);
error.resize(edNode);
FORALL(error,i) error[i] = 0.0;
Vector<Num> uPrevOnNewMesh(1);
TransientProblem* trP = 0;
if ((trP = problem.castToTransientProblem()))
{
uPrevOnNewMesh.resize(1,edNode);
trP->solutionToNewMesh(uPrevOnNewMesh, interfaceDLY);
}
// -- full assembling procedure for new lin. system:
AQ.extend(edNode);
LinSystem AbQ(problem.spaceDim, AQ.symmetry, problem.precond,
problem.dirichletBCsDLY);
AbQ.update(&AQ);
AbQ.reset();
const int high = lagrangeElementDLY->NoOfNodes();
Matrix<Num> AElem(high,high);
Vector<Num> bElem(high);
Vector<int> globalNodes(high);
Jacobian Jac;
Vector<Real> LNorm(high);
Vector<Real> lagrangeAreaWeights(error.high());
FORALL(lagrangeAreaWeights,i) lagrangeAreaWeights[i] = 0.0;
mesh->resetElemIter();
while ((patch = mesh->elemIterAll()))
{
patch->compJinv(Jac);
lagrangeElementDLY->initAb(AElem,bElem);
problem.assemble(*lagrangeElementDLY,*patch,Jac,AElem,bElem,0,True);
// 0: use default symm. pattern
interfaceDLY->getGlobalNodes(patch, globalNodes);
AbQ.store(AElem, bElem, globalNodes);
FORALL(LNorm,i) LNorm[i] = 0.0;
lagrangeElementDLY->assembleLNorm(*patch,LNorm,Jac,0);
for (i=1; i<=high; ++i)
if (globalNodes[i]) lagrangeAreaWeights[globalNodes[i]] += LNorm[i];
}
Real time = 0.0;
if (trP) // transient problem: contrib. of previous step to rhs
{
const NonLinearity* nonLin = problem.getNonLinearity();
FORALL(uPrevOnNewMesh,i)
AbQ.b[i] += lagrangeAreaWeights[i] * nonLin->HValue(uPrevOnNewMesh[i]);
time = trP->time;
}
interfaceDLY->updateDirichletBCs(time);
AbQ.setDirichletBCs(*problem.dirichletBCsDLY);
// -- end of assembling proc.
FORALL(uQ,i) uQ[i] = 0.0;
FORALL(u,i) uQ[i] = u[i];
interfaceDLY->solToNB(uQ); // interpolate solution from lin. to
// quad nodal Basis
AbQ.Residual(error,uQ);
nonLinCorrection(error, uQ, AQ, lagrangeAreaWeights, problem);
}
//-------------------------------------------------------------------------
void RKA:: nonLinCorrection(Vector<Num>& error, const Vector<Num>& uQ,
MLSparseMatrix& AQ,
const Vector<Real>& lagrangeAreaWeight,
Problem& problem)
{
int i;
Num eNew, res, a;
Real uppDefO, lowDefO;
Bool critical;
const int quadDim = error.high();
Vector<Real> extLowO(quadDim), extUppO(quadDim);
const NonLinearity* nonLin = problem.getNonLinearity();
nonLin->updateObstacles(extUppO, extLowO, *interfaceDLY);
if (Cmd.isSet("RKSmooth","GS")) // non-linear Gauss-Seidel steps:
{
AQ.symmetricExpansion();
Vector<Num> r(quadDim);
FORALL(r,i) { r[i] = error[i]; error[i] = 0.0; }
for (i=1; i<=quadDim; ++i)
{
if (problem.dirichletBCsDLY->isSet(i)) continue;
res = AQ.MultRow(i,error);
res = r[i] - res;
nonLin->defectCorrection(error[i], res, uQ[i], AQ.Diag(i), a,
uppDefO, lowDefO, critical, extUppO[i],
extLowO[i], lagrangeAreaWeight[i]);
}
for (i=quadDim; i>=1; --i)
{
if (problem.dirichletBCsDLY->isSet(i)) continue;
res = AQ.MultRow(i,error);
res = r[i] - res;
nonLin->defectCorrection(eNew, res, uQ[i], AQ.Diag(i), a,
uppDefO, lowDefO, critical, extUppO[i],
extLowO[i], lagrangeAreaWeight[i]);
error[i] += eNew;
}
}
else // non-linear Jacobi-step:
{
for (i=1; i<=quadDim; ++i)
{
if (problem.dirichletBCsDLY->isSet(i)) continue;
res = error[i];
nonLin->defectCorrection(error[i], res, uQ[i], AQ.Diag(i), a,
uppDefO, lowDefO, critical, extUppO[i],
extLowO[i], lagrangeAreaWeight[i]);
}
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
RK2:: RK2(Element* elementDLY0, Element* elementLag0, Interface* interface0)
: RKA(elementDLY0, elementLag0, interface0)
{ }
//-------------------------------------------------------------------------
Bool RK2:: estimateError(Problem& problem)
{
if (infoErrorEstimator) cout << "\n\tError Est.: RK2";
if (problem.precond->mode() != multiLevel)
{
cout << "\n*** RKC:: Preconditioner must be of multi-level-type:\n"
<< "NonLinearMLSG or TrcNonLinearMLSG\n";
cout.flush(); abort();
}
distributeError(problem);
enforcedPMediaRefinement(problem);
return True;
}
//-------------------------------------------------------------------------
void RK2:: nonLinCorrection(Vector<Num>& error, const Vector<Num>& uQ,
MLSparseMatrix& AQ,
const Vector<Real>& lagrangeAreaWeight,
Problem& problem)
{
int i;
Num eNew, res, a;
Real uppDefO, lowDefO;
Bool critical;
const int quadDim = error.high();
Vector<Real> extLowO(quadDim), extUppO(quadDim);
const NonLinearity* nonLin = problem.getNonLinearity();
nonLin->updateObstacles(extUppO, extLowO, *interfaceDLY);
const int linDim = nonLin->lowObstacle.high();
Vector<Num> r(quadDim), Diag0(quadDim), aux(quadDim);
FORALL(r,i)
{
r[i] = error[i];
error[i] = 0.0;
Diag0[i] = AQ.Diag(i);
}
NonLinearMLGS* precond = (NonLinearMLGS*) problem.precond;
Vector<Real>& uppDefObstacle = *precond->uppDefObstacle.Top();
Vector<Real>& lowDefObstacle = *precond->lowDefObstacle.Top();
FORALL(uppDefObstacle,i)
{
uppDefObstacle[i] = lowDefObstacle[i] = 0.0;
precond->setUnCritical(i);
}
if (Cmd.isSet("RKSmooth","GS")) // non-linear Gauss-Seidel steps:
{
AQ.symmetricExpansion();
for (i=1; i<=quadDim; ++i)
{
if (problem.dirichletBCsDLY->isSet(i)) continue;
res = AQ.MultRow(i,error);
res = r[i] - res;
nonLin->defectCorrection(error[i], res, uQ[i], AQ.Diag(i), a,
uppDefO, lowDefO, critical, extUppO[i],
extLowO[i], lagrangeAreaWeight[i]);
}
for (i=quadDim; i>=1; --i)
{
if (problem.dirichletBCsDLY->isSet(i)) continue;
res = AQ.MultRow(i,error);
res = r[i] - res;
nonLin->defectCorrection(eNew, res, uQ[i], AQ.Diag(i), a,
uppDefO, lowDefO, critical, extUppO[i],
extLowO[i], lagrangeAreaWeight[i]);
r[i] -= res;
error[i] += eNew;
AQ.Diag(i) += a;
if (i <= linDim)
{
uppDefObstacle[i] = uppDefO - real(error[i]);
lowDefObstacle[i] = lowDefO - real(error[i]);
if (critical) precond->setCritical(i);
else precond->setUnCritical(i);
}
}
}
else // non-linear Jacobi-step:
{
for (i=1; i<=quadDim; ++i)
{
if (problem.dirichletBCsDLY->isSet(i)) continue;
res = r[i];
nonLin->defectCorrection(error[i], res, uQ[i], AQ.Diag(i), a,
uppDefO, lowDefO, critical, extUppO[i],
extLowO[i], lagrangeAreaWeight[i]);
AQ.Diag(i) += a;
r[i] -= res;
if (i <= linDim)
{
uppDefObstacle[i] = uppDefO - real(error[i]);
lowDefObstacle[i] = lowDefO - real(error[i]);
if (critical) precond->setCritical(i);
else precond->setUnCritical(i);
}
}
}
// -- restrict Matrix and residual to hierarchical linear basis:
Generation gen(quadDim);
interfaceDLY->setHBGeneration(gen); // for Galerkin-Restriction
MLMatrix* ALin = precond->Al.Top()->castToMLMatrix();
ALin->reset();
ALin->markNodes();
ALin->GalerkinRestriction(AQ, gen);
AQ.Mult(aux,error);
FORALL(r,i) aux[i] = r[i] - aux[i]; // the new residual
interfaceDLY->rhsToHB(aux);
const int maxLevel = precond->maxLevel;
Vector<Num> rLin(linDim);
FORALL(rLin,i)
{
if (precond->isCritical(i,maxLevel)) rLin[i] = 0.0;
else if (problem.dirichletBCs->isSet(i)) rLin[i] = 0.0;
else rLin[i] = aux[i];
}
// initialize the multi-level preconditioner and call one V-Cycle:
const int nPreSmooth0 = precond->nPreSmooth;
const int nPostSmooth0 = precond->nPostSmooth;
precond->nPreSmooth = 1;
precond->nPostSmooth = 1;
Vector<Num> b(1); // a dummy
precond->initialize(ALin, problem.u, b);
Vector<Num>& eMG = precond->eMG;
precond->el.push(&eMG);
precond->MGCycle(eMG, *ALin, rLin, precond->maxLevel);
precond->el.pop();
FORALL(aux,i) aux[i] = 0.0;
FORALL(eMG,i)
{
if (!precond->isCritical(i, maxLevel) &&
!problem.dirichletBCs->isSet(i)) aux[i] = eMG[i];
}
interfaceDLY->solToNB(aux);
FORALL(error,i) error[i] += aux[i];
// -- restore preconditioner settings:
precond->nPreSmooth = nPreSmooth0;
precond->nPostSmooth = nPostSmooth0;
FORALL(precond->Al,i) precond->Al[i]->removeSymmetricExpansion();
FORALL(Diag0,i) AQ.Diag(i) = Diag0[i]; // restore diagonal
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
syntax highlighted by Code2HTML, v. 0.9.1