/*
$Id: nonlin.cc,v 1.2 1996/10/04 15:07:10 roitzsch Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "nonlin.h"
#include "numerics.h"
#include "int.h"
#include "nodeco.h"
#include "cmdpars.h"
extern CmdPars Cmd;
//-------------------------------------------------------------------------
static const Real MachMax = machMax(Real(0.0));
//-------------------------------------------------------------------------
NonLinearity:: NonLinearity() : u(1), H(1), areaWeight(1),
uppObstacle(1), lowObstacle(1)
{ }
NonLinearity:: ~NonLinearity() { }
//-------------------------------------------------------------------------
Real NonLinearity:: lowerObstacle(Vector<Real>& /*x*/) const { return -MachMax; }
Real NonLinearity:: upperObstacle(Vector<Real>& /*x*/) const { return MachMax; }
//-------------------------------------------------------------------------
void NonLinearity:: addAreaWeight(const Vector<Real>& LNorm,
const Vector<int>& nodes)
{
int i;
FORALL(nodes,i) { if (nodes[i]) areaWeight[nodes[i]] += LNorm[i]; }
}
//-------------------------------------------------------------------------
void NonLinearity:: update(const Interface& interface)
{
int i;
areaWeight.resize(interface.Dim());
FORALL(areaWeight,i) areaWeight[i] = 0.0;
updateObstacles(uppObstacle, lowObstacle, interface);
}
//-------------------------------------------------------------------------
void NonLinearity:: defaultObstacleUpdate(Vector<Real>& uppO,Vector<Real>& lowO,
const Interface& interface) const
{
int i;
uppO.resize(interface.Dim());
lowO.resize(interface.Dim());
FORALL(uppO,i)
{
lowO[i] = -MachMax;
uppO[i] = MachMax;
}
}
//-------------------------------------------------------------------------
void NonLinearity:: fctObstacleUpdate(Vector<Real>& uppO, Vector<Real>& lowO,
const Interface& interface) const
{
int i;
const int dim = interface.Dim();
const int spaceDim = interface.SpaceDim();
uppO.resize(dim);
lowO.resize(dim);
Vector<Real> x(spaceDim);
NodeCoordinates nc(spaceDim, 1, dim);
interface.getNodeCoordinates(nc);
FORALL(uppO,i)
{
nc.getCoordinate(i,x);
lowO[i] = lowerObstacle(x);
uppO[i] = upperObstacle(x);
}
}
//-------------------------------------------------------------------------
Real NonLinearity:: HValue(int node, Real uVal) const
{
if (node <= lowObstacle.high())
{
if (uVal < lowObstacle[node])
{
cout << "\n** NonLinearity::HValue(...): value " << uVal
<< " of node " << node << " below lower obstacle: "
<< lowObstacle[node] << "\n";
//cout.flush(); abort();
}
if (uVal > uppObstacle[node])
{
cout << "\n** NonLinearity::HValue(...): value " << uVal
<< " of node " << node << " above upper obstacle: "
<< uppObstacle[node] << "\n";
//cout.flush(); abort();
}
}
return HValue(uVal);
}
//-------------------------------------------------------------------------
Real NonLinearity:: HValue(Real uVal) const
{
int i;
if (uVal <= u[1]) return linearInterp(uVal, u[1], u[2], H[1], H[2]);
else if (uVal >= u[M]) return linearInterp(uVal, u[M-1], u[M], H[M-1], H[M]);
else
{
for (i=2; i<=M; ++i)
{
if (uVal <= u[i])
{
if (equal(uVal,u[i]) && equal(u[i-1],u[i]))
return 0.5 * (H[i-1] + H[i]);
else return linearInterp(uVal, u[i-1], u[i], H[i-1], H[i]);
}
}
}
cout << "\n*** NonLinearity::HValue(...): uVal " << uVal << " not located\n";
cout.flush(); abort();
return 0.0;
}
//-------------------------------------------------------------------------
// -- interpolate H = H(u), H1 = H(u1), H2 = H(u2)
Real NonLinearity:: linearInterp(Real ui, Real u1, Real u2, Real H1, Real H2)
const
{
Real x = (ui-u1)/(u2-u1);
return H1 + x*(H2-H1);
}
//-------------------------------------------------------------------------
Real PorousMedia:: HValue(Real uVal) const
{
int i;
if (uVal <= u[1])
{
if (uVal > 0.0) return pow(uVal, 1.0/mm);
else return 0.0;
}
else if (uVal >= u[M]) return linearInterp(uVal, u[M-1], u[M], H[M-1], H[M]);
else
{
for (i=2; i<=M; ++i)
{
if (uVal <= u[i])
{
if (equal(uVal,u[i]) && equal(u[i-1],u[i]))
return 0.5 * (H[i-1] + H[i]);
else return linearInterp(uVal, u[i-1], u[i], H[i-1], H[i]);
}
}
}
cout << "\n*** NonLinearity::HValue(...): uVal " << uVal << " not located\n";
cout.flush(); abort();
return 0.0;
}
//-------------------------------------------------------------------------
Casting:: Casting() : NonLinearity()
{
kappa0 = 7.0/20000.0; Cmd.get("kappa0", &kappa0);
kappa1 = 7.0/4000.0; Cmd.get("kappa1", &kappa1);
M=8;
u.resize(M);
H.resize(M);
/*
u[1]= -0.315; H[1]= 3.27;
u[2]= -0.0105; H[2]= 6.94;
u[3]= -0.0105; H[3]= 7.15;
u[4]= -0.007; H[4]= 7.30;
u[5]= -0.007; H[5]= 7.64;
u[6]= 0.0; H[6]= 7.92;
u[7]= 0.0; H[7]= 9.44;
u[8]= 0.175; H[8]= 9.84;
*/
u[1]= -0.45; H[1]= 3.27;
u[2]= -0.015; H[2]= 6.94;
u[3]= -0.015; H[3]= 7.15;
u[4]= -0.01; H[4]= 7.30;
u[5]= -0.01; H[5]= 7.64;
u[6]= 0.0; H[6]= 7.92;
u[7]= 0.0; H[7]= 9.44;
u[8]= 0.175; H[8]= 9.84;
}
//-------------------------------------------------------------------------
Stefan:: Stefan() : NonLinearity()
{
kappa0 = 1.0; Cmd.get("kappa0", &kappa0);
kappa1 = 1.0; Cmd.get("kappa1", &kappa1);
c0 = 1.0; Cmd.get("c0", &c0);
c1 = 1.0; Cmd.get("c1", &c1);
s0 = 0.0; Cmd.get("s0", &s0);
s1 = 0.0; Cmd.get("s1", &s1);
theta1 = 0.0; Cmd.get("theta1", &theta1);
Real a0 = c0 / kappa0;
Real a1 = c1 / kappa1;
M=4;
u.resize(M);
H.resize(M);
u[1]= theta1 - 1.0; H[1]= -(a0+s0);
u[2]= theta1; H[2]= -s0;
u[3]= theta1; H[3]= s1;
u[4]= theta1 + 1.0; H[4]= a1+s1;
cout << "\nStefan constructed\n";
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
Obstacle:: Obstacle() : NonLinearity()
{
M=2;
u.resize(M);
H.resize(M);
u[1]=-1.0; H[1]= 0.0;
u[2]= 0.0; H[2]= 0.0;
cout << "\nObstacle constructed\n";
}
Real Obstacle:: lowerObstacle(Vector<Real>& /*xvec*/) const { return -MachMax; }
Real Obstacle:: upperObstacle(Vector<Real>& xvec) const
{
Real x = xvec[1];
Real y = xvec[2];
if (y < 0.5)
{
if (x < 0.5) return (x<y)? x:y;
else return ((1-x)<y)? (1-x):y;
}
else
{
if (x < 0.5) return (x<1.-y)? x:(1.-y);
else return ((1-x)< (1.-y) )? (1.-x):(1.-y);
}
//return 1000.;
}
//-------------------------------------------------------------------------
Real Obstacle:: HValue(int /*node*/, Real /*uVal*/) const
{
cout << "\n*** Obstacle::HValue(...): H makes no sense\n";
cout.flush(); abort();
return 0.0;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
PorousMedia:: PorousMedia() : NonLinearity(), critic(1)
{
int i;
Real hM, rightEnd=1.;
mm = 2.;
M=100;
u.resize(M);
H.resize(M);
hM = pow(rightEnd,1./mm)/M;
for (i=1; i<=M; i++)
{
u[i] = pow(i*hM,mm);
H[i] = pow(u[i],1./mm);
}
u[1] = u[2]/2.;
H[1] = pow(u[1],1./mm);
cout << "\nPorousMedia constructed\n";
}
//-------------------------------------------------------------------------
Real PorousMedia:: lowerObstacle(Vector<Real>& /*xvec*/) const { return 0.0; }
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void PorousMedia:: update(const Interface& interface)
{
int i;
NonLinearity::update(interface);
critic.resize(interface.Dim());
FORALL(critic,i) critic[i] = True; // to include dirichlet nodes
}
//-------------------------------------------------------------------------
void Obstacle:: defectCorrection(Real& e, Real& r, Real oldSol,
Real diag, Real& a, Real& uppDefO,
Real& lowDefO, Bool& critical,
Real uppO, Real lowO, Real /*arWeight*/) const
{
Real newSol = oldSol + r/diag;
critical = False;
if (newSol <= lowO)
{
newSol = lowO;
lowDefO = 0.0;
uppDefO = 0.0;
critical = True;
}
else if (newSol >= uppO)
{
newSol = uppO;
uppDefO = 0.0;
lowDefO = 0.0;
critical = True;
}
else
{
lowDefO = lowO-newSol;
uppDefO = uppO-newSol;
}
a = 0.0;
r = 0.0;
e = newSol - oldSol;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void NonLinearity:: defectCorrection(int node, Real& e, Real& r, Real oldSol,
Real diag)
{
Real a, uppDefO, lowDefO; // dummies
Bool critical = True;
defectCorrection(e, r, oldSol, diag, a, uppDefO, lowDefO, critical,
uppObstacle[node], lowObstacle[node], areaWeight[node]);
if (critical) setCritical(node);
else setUnCritical(node);
}
//-------------------------------------------------------------------------
void NonLinearity:: defectCorrection(int node, Real& e, Real& r, Real oldSol,
Real diag, Real& a, Real& uppDefO,
Real& lowDefO, Bool& critical)
{
defectCorrection(e, r, oldSol, diag, a, uppDefO, lowDefO, critical,
uppObstacle[node], lowObstacle[node], areaWeight[node]);
if (critical) setCritical(node);
else setUnCritical(node);
}
//-------------------------------------------------------------------------
// -- apply nonlinear contributions:
// -> defect correction e, residual correction r;
// for multi-grid: correction a on diagonal, obstacles, and critical flag
//
void NonLinearity:: defectCorrection(Real& e, Real& r, Real oldSol,
Real diag, Real& a, Real& uppDefO,
Real& lowDefO, Bool& critical,
Real uppO, Real lowO, Real arWeight) const
{
int k;
Real newSol, b, u0, u1, H0, H1, hSquare, s=0.0;
critical = False;
a = 0.0;
b = r + diag*oldSol;
hSquare = arWeight;
H0 = H[1]; u0 = u[1];
H1 = H[2]; u1 = u[2];
if (b < H1*hSquare+diag*u1) // the first interval
{
a = (H1-H0)/(u1-u0);
if (lowO > u0) cout << "Help!!!\n";
u0 = lowO;
H0 = a*(u0-u1) + H1;
if (b < H0*hSquare + diag*u0)
{
a=0.;
s=0.;
newSol = u0;
lowDefO = 0.0;
uppDefO = 0.0;
critical = True;
}
else
{
a = a*hSquare; s = H1*hSquare-a*u1;
newSol = (b-s)/(diag+a);
lowDefO = u0-newSol;
uppDefO = u1-newSol;
}
e = newSol - oldSol;
}
else
{
for (k=2; k < M-1; k++)
{
H0 = H[k]; u0 = u[k];
H1 = H[k+1]; u1 = u[k+1];
if ((H0*hSquare + diag*u0 <= b)&&(b < H1*hSquare+diag*u1) )
{
if (equal(u0,u1))
{
a=0.; s=0.;
newSol = u0;
lowDefO = 0.0;
uppDefO = 0.0;
critical = True;
}
else
{
a = (H0-H1)*hSquare/(u0-u1); s = H0*hSquare-a*u0;
newSol = (b-s)/(diag+a);
lowDefO = u0 - newSol;
uppDefO = u1 - newSol;
}
e = newSol - oldSol;
//if (a < theta)
//{
//lowDefO = 0.0;
//uppDefO = 0.0;
//critical = True;
//}
break;
}
}
}
H0 = H[M-1]; u0 = u[M-1];
H1 = H[M]; u1 = u[M];
if (H0*hSquare+diag*u0 <= b) // the last interval
{
a = (H1-H0)/(u1-u0);
if (uppO < u1) cout << "MayDay\n";
u1 = uppO;
H1 = a*(u1 - u0) + H0;
if (H1*hSquare+diag*u1 <= b)
{
a=0.; s=0.;
newSol = u1;
lowDefO = 0.0;
uppDefO = 0.0;
critical = True;
}
else
{
a = a*hSquare; s = H0*hSquare-a*u0;
newSol = (b-s)/(diag+a);
lowDefO = u0-newSol;
uppDefO = u1-newSol;
}
e = newSol - oldSol;
}
r = s + a*oldSol;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
NLTest:: NLTest() : NonLinearity()
{
M=2;
u.resize(M);
H.resize(M);
u[1]=-1.0; H[1]= 0.0;
u[2]= 0.0; H[2]= 0.0;
m = 1000.0;
cout << "\nNLTest constructed\n";
Cmd.get("MValue", &m);
}
//-------------------------------------------------------------------------
void NLTest:: updateObstacles(Vector<Real>& uppO, Vector<Real>& lowO,
const Interface& interface) const
{
int i;
uppO.resize(interface.Dim());
lowO.resize(interface.Dim());
FORALL(uppO,i) { lowO[i] = -10000.0; uppO[i] = 10000.0; }
}
//-------------------------------------------------------------------------
Real NLTest:: HValue(int /*node*/, Real uVal) const { return HValue(uVal); }
Real NLTest:: HValue(Real uVal) const
{
return m * uVal * uVal;
}
//-------------------------------------------------------------------------
void NLTest:: defectCorrection(Real& e, Real& r, Real oldSol,
Real diag, Real& a, Real& uppDefO,
Real& lowDefO, Bool& critical,
Real uppO, Real lowO, Real arWeight) const
{
Real newSol;
Real dH = m*arWeight * oldSol;
Real H0 = dH * oldSol;
dH *= 2.0;
newSol = oldSol + (r - H0)/(diag + dH);
critical = False;
if (newSol <= lowO)
{
newSol = lowO;
lowDefO = 0.0;
uppDefO = 0.0;
critical = True;
}
else if (newSol >= uppO)
{
newSol = uppO;
uppDefO = 0.0;
lowDefO = 0.0;
critical = True;
}
else
{
lowDefO = lowO-newSol;
uppDefO = uppO-newSol;
}
e = newSol - oldSol;
a = dH;
r = H0;
}
//-------------------------------------------------------------------------
syntax highlighted by Code2HTML, v. 0.9.1