/*
$Id: problem.cc,v 1.4 1996/10/11 09:10:05 roitzsch Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "problem.h"
#include "physics.h"
#include "linsystem.h"
#include "precondsg.h"
#include "precondmg.h"
#include "triang.h"
#include "elements.h"
#include "int.h"
#include "dirichlet.h"
#include "materials.h"
#include "adapt.h"
#include "nodeco.h"
#include "cmdpars.h"
extern CmdPars Cmd;
extern ostream* infoFile;
//-------------------------------------------------------------------------
Problem:: Problem (char* fileName0, symmetryType symmetry0, int spaceDim0,
int nComp0) :
symmetry(symmetry0), spaceDim(spaceDim0), nComp(nComp0),
Fct(0,0), Energy(0,0), Error(0,0)
{
fileName = 0;
material = 0;
dirichletBCs = dirichletBCsDLY = 0;
precond = 0;
Ab = 0;
element = 0;
interface = 0;
errorEstimator= 0;
statistic = new Statistic();
fileName = new char[strlen(fileName0)+10];
strcpy(fileName, fileName0);
// inquire the problem-specific parameters:
biSection = Cmd.isTrue("biSection");
infoSolution = True; Cmd.get("infoSolution", &infoSolution);
plotSolution = False; Cmd.get("plotSolution", &plotSolution);
postScript = False; Cmd.get("postScript", &postScript);
writeSolution = False; Cmd.get("writeSolution",&writeSolution);
compareSolution = False; Cmd.get("compareSolution",&compareSolution);
timeAssembler = 0; Cmd.get("timeAssembler", &timeAssembler);
accTimeAssembler = 0; Cmd.get("accTimeAssembler", &accTimeAssembler);
maxRefSteps = 50; Cmd.get("maxRefSteps", &maxRefSteps);
minRefSteps = 0; Cmd.get("minRefSteps", &minRefSteps);
maxRefDepth = 50; Cmd.get("maxRefDepth", &maxRefDepth);
globConvTest = estimRelError;
}
//-------------------------------------------------------------------------
Problem:: ~Problem()
{
delete fileName;
delete precond;
delete material;
delete dirichletBCs; delete dirichletBCsDLY;
delete Ab;
delete element;
delete interface;
delete errorEstimator;
}
//-------------------------------------------------------------------------
// NonLinearity* Problem:: getNonLinearity() const { return 0; }
//-------------------------------------------------------------------------
void Problem:: newLinSystem()
{
if (precond == 0) missingObject("newLinSystem","precond");
if (dirichletBCs == 0) missingObject("newLinSystem","dirichletBCs");
delete Ab;
Ab = new LinSystem(spaceDim, symmetry, precond, dirichletBCs);
}
//-------------------------------------------------------------------------
void Problem:: newPreconditioner()
{
delete precond;
if (Cmd.isSet("precond","jacobi")) precond = new Jacobi();
else if (Cmd.isSet("precond","sgs")) precond = new SGS();
else if (Cmd.isSet("precond","ssor")) precond = new SSOR();
else if (Cmd.isSet("precond","TRssor")) precond = new TrSSOR();
else if (Cmd.isSet("precond","ILU")) precond = new ILU();
else if (Cmd.isSet("precond","MGsgs")) precond = new MGSGS();
else if (Cmd.isSet("precond","MGssor")) precond = new MGSSOR();
else if (Cmd.isSet("precond","MGjac")) precond = new MGJacobi();
else if (Cmd.isSet("precond","MGCG")) precond = new MGCG();
else if (Cmd.isSet("precond","HB")) precond = new HB();
else if (Cmd.isSet("precond","AMGjac")) precond = new AMGJacobi();
else if (Cmd.isSet("precond","AMGsgs")) precond = new AMGSGS();
else if (Cmd.isSet("precond","MLJacobi")) precond = new MLJacobi();
else if (Cmd.isSet("precond","MLsgs")) precond = new MLSGS();
else if (Cmd.isSet("precond","MLssor")) precond = new MLSSOR();
else if (Cmd.isSet("precond","AMLJacobi")) precond = new AMLJacobi();
else if (Cmd.isSet("precond","BPX")) precond = new AMLJacobi();
else if (Cmd.isSet("precond","AMLsgs")) precond = new AMLSGS();
else if (Cmd.isSet("precond","none")) precond = new DummyPreconditioner();
else MissingParameter("Preconditioner");
}
//-------------------------------------------------------------------------
void Problem:: notImplemented(const char* s) const
{
cout << "\n*** class Problem: " << s << " not implemented\n";
cout.flush(); abort();
}
//-------------------------------------------------------------------------
void Problem:: missingObject(const char* function, const char* object) const
{
cout << "\n*** " << function << ": " << object << " missing\n";
cout.flush(); abort();
}
//-------------------------------------------------------------------------
void Problem:: missingMaterialTerm(const char* name, int classVal) const
{
cout << "\n*** assemble: " << name << " value for material class "
<< classVal << " missing\n";
cout.flush(); abort();
}
//-------------------------------------------------------------------------
Bool Problem:: globalConvergenceTest(Real globalPrecision, int step)
{
switch (globConvTest)
{
case estimRelError: return estimRelErrorConvTest(globalPrecision, step);
case estimAbsError: return estimAbsErrorConvTest(globalPrecision, step);
//case extrapolError: return extrapolatedErrorConvTest(globalPrecision,
// step);
default: cout << "\n*** Problem: no valid convergence Test specified\n";
cout.flush(); abort(); return 0;
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void Problem:: contractVector(Vector<Num>& newVec, Vector<int>& oldPos,
Vector<Num>& oldVec)
{
int i;
if (newVec.high() < oldVec.high())
{ cout << "\nProblem.contract: inconsistent vector sizes \n"; cout.flush(); abort(); }
FORALL(newVec,i) { if (oldPos[i]) oldVec[oldPos[i]] = newVec[i]; }
}
void Problem:: expandVector(Vector<Num>& oldVec, Vector<int>& oldPos,
Vector<Num>& newVec, Num fill)
{
int i;
if (newVec.high() < oldVec.high())
{ cout << "\nProblem.expand: inconsistent vector sizes \n"; cout.flush(); abort(); }
FORALL(newVec,i)
{
if (oldPos[i]) newVec[i] = oldVec[oldPos[i]];
else newVec[i] = fill;
}
}
//-------------------------------------------------------------------------
Bool Problem:: estimRelErrorConvTest(Real globalPrecision, int step)
{
/*
if (equal(Energy[step],0.0)) return False;
Real relError = Error[step]/Energy[step];
if (Cmd.isTrue("writeCpuTime")) *infoFile << "\t " << relError << "\n";
if (relError < globalPrecision) return True;
else return False;
*/
Bool ConvOk = (Error[step] < globalPrecision * Energy[step])||
(Error[step]==0);
if (Cmd.isTrue("writeCpuTime"))
{ if (equal(Energy[step],0.0))
{
if (ConvOk) *infoFile << "\t < " << globalPrecision << "\n";
else *infoFile << "\t " << machMax(Real(0.0)) << "\n";
}
else
{
Real relError = Error[step]/Energy[step];
*infoFile << "\t " << relError << "\n";
}
}
return ConvOk;
}
//-------------------------------------------------------------------------
Bool Problem:: estimAbsErrorConvTest(Real globalPrecision, int step)
{
if ( (Error[step] < globalPrecision) || (Error[step]==0) )
return True;
else return False;
}
//-------------------------------------------------------------------------
Bool Problem:: extrapolatedErrorConvTest(Real globalPrecision, int step)
{
Real ratio, expo, extrapolatedError = machMax(Real(0));
static int dimM1;
Bool status = False;
if (step > 0)
{
extrapolatedError = Abs(Energy[step]-Energy[step-1]);
expo = 2.0/Real(spaceDim);
ratio = Real(dimM1)/Real(u.size());
ratio = pow(ratio,expo);
extrapolatedError *= ratio/(1.0-ratio);
if ( (extrapolatedError < globalPrecision*Energy[step]) ||
(extrapolatedError ==0) ) status = True;
}
if (Cmd.isTrue("writeCpuTime"))
*infoFile << "\t " << quot(extrapolatedError,Energy[step]) << "\n";
dimM1 = u.size();
return status;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void Problem:: SolutionInfo(int step)
{
static int i, dimM1;
static Timer accTimer;
static Real accTime = 0.0;
Real error, extrapolatedError=0.0, ratio=0.0, expo;
if (step > 0)
{
extrapolatedError =
Abs(quot(Energy[step]-Energy[step-1],Energy[step]));
expo = 2.0/Real(spaceDim);
ratio = Real(dimM1)/Real(u.size());
ratio = pow(ratio,expo);
extrapolatedError *= 100.0 * ratio/(1.0-ratio);
ratio = quot(Error[step],Error[step-1]);
}
if (!equal(Energy[step],0.0)) error = 100.0*Error.Top()/Energy.Top();
else error = Error.Top();
cout << Form("\n %7s %5s %6s %11s","RefStep", "level", "nodes", "Energy");
if (!equal(Energy[step],0.0)) cout << Form(" %8s","Error(%)");
else cout << Form(" %8s","Error");
if (step > 0) cout << Form(" %7s %8s", "ratio", "dE(%)");
cout << Form(" %9s", "cpu(sec)");
// cout << Form(" %9s", "mem(MB)");
cout << Form("\n %7d %5d %6d %11.6g", step, interface->MaxDepth(),
interface->Dim(), Energy.Top());
cout << Form(" %8.3g", error);
if (step > 0)
{
cout << Form(" %7.2g %8.3g", ratio,
quot(100.0*(Energy[step]-Energy[step-1]),Energy[step]));
}
accTime = accTimer.cpu(False);
cout << Form(" %9.2f", accTime);
//int space = mesh->MemSpace(); // + Ab->MemSpace() + precond->MemSpace();
//cout << Form(" %9.6f", float(space)/1e6);
/*
cout << "\n Step " << step << ": E=" << Form("%1.5g", Energy[step]);
if (!equal(Energy[step],0.0))
cout << " Err=" << Form("%1.3g %%", 100.*Error[step]/Energy[step]);
cout << Form(" (ext=%1.3g)", 100*extrapolatedError);
cout << Form(" ratio(Err)=%1.2f",(Error[step] / Error[step-1]));
cout << Form(" dE(%1d)=%1.3g", step-1,
100*(Energy[step]-Energy[step-1])/Energy[step]) << " %";
*/
//---------------------------------------------------------------
cout << "\n "; for (i=1; i<=76; ++i) cout << '-'; cout << "\n";
//---------------------------------------------------------------
if (Cmd.isTrue("EnergyNorm")) cout << "\n E-Norm: " << EnergyNorm(u) <<"\n";
if (Cmd.isTrue("L2Norm")) cout << "\n L2-Norm: " << L2Norm(u) << "\n";
//---------------------------------------------------------------
dimM1 = u.size();
while ( Pause() == picture )
interface->post (u, step, fileName);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void Problem:: compEnergy(Vector<Num>& sol, Real* energy, Num* fct)
{
Bool print=True;
Ab->compEnergy(sol, energy, fct, print);
}
//-------------------------------------------------------------------------
Real Problem:: L2Norm(const Vector<Num>& x) const
{
int dim;
Real L2 = 0.0;
PATCH* patch;
dim = element->NoOfNodes();
Matrix<Num> AElem(dim,dim);
Vector<int> globalNodes(dim);
Jacobian Jac;
Mesh()->resetElemIter();
while ((patch = Mesh()->elemIterAll()))
{
patch->compJinv(Jac);
element->initAb(AElem);
element->assembleL2Norm(*patch, AElem, Jac, 0);
interface->getGlobalNodes(patch, globalNodes);
L2 += energyProd(AElem, x, globalNodes);
}
// L2 = sqrt(2.0*L2); // annihilate factor 1/2 in energy
return L2;
}
//-------------------------------------------------------------------------
Real Problem:: L2MassNorm(const Vector<Num>& x) const
{
int dim;
Real L2 = 0.0;
PATCH* patch;
dim = element->NoOfNodes();
Matrix<Num> AElem(dim,dim);
Vector<int> globalNodes(dim);
Jacobian Jac;
Mesh()->resetElemIter();
while ((patch = Mesh()->elemIterAll()))
{
patch->compJinv(Jac);
element->initAb(AElem);
element->assembleMass(*patch, AElem, Jac, 0);
interface->getGlobalNodes(patch, globalNodes);
L2 += energyProd(AElem, x, globalNodes);
}
// L2 = sqrt(2.0*L2); // annihilate factor 1/2 in energy
return L2;
}
//-------------------------------------------------------------------------
Real Problem:: EllipticNorm(const Vector<Num>& x) const
{
int dim;
Real E = 0.0;
PATCH* patch;
dim = element->NoOfNodes();
Matrix<Num> AElem(dim,dim);
Vector<Num> bElem(dim);
Vector<int> globalNodes(dim);
Jacobian Jac;
Mesh()->resetElemIter();
while ((patch = Mesh()->elemIterAll()))
{
patch->compJinv(Jac);
element->initAb(AElem,bElem);
element->assembleEllip(*patch, AElem, Jac, 0);
interface->getGlobalNodes(patch, globalNodes);
E += energyProd(AElem, x, globalNodes);
}
return E;
}
//-------------------------------------------------------------------------
Real Problem:: EnergyNorm(const Vector<Num>& x)
{
int dim;
Real E = 0.0;
PATCH* patch;
dim = element->NoOfNodes();
Matrix<Num> AElem(dim,dim);
Vector<Num> bElem(dim);
Vector<int> globalNodes(dim);
Jacobian Jac;
Mesh()->resetElemIter();
while ((patch = Mesh()->elemIterAll()))
{
patch->compJinv(Jac);
element->initAb(AElem,bElem);
assemble(*element, *patch, Jac, AElem, bElem, 0, True);
interface->getGlobalNodes(patch, globalNodes);
E += energyProd(AElem, x, globalNodes);
}
return E;
}
//-------------------------------------------------------------------------
Real Problem:: energyProd(Matrix<Num>& AElem, const Vector<Num>& sol,
const Vector<int>& globalNodes) const
{
int i, j, row, col;
Num sum;
Num E = 0.0;
const int dim = globalNodes.high();
if (symmetry == sym) // upper triangle
{
for (i=1; i<=dim; ++i)
for (j=1; j<i; ++j) AElem(j,i) = AElem(i,j);
}
for (i=1; i<=dim; ++i)
{
if ((row=globalNodes[i]))
{
sum = 0.0;
for (j=1; j<=dim; ++j)
if ((col=globalNodes[j])) sum += AElem(i,j) * sol[col];
E += conj(sol[row]) * sum;
}
}
return 0.5*Abs(E);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Bool Problem:: solve(Real relGlobalPrecision)
{
Bool absPrecision=False;
if (Cmd.isTrue("absPrecision")) absPrecision=True;
return staticSolution(relGlobalPrecision, absPrecision);
}
//-------------------------------------------------------------------------
void Problem:: compare(Vector<Num>& v, Real time)
{
if (material->trueSolKnown())
{
Real maxError=0.0, error=0.0;
Num trueSol=0.0;
int i;
const int dim = interface->Dim();
// const spaceDim = interface->SpaceDim();
Vector<Real> x(spaceDim);
NodeCoordinates nc(spaceDim, 1, dim);
interface->getNodeCoordinates(nc);
FORALL(v,i)
{
nc.getCoordinate(i,x);
trueSol = material->trueSolInPoint(x,time);
error = Abs(v[i] - trueSol);
if (error > maxError) maxError = error;
}
cout << "\n\tCompare with true solution: Max. error at mesh points="
<< maxError << "\n";
}
else cout << "\n Can not compare, since true solution is not known\n";
}
//-------------------------------------------------------------------------
Bool Problem:: staticSolution(Real globalPrecision, Bool absPrecision, Real time)
{
int i, step;
Num fct;
Real error=0.0, energy=1.0, relPrecision;
if (absPrecision) globConvTest = estimAbsError;
else globConvTest = estimRelError;
Energy.resize(0,DefaultStackSize);
Error.resize (0,DefaultStackSize);
Fct.resize (0,DefaultStackSize);
newPreconditioner();
newLinSystem();
newMesh();
newInterface();
newErrorEstimator();
interface->updateLevel0(u);
for (step=0; True; ++step)
{
if (Cmd.isTrue("uReset")) FORALL(u,i) u[i] = 0.0;
assembleGlobal(time);
// cout << "\n\tNo of Dirichlet points = "
// << dirichletBCs->noOfConstraints() << "\n";
interface->updatePrecond();
if (absPrecision && !(energy==0)) relPrecision = globalPrecision/energy;
else relPrecision = globalPrecision;
if (statistic->active)
{
statistic->ZD_IntWrite(statistic->idRefStep,step);
statistic->ZD_IntWrite(statistic->idDepth,interface->MaxDepth());
statistic->ZD_IntWrite(statistic->idNodes,interface->Dim());
}
Ab->statistic = statistic;
Ab->solve(u, energy, error, relPrecision, step);
compEnergy(u, &energy, &fct);
Energy.push(energy);
Fct.push(fct);
if (plotSolution) interface->plot (u, step, fileName);
if (postScript) interface->post (u, step, fileName);
if (writeSolution) interface->write(u, step, fileName);
if (compareSolution) compare(u,time);
if (absPrecision)
errorEstimator->adapt(*this, &error, globalPrecision);
else
errorEstimator->adapt(*this, &error, relPrecision*energy);
Error.push(error);
if (statistic->active)
{
statistic->ZD_RealWrite(statistic->idDiscErr,error);
}
if (infoSolution) SolutionInfo(step);
if (step >= minRefSteps)
{
if (globalConvergenceTest(globalPrecision, step)) return True;
if (step == maxRefSteps)
{
cout << "\n* maxRefStep " << maxRefSteps << " reached\n";
return False;
}
if (interface->MaxDepth() > maxRefDepth)
{
cout << "\n* maxRefDepth " << maxRefDepth << " exceeded\n";
return False;
}
}
interface->refine(u);
}
}
//-------------------------------------------------------------------------
void Problem:: assembleGlobal(Real time)
{
static Real accTime = 0.0;
Timer timer, accTimer;
PATCH* patch;
const int dim = element->NoOfNodes();
Matrix<Num> AElem(dim,dim);
Vector<Num> bElem(dim);
Vector<int> nodes(dim);
Jacobian Jac;
Ab->reset();
Mesh()->resetElemIter();
while ((patch = Mesh()->elemIterAll()))
{
patch->compJinv(Jac);
element->initAb(AElem,bElem);
assemble(*element, *patch, Jac, AElem, bElem, 0);
interface->getGlobalNodes(patch, nodes);
Ab->store(AElem, bElem, nodes);
}
interface->updateDirichletBCs(time);
Ab->setDirichletBCs(*dirichletBCs);
// -- infos:
if (timeAssembler) { cout << "\n\tGlobal Assembling: "; timer.cpu(); }
if (accTimeAssembler)
{
accTime += accTimer.cpu(False);
cout << Form("\tAccumulated time: %1.2f sec.\n", accTime);
}
}
//-------------------------------------------------------------------------
void Problem:: resolveMesh(Real lambda0)
{
Timer timer;
int nElem;
PATCH* patch;
Real lambda, h, eps, mu;
lambda0 *= 0.999;
while(1)
{
nElem = 0;
Mesh()->resetElemIter();
while ((patch = Mesh()->elemIterAll()))
{
mu = 1/material->E(patch->Class());
eps = material->M(patch->Class());
h = patch->hMax();
lambda = lambda0 / sqrt(eps*mu);
if (h > lambda) { patch->setMark(); ++nElem; }
else patch->unMark();
}
if (nElem == 0) break;
Mesh()->Refine();
}
if (spaceDim != 3) Mesh()->Flatten(); // !!!
if (timeAssembler) { cout << "\n\tResolve Mesh: "; timer.cpu(); }
}
syntax highlighted by Code2HTML, v. 0.9.1