/*
$Id: problemtr.cc,v 1.5 1996/11/20 10:00:10 roitzsch Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "problemtr.h"
#include "nodeco.h"
#include "physics.h"
#include "triang.h"
#include "elements.h"
#include "dirichlettr.h"
#include "materialstr.h"
#include "adapt.h"
#include "linsystem.h"
#include "sysmatml.h"
#include "sysmatbl.h"
#include "int.h"
#include "cmdpars.h"
extern CmdPars Cmd;
//-------------------------------------------------------------------------
TransientProblem:: TransientProblem() : uPrev(1), uPrevOnNewMesh(1)
{
transientMode = True;
lumpedMass = True; Cmd.get("lumpedMass", &lumpedMass);
plotTimeStep = False; Cmd.get("plotTimeStep", &plotTimeStep);
postTimeStep = False; Cmd.get("postTimeStep", &postTimeStep);
writeTimeStep = False; Cmd.get("writeTimeStep",&writeTimeStep);
maxTimeSteps=1000; Cmd.get("maxTimeSteps", &maxTimeSteps);
startTime = 0.0; Cmd.get("startTime", &startTime);
if (!Cmd.get("endTime",&endTime)) MissingParameter("endTime");
theta = 1.0; Cmd.get("theta",&theta);
if (theta < 0.5 || theta > 1)
{
cout << "\n*** Wrong value for theta given (" << theta
<< ") -- must be between 0.5 and 1.0\n";
cout.flush(); abort();
}
if (equal(theta,1.0)) implicitEuler = True;
else implicitEuler = False;
midPointRule = False; Cmd.get("midPointRule", &midPointRule);
if (midPointRule) { theta = 0.5; implicitEuler = False; }
if ((constTimeStep = Cmd.isTrue("constTimeStep")))
{
if (!Cmd.get("tau", &tau)) missingTau();
if (tau >= endTime) endTime = tau;
}
if (!constTimeStep)
{
tau = 0.1*(endTime-startTime); Cmd.get("Tau",&tau); // init. guess
staticFirstStep = Cmd.isTrue("staticFirstStep");
if ((fixedFirstStep = Cmd.isTrue("fixedFirstStep")))
{
if (!Cmd.get("tau", &tau)) missingTau();
if (tau >= endTime) tau = 0.99*endTime;
}
maxTau = 0.1*(endTime-startTime); Cmd.get("maxTau", &maxTau);
rhoTime = 0.5; Cmd.get("rhoTime", &rhoTime); // 0.18
rhoSpace = 1.0-rhoTime-0.1; Cmd.get("rhoSpace", &rhoSpace);
fSpace = 1.0; Cmd.get("fSpace", &fSpace);
rhoSpace *= fSpace;
maxReductions = 25; Cmd.get("maxReductions", &maxReductions);
dynamicScaling = Cmd.isSet("scaling","dynamic");
}
}
//-------------------------------------------------------------------------
void TransientProblem:: missingTau()
{
cout << "\n*** time step tau not specified\n"; cout.flush(); abort();
}
//-------------------------------------------------------------------------
void TransientProblem:: notImplemented(const char* name) const
{
cout << "\n*** class TransientProblem: function " << name
<< " not implemented\n";
cout.flush(); abort();
}
//-------------------------------------------------------------------------
void TransientProblem:: assembleGlobalDefect(MLMatrix& /*AD*/,
Vector<Num>& /*MD*/, Vector<Num>& /*bD*/)
{
notImplemented("assembleGlobalDefect");
}
//-------------------------------------------------------------------------
Bool TransientProblem:: solve(Real relGlobalPrecision)
{
if (constTimeStep) return fixedTimeSteps (relGlobalPrecision);
else return adaptiveTimeSteps(relGlobalPrecision);
}
//-------------------------------------------------------------------------
Bool TransientProblem:: fixedTimeSteps(Real relGlobalPrecision)
{
Bool lastStep = False, absPrecision = False;
int i, timeStep;
time = startTime + tau;
for (timeStep=1; timeStep<=maxTimeSteps; ++timeStep)
{
staticSolution(relGlobalPrecision, absPrecision, time);
TransientSolutionInfo(timeStep, time, tau);
if (lastStep) break;
else if (time+tau >= endTime) { tau = endTime-time; lastStep = True; }
if (equal(endTime+tau, endTime)) break;
time += tau;
uPrev.resize(u.high());
FORALL(u,i) uPrev[i] = u[i];
shiftMesh();
}
if (timeStep > maxTimeSteps)
{
cout << "\n** Time step limit reached: " << timeStep << "\n";
return False;
}
else return True;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Bool TransientProblem:: adaptiveTimeSteps(Real globalPrecision)
{
int i, timeStep, reductions;
Real statPrecision, newTau;
Real UNorm, timeError, newUNorm;
Bool lastStep = False;
const Bool absPrecision = True;
const Real TOL = sqrt(globalPrecision);
time = startTime + tau;
if (staticFirstStep) staticInitialSolution(&UNorm, TOL);
else if (fixedFirstStep) initialUNorm(&UNorm);
else adaptiveInitialSolution(&UNorm, TOL);
newUNorm = UNorm;
for (timeStep=1; timeStep<=maxTimeSteps; ++timeStep)
{
for (reductions=0; True; ++reductions)
{
if (dynamicScaling) UNorm = newUNorm;
else if (timeStep == 1) UNorm = Max(UNorm,newUNorm);
statPrecision = sqr(rhoSpace*TOL*UNorm);
staticSolution(statPrecision, absPrecision, time);
compTimeError(&timeError, &newUNorm, TOL);
if (timeConvergence(TOL, timeError, UNorm, &newTau, timeStep,
reductions)) break;
if (fixedFirstStep && timeStep==1) break;
time -= tau;
tau = newTau;
time += tau;
}
tau = newTau;
TransientSolutionInfo(timeStep, time, tau);
if (lastStep) break;
else if (time+tau >= endTime) { tau = endTime-time; lastStep = True; }
if (equal(endTime+tau, endTime)) break;
time += tau;
uPrev.resize(u.high());
FORALL(u,i) uPrev[i] = u[i]; // save solution for next step
shiftMesh();
}
if (timeStep > maxTimeSteps)
{
cout << "\n** Time step limit reached: " << timeStep << "\n";
return False;
}
else return True;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
// -- one coarse initial solution step to compute UNorm for scaling:
void TransientProblem:: initialUNorm(Real* UNorm)
{
staticSolution(0.05, False, time); // 5 % are sufficient
*UNorm = compUNorm(u);
cout << Form("\n Initial UNorm = %1.3g\n", *UNorm);
TransientSolutionInfo(0, time, tau);
}
//-------------------------------------------------------------------------
// -- try to find a good guess for the first time step by successively
// reducing the tolerance
void TransientProblem:: adaptiveInitialSolution(Real* UNorm, Real TOL)
{
Real tol, initTol, statPrecision, timeError, newTau;
Bool absPrecision = True;
int reductions;
initialUNorm(UNorm);
Real newUNorm = *UNorm;
if (TOL >= 0.1) initTol = 2.1*TOL;
else initTol = 0.1;
cout << Form("\n Search for initial time step")
<< Form(" - start with TOL = %1.3g, UNorm = %1.3g\n",initTol,*UNorm);
for (tol=initTol; tol >= 2.001*TOL; tol=0.5*tol)
{
for (reductions=0; True; ++reductions)
{
*UNorm = Max(*UNorm,newUNorm);
statPrecision = sqr(rhoSpace*tol*(*UNorm));
staticSolution(statPrecision, absPrecision, time);
compTimeError(&timeError, &newUNorm, tol);
if (timeConvergence(tol, timeError, *UNorm, &newTau, 0, reductions))
break;
tau = newTau;
time = startTime + tau;
}
tau = newTau;
TransientSolutionInfo(0, time, tau);
time = startTime + tau;
}
}
//-------------------------------------------------------------------------
// -- static solution for tau = 0 to supply initial values
void TransientProblem:: staticInitialSolution(Real* UNorm, Real TOL)
{
int i;
Bool absPrecision = False;
Real precision;
time = startTime;
transientMode = False;
// -- one first coarse solution to get UNorm for absolute precision
staticSolution(0.05, absPrecision, time);
*UNorm = compUNorm(u);
cout << Form("\n Static Initial Solution: UNorm = %1.3g\n", *UNorm);
absPrecision = True;
precision = sqr(rhoSpace*TOL*(*UNorm));
staticSolution(precision, absPrecision, time);
TransientSolutionInfo(0, time, 0.0);
time = startTime + tau;
uPrev.resize(u.high());
FORALL(u,i) uPrev[i] = u[i]; // save solution for next step
shiftMesh();
transientMode = True;
}
//-------------------------------------------------------------------------
Real TransientProblem:: compUNorm(Vector<Num>& v, Vector<Num>* MDiag)
{
int i;
Real minU = ::real(v[1]);
FORALL(v,i) { if (::real(v[i]) < minU) minU = ::real(v[i]); }
// FORALL(v,i) { if (::real(v[1]) < minU) minU = ::real(v[1]); }
FORALL(v,i) v[i] -= minU;
Real UNorm = 0.0;
if (MDiag) FORALL(v,i) UNorm += Abs((*MDiag)[i] * v[i]*v[i]);
else UNorm = 2.0*L2MassNorm(v);
FORALL(v,i) v[i] += minU;
return sqrt(UNorm);
}
//-------------------------------------------------------------------------
Real TransientProblem:: 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->assembleP(*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;
}
//-------------------------------------------------------------------------
Bool TransientProblem:: timeConvergence(Real TOL, Real timeError, Real UNorm,
Real* newTau, int timeStep,
int reductions)
{
Bool status;
Real totalError, spaceError;
if (timeError==0) *newTau = maxTau;
else *newTau = tau * sqrt(rhoTime*TOL*UNorm/timeError);
if (*newTau > maxTau) *newTau = maxTau;
spaceError = sqrt(2.0*Error.Top());
totalError = timeError + spaceError;
if ( (totalError < TOL*UNorm) || (totalError==0) ) status = True;
else
{
status = False;
*newTau *= 0.5;
if (reductions > maxReductions)
{
cout << "\n*** Time Step reductions > maxReductions ("
<< maxReductions << ")\n"; cout.flush(); abort();
}
}
timeErrorInfo(totalError, timeError, spaceError, UNorm, TOL, tau, *newTau,
timeStep, reductions);
return status;
}
//-------------------------------------------------------------------------
void TransientProblem:: timeErrorInfo(Real /*totalError*/, Real timeError,
Real spaceError, Real UNorm, Real TOL,
Real tau0, Real newTau, int timeStep,
int reductions)
{
const char *format0 = "\n% 5s %10s %10s %7s %7s %8s %8s %3s %6s";
const char *format =
"\n %5d %8.3g %% %8.3g %% %7.3g %7.3g %8.3g %8.3g %3d %6d\n";
Real absTOL = UNorm*TOL;
cout << Form(format0, " TStep", "TimeError", "SpaceError", "TOL",
"Norm", "tau", "new/tau", "red", "nodes");
cout << Form(format, timeStep, quot(100*timeError,absTOL),
quot(100*spaceError,absTOL), TOL, UNorm, tau0,
quot(newTau,tau0), reductions, interface->Dim());
cout << ' '; for (int i=1; i<=76; ++i) cout << '~'; cout << "\n";
// Pause();
}
//-------------------------------------------------------------------------
void TransientProblem:: TransientSolutionInfo(int timeStep, Real time0,
Real tau0)
{
const int step = Error.high();
Real error = sqrt(quot(Error.Top(),Energy.Top()));
const char *format0 = "\n %8s %10s %8s %5s %5s %9s %10s %10s";
const char *format = "\n %8d %10.3g %8.3g %5d %5d %9d %10.6g %10.3g\n";
cout << Form(format0,"TimeStep", "time(sec)", " new tau",
"Steps", "level", "nodes", "ENorm", "rel.Error");
cout << Form(format, timeStep, time0, tau0,
step, interface->MaxDepth(), interface->Dim(),
sqrt(2.0*Energy.Top()), error);
cout << ' '; for (int i=1; i<=76; ++i) cout << '='; cout << "\n";
if (plotTimeStep) interface->plot (u, step, fileName, timeStep);
if (writeTimeStep) interface->write(u, step, fileName);
if (postTimeStep > 0)
if (timeStep%postTimeStep == 0)
interface->autoPost(u, step, fileName, timeStep);
while ( Pause() == picture )
interface->post (u, step, fileName);
// Pause();
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
// -- transport Solution of previous time step to the nodes of the new mesh
void TransientProblem:: solutionToNewMesh(Vector<Num>& uPrevOnNewMesh0,
const Interface* interfaceDLY) const
{
Timer timer;
int i, k, n, node, comp;
Vector<Real> x(spaceDim);
const Interface* interf = interface;
if (interfaceDLY) interf = interfaceDLY;
const int lowNode = uPrevOnNewMesh0.low();
const int highNode = uPrevOnNewMesh0.high();
NodeCoordinates nc(spaceDim, lowNode, highNode);
interf->getNodeCoordinates(nc);
if (PrevMesh() == 0) // the first time step
{
node = lowNode;
for (n=lowNode; n<=highNode; n+=nComp)
{
nc.getCoordinate(n,x);
for (comp=1; comp<=nComp; ++comp)
uPrevOnNewMesh0[node++] = dirichletBCs->initialValue(x, startTime, comp);
}
}
else
{
Vector<Real> xUnit(spaceDim);
Vector<int> prevNodes(element->NoOfNodes());
Vector<Num> uLocal(element->NoOfNodes());
Vector<int> nodes(interf->element->NoOfNodes());
const PATCH* patch, *prevPatch;
Vector<Bool> done(lowNode, highNode);
FORALL(done,i) done[i] = False;
MESH* mesh = Mesh();
mesh->resetElemIter();
while ((patch = mesh->elemIterAll()))
{
interf->getGlobalNodes(patch, nodes);
for (n=1; n<=nodes.high(); n+=nComp)
{
if ((node=nodes[n]) >= lowNode)
{
if (!done[node])
{
nc.getCoordinate(node, x);
prevPatch = PrevMesh()->findPatch(x, xUnit, patch);
interface->getGlobalNodes(prevPatch, prevNodes);
FORALL(uLocal,k) uLocal[k] = uPrev[prevNodes[k]];
if (nComp==1) uPrevOnNewMesh0[node] =
element->valueAt(xUnit, uLocal);
else element->valueAt(xUnit, uLocal, uPrevOnNewMesh0,
node, nComp);
done[node] = True;
}
}
}
}
}
if (Cmd.isTrue("TimeTransport"))
{ cout << "\n\tSolution Transport: "; timer.cpu(); }
/*
FORALL(uPrevOnNewMesh0,node)
{
nc.getCoordinate(node,x);
const PATCH* patch = PrevMesh()->findPatch(x, xUnit, 0);
interface->getGlobalNodes(patch, prevNodes);
FORALL(prevNodes,i) uLocal[i] = uPrev[prevNodes[i]];
uPrevOnNewMesh0[node] = element->valueAt(xUnit, uLocal);
}
*/
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
TransientHeatConduction:: TransientHeatConduction() { }
//-------------------------------------------------------------------------
DirichletBCs* TransientHeatConduction:: newDirichletBCs()
{
if (Cmd.isSet("DirichletBCs","ConstDirichlet"))
return new ConstDirichletBCs(fileName);
else if (Cmd.isSet("DirichletBCs","ConstTransDirichlet"))
return new ConstTransDirichlet(fileName);
else if (Cmd.isSet("DirichletBCs","StepDirichlet"))
return new StepDirichlet();
else if (Cmd.isSet("DirichletBCs","UserTransient"))
return new UserTransDirichlet();
else if (Cmd.isSet("DirichletBCs","JumpDirichlet"))
return new JumpDirichlet();
else if (Cmd.isSet("DirichletBCs","TransDirichlet"))
return new TransDirichlet();
else MissingParameter("DirichletBCs");
return 0;
}
//-------------------------------------------------------------------------
void TransientHeatConduction:: newMaterial()
{
delete material;
if (Cmd.isSet("Material","UserTransient"))
material = new UserTransMaterial(fileName, spaceDim);
else if (Cmd.isSet("Material","StepSource"))
material = new StepSource(fileName, spaceDim);
else if (Cmd.isSet("Material","DefaultMaterial"))
material = new DefaultMaterial(fileName, spaceDim);
else if (Cmd.isSet("Material","TransPeakSource"))
material = new TransPeakSource(fileName, spaceDim);
else MissingParameter("Material");
}
//-------------------------------------------------------------------------
void TransientHeatConduction:: assembleGlobal(Real time0)
{
if (transientMode) // initial static solution possible
{
uPrevOnNewMesh.resize(1,interface->Dim());
solutionToNewMesh(uPrevOnNewMesh);
}
Problem::assembleGlobal(time0);
}
//-------------------------------------------------------------------------
void TransientHeatConduction:: assemble(const Element& elem, const PATCH& t,
const Jacobian& Jac, Matrix<Num>& A,
Vector<Num>& b,
const Matrix<Bool>* pattern,
Bool errorEstimatorCall)
{
if (!transientMode) // initial static Solution
{
elem.assembleEllip (t,A,Jac,pattern);
if (lumpedMass) elem.assembleLumpedMass(t,A,Jac,pattern);
else elem.assembleMass(t,A,Jac,pattern);
elem.assembleSource(t,b,Jac,pattern,time);
if (t.onBoundary())
{
elem.assembleCauchyBCs(t,A,b,pattern,time);
elem.assembleNeumannBCs(t,b,pattern,time);
}
//if (Cmd.isTrue("innerBoundary")) elem.assembleInnerBCs(t,b,pattern,time);
if (Mesh()->innerBoundary) elem.assembleInnerBCs(t,b,pattern,time);
return;
}
int i, k;
Real factor;
const int dim = elem.NoOfNodes();
Matrix<Real> E(dim,dim); // 'elliptic' sub-matrix
Matrix<Real> P(dim,dim); // 'parabolic' sub-matrix (mass-matrix)
Vector<Real> s(dim), s0(dim);
for (i=1; i<=dim; ++i)
for (k=1; k<=i; ++k) E(i,k) = P(i,k) = 0.0;
if (!elem.assembleEllip(t,E,Jac,pattern))
missingMaterialTerm("elliptic", t.Class());
if (lumpedMass) elem.assembleLumpedMass(t,E,Jac,pattern);
else elem.assembleMass(t,E,Jac,pattern);
if (!errorEstimatorCall)
{
Bool pTerm = False;
if (lumpedMass) pTerm = elem.assembleLumpedP(t,P,Jac,pattern);
else pTerm = elem.assembleP(t,P,Jac,pattern);
if (!pTerm) missingMaterialTerm("mass (heat capacitance)",t.Class());
}
else elem.assembleP(t,P,Jac,pattern); // no lumping for error estimator!
// -- the contribution of the previous step to the right-hand side:
if (!errorEstimatorCall) // ! A is used as working space
{
Vector<int> nodes(dim);
interface->getGlobalNodes(&t,nodes);
if (implicitEuler)
{
if (lumpedMass) FORALL(b,i) b[i] += P(i,i)*uPrevOnNewMesh[nodes[i]];
else
{
for (i=1; i<=dim; ++i)
for (k=1; k< i; ++k) P(k,i) = P(i,k);
for (i=1; i<=dim; ++i)
for (k=1; k<=dim; ++k) b[i] += P(i,k)*uPrevOnNewMesh[nodes[k]];
}
}
else
{
factor = tau*(1.0-theta);
for (i=1; i<=dim; ++i)
for (k=1; k<=i; ++k)
{
A(i,k) = P(i,k) - factor*E(i,k);
A(k,i) = A(i,k);
}
for (i=1; i<=dim; ++i)
for (k=1; k<=dim; ++k) b[i] += A(i,k) * uPrevOnNewMesh[nodes[k]];
}
}
// -- the stiffness matrix:
factor = tau*theta;
for (i=1; i<=dim; ++i)
for (k=1; k<=i; ++k) A(i,k) = P(i,k) + factor*E(i,k);
// -- the source term:
if (material->SourceTerm(t.Class()))
{
FORALL(s,i) s[i] = 0.0;
if (implicitEuler) elem.assembleSource(t,s,Jac,pattern,time);
else if (midPointRule) elem.assembleSource(t,s,Jac,pattern,time-0.5*tau);
else
{
FORALL(s0,i) s0[i] = 0.0;
elem.assembleSource(t,s0,Jac,pattern,time-tau);
elem.assembleSource(t,s, Jac,pattern,time);
factor = 1.0-theta;
FORALL(s,i) s[i] = theta*s[i] + factor*s0[i];
}
FORALL(s,i) b[i] += tau*s[i];
}
// -- the boundary conditions:
if (t.onBoundary())
{
Bool cauchyBCs = False;
Bool neumBCs = False;
for (i=1; i<=dim; ++i)
{
s[i] = 0.0;
for (k=1; k<=i; ++k) P(i,k) = E(i,k) = 0.0;
}
if (implicitEuler)
{
cauchyBCs = elem.assembleCauchyBCs(t,P,s,pattern,time);
neumBCs = elem.assembleNeumannBCs(t,s,pattern, time);
}
else if (midPointRule)
{
Real midTime = time - 0.5*tau;
cauchyBCs = elem.assembleCauchyBCs(t,P,s,pattern,midTime);
neumBCs = elem.assembleNeumannBCs(t,s,pattern, midTime);
}
else
{
FORALL(s0,i) s0[i] = 0.0;
neumBCs = elem.assembleNeumannBCs(t,s0,pattern,time-tau);
neumBCs = elem.assembleNeumannBCs(t,s,pattern, time);
cauchyBCs = elem.assembleCauchyBCs(t,E,s0,pattern,time-tau);
cauchyBCs = elem.assembleCauchyBCs(t,P,s,pattern, time);
factor = 1.0-theta;
if (cauchyBCs || neumBCs) s[i] = theta*s[i] + factor*s0[i] ;
if (cauchyBCs)
{
for (i=1; i<=dim; ++i)
for (k=1; k<=i; ++k) P(i,k) = theta*P(i,k) + factor*E(i,k);
}
}
if (cauchyBCs || neumBCs) FORALL(b,i) b[i] += tau*s[i];
if (cauchyBCs)
{
for (i=1; i<=dim; ++i)
for (k=1; k<=i; ++k) A(i,k) += tau*P(i,k);
}
}
if (Mesh()->innerBoundary)
{
FORALL(s,i) s[i] = 0.0;
if (implicitEuler) elem.assembleInnerBCs(t,s,pattern,time);
else if (midPointRule) elem.assembleInnerBCs(t,s,pattern,time-0.5*tau);
else
{
FORALL(s0,i) s0[i] = 0.0;
elem.assembleInnerBCs(t,s0,pattern,time-tau);
elem.assembleInnerBCs(t,s, pattern,time);
factor = 1.0-theta;
FORALL(s,i) s[i] = theta*s[i] + factor*s0[i];
}
FORALL(s,i) b[i] += tau*s[i];
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
// -- assemble right-hand side for time error estimation and
// multiplicative defect correction
void TransientHeatConduction:: assembleGlobalDefect(MLMatrix& AD,
Vector<Num>& MD,
Vector<Num>& bD)
{
int i, node;
PATCH* patch;
const int dim = element->NoOfNodes();
Matrix<Num> AElem(dim,dim);
Vector<Num> bElem(dim), MElem(dim);
Vector<int> nodes(dim);
Jacobian Jac;
AD.reset();
FORALL(bD,i) bD[i] = MD[i] = 0.0;
Mesh()->resetElemIter();
while ((patch = Mesh()->elemIterAll()))
{
patch->compJinv(Jac);
element->initAb(AElem,bElem);
assembleDefect(*element, *patch, Jac, AElem, MElem, bElem);
interface->getGlobalNodes(patch, nodes);
AD.store(AElem, nodes);
FORALL(nodes,i)
if ((node=nodes[i]))
{
bD[node] += bElem[i];
MD[node] += MElem[i];
}
}
}
//-------------------------------------------------------------------------
void TransientHeatConduction:: assembleDefect(const Element& elem,
const PATCH& t,
const Jacobian& Jac,
Matrix<Num>& AD, Vector<Num>& MD,
Vector<Num>& bD)
{
int i, k;
const int dim = elem.NoOfNodes();
Matrix<Real> P(dim,dim); // 'parabolic' sub-matrix (mass-matrix)
Vector<Real> s0(dim);
for (i=1; i<=dim; ++i)
for (k=1; k<=i; ++k) P(i,k) = 0.0;
elem.assembleEllip(t,AD,Jac,0);
if (lumpedMass)
{
elem.assembleLumpedP(t,P,Jac,0);
FORALL(MD,i) MD[i] = P(i,i);
}
else
{
elem.assembleP(t,P,Jac,0);
for (i=1; i<=dim; ++i)
for (k=1; k< i; ++k) P(k,i) = P(i,k);
for (i=1; i<=dim; ++i) // lump Mass by row sum
{
MD[i] = 0.0;
for (k=1; k<=dim; ++k) MD[i] += P(i,k);
}
}
if (material->SourceTerm(t.Class()))
{
FORALL(s0,i) s0[i] = 0.0;
elem.assembleSource(t,bD,Jac,0,time);
elem.assembleSource(t,s0,Jac,0,time-tau);
FORALL(bD,i) bD[i] -= s0[i];
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
//#include "feplot1.h"
//#include "triang1.h"
void TransientProblem:: timeErrorByMidPointRule(Real* timeError,
Real* UNorm, Real TOL)
{
int i;
*UNorm = compUNorm(u);
Vector<Num> u2(u.high());
FORALL(u,i) u2[i] = u[i];
const Bool midPointRule0 = midPointRule;
const Real theta0 = theta;
midPointRule = True;
theta = 0.5;
Problem::assembleGlobal(time);
midPointRule = midPointRule0;
theta = theta0;
Real precision = sqr(rhoSpace*TOL);
precision *= quot(sqr(*UNorm),Energy.Top());
int step = Energy.high();
Ab->solve(u2, 0.0, 0.0, precision, step);
FORALL(u,i) u2[i] -= u[i];
*timeError = compUNorm(u2);
}
//-------------------------------------------------------------------------
// -- this routine uses a multiplicative defect correction scheme
// for time error estimation
void TransientProblem:: compTimeError(Real* timeError, Real* UNorm, Real TOL)
{
if (Cmd.isSet("TError", "midPointRule")) // for test purposes
{
timeErrorByMidPointRule(timeError, UNorm, TOL);
return;
}
static Real accTime = 0.0;
Timer timer, accTimer;
int i;
const int dim = u.high();
MLMatrix* AD;
if (nComp == 1) AD = new MLSparseMatrix(sym, SpaceDim(), dim);
else AD = new MLBlockMatrix (sym, SpaceDim(), dim, nComp);
Vector<Num> u2(dim), MD(dim), bD(dim), eta(dim);
assembleGlobalDefect(*AD,MD,bD);
*UNorm = compUNorm(u, &MD);
FORALL(bD,i) u2[i] = u[i] - uPrevOnNewMesh[i];
AD->Mult(eta,u2);
FORALL(bD,i) bD[i] = 0.5*tau * (eta[i] - bD[i]); // new right hand side
FORALL(bD,i) if (dirichletBCs->isSet(i)) bD[i] = 0.0;
dirichletBCs->setValuesToZero();
Ab->setRhs(bD);
Real precision = sqr(rhoSpace*TOL);
precision *= quot(sqr(*UNorm) , Energy.Top());
int step = Energy.high();
FORALL(eta,i) eta[i] = 0.0;
Ab->solve(eta, 0.0, 0.0, precision, step);
*timeError = compUNorm(eta, &MD);
if (Cmd.isTrue("timeTimeError"))
{ cout << "\n\tTimeError Estimation: "; timer.cpu(); }
if (Cmd.isTrue("accTimeTimeError"))
{
accTime += accTimer.cpu(False);
cout << Form("\tAccumulated time: %1.2f sec.\n", accTime);
}
delete AD;
//------------------------- test output ------------------------
/*
MESH1 *mesh = (MESH1*) Mesh();
FEPlotMESH1 Plot(mesh, X11, "DefectCorrection", 0.4);
Plot.plotElements();
Plot.plotBoundary();
FORALL(bD,i) u2[i] = u[i] + eta[i];
Plot.plotSolution(eta);
Pause();
Real eNorm = sqrt(2.0*EnergyNorm(eta));
cout << "\n--- Time Error: eNorm=" << eNorm << " L2Mass= "
<< *timeError << " E/M = " << quot(eNorm,(*timeError)) << "\n";
AD.Mult(bD,eta);
eNorm = sqrt(cdot(bD,eta));
Real tError = sqrt(2.0*L2MassNorm(eta));
cout << "--- Time Error: ellipNorm=" << eNorm << " L2Mass= " << tError
<< " Ell/M = " << quot(eNorm,tError) << "\n";
AD.Mult(bD,u);
eNorm = sqrt(2.0*EllipticNorm(u));
cout << "--- Time Sol: eNorm=" << eNorm << " UNorm= " << *UNorm
<< " E/M = " << quot(eNorm,(*UNorm)) << "\n";
*/
}
syntax highlighted by Code2HTML, v. 0.9.1