Scalar code example 1.
Jean Marshall
jcm at lanl.gov
Sun Apr 22 00:10:05 UTC 2001
An HTML attachment was scrubbed...
URL: <http://sourcerytools.com/pipermail/pooma-dev/attachments/20010421/a7c1eb42/attachment.html>
-------------- next part --------------
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// CCSource File: QRelations
// Author: jcm
// Date: Sat - Nov 18, 2000
// Namespace: conejo
// Framework: Tecolote
// Copyright: Los Alamos National Laboratory
// Full Copyright=$(TECOLOTE_ROOT)/Doc/Copyright
// RCS_VERSION_ID: $Id:
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#if !defined (__MWERKS__)
#include "Demo_src/include/DemoPCH.hh"
#pragma hdrstop
#endif // !__MWERKS__
#include "Demo_src/Model/CompatibleHydro/QRelations.t.hh"
#include "TecFramework_src/MetaTypes/MetaTypes.hh"
#include "TecFramework_src/Foundation/LoadObject.hh"
namespace conejo
{
using namespace TecFramework;
using namespace poomalote;
using namespace Hydrodynamics;
using namespace PhysicsBaseClasses;
using namespace std;
static MetaClass<ITecoloteTraits<QRelations<ThreeDF<DefaultTraits> >, RelationPkg> >
QRelations3DMeta("QRelations3D", QRelations<ThreeDF<DefaultTraits> >::MakePersistents());
LoadObjectGroup QRelations3DBase_cc = { &QRelations3DMeta };
} // end namespace conejo
namespace TecFramework
{
using namespace conejo;
LoadObjectGroup QRelations3D_cc = { &QRelations3DBase_cc };
} // end namespace TecFramework
-------------- next part --------------
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// Source File: QRelations
// Author: jcm
// Date: Sat - Nov 18, 2000
// Namespace: conejo
// Framework: Tecolote
// Copyright: Los Alamos National Laboratory
// Full Copyright=$(TECOLOTE_ROOT)/Doc/Copyright
// RCS_VERSION_ID: $Id:
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#ifndef __conejo_QRelations_t_hh
#define __conejo_QRelations_t_hh
#include "PhysicsBaseClasses_src/HelperClasses/MMField.hh"
#include "Demo_src/Model/CompatibleHydro/QRelations.hh"
#include "Evaluator/ScalarCode.h"
//#define ENTER(a)
#define ENTER(a) tecout << "Entering " << a << endl;
namespace conejo
{
using namespace TecFramework;
using namespace poomalote;
using namespace PhysicsBaseClasses;
using namespace Hydrodynamics;
using namespace std;
const Real spokeCutoff = 1.0e-12; // We get this value from Ed's code
// $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
template <class Traits>
BEGIN_PERSISTENT(QRelations<Traits>)
PERSISTENT( linearQ, "LinearQ" )
PERSISTENT( quadQ, "QuadQ" )
END_PERSISTENT
//$ linearQ : Real - linear coefficient for Q
//
//$ quadQ : Real - quadratic coefficient for Q
//
// $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
//-----------------------------------------------------------------------------
// Test function
//-----------------------------------------------------------------------------
template<class Traits>
struct EdgeQPressureInfo
{
void scalarCodeInfo(ScalarCodeInfo &info) const
{
info.arguments(5);
info.write(0, true );
info.write(1, false);
info.write(2, false);
info.write(3, false);
info.write(4, false);
// info.write(5, false);
// info.write(6, false);
// info.write(7, false);
info.useGuards(0, false);
info.useGuards(1, true);
info.useGuards(2, true);
info.useGuards(3, true);
info.useGuards(4, true);
// info.useGuards(5, true);
// info.useGuards(6, true);
// info.useGuards(7, true);
info.dimensions(Traits::Dim);
for (int i = 0; i < Traits::Dim; ++i)
{
info.lowerExtent(i) = 2;
info.upperExtent(i) = 2;
}
}
};
template<class Traits>
struct ScalarEdgeQPressure
: public EdgeQPressureInfo<Traits>
{
// Typedefs
FIELD_TYPEDEFS(Traits)
ScalarEdgeQPressure(const Real& inLinearQ )
: EdgeQPressureInfo<Traits>(),
linearQ(inLinearQ)
{
}
template<class F1, class F2, class F3, class F4, class F5>
void operator()(const F1& EdgeQPressure, const F2& EdgeGammaConstant, const F3& EdgeSoundSpeed,
const F4& EdgeVelocity, const F5& EdgePsiLimiter, const Loc<Traits::Dim> &loc)
{
if( EdgePsiLimiter(loc) < eps ) {
EdgeQPressure(loc) = 0.0;
return;
}
Real edgeVelocityMagnitude = sqrt(dot(EdgeVelocity(loc),EdgeVelocity(loc)));
EdgeQPressure(loc) = edgeVelocityMagnitude * EdgePsiLimiter(loc) *
(EdgeGammaConstant(loc) * edgeVelocityMagnitude +
sqrt( linearQ * linearQ * EdgeSoundSpeed(loc) * EdgeSoundSpeed(loc) +
EdgeGammaConstant(loc) * EdgeGammaConstant(loc) * edgeVelocityMagnitude * edgeVelocityMagnitude));
}
private:
Real linearQ;
};
//======================================================================
// Constructor -- QRelations<Traits>::QRelations
//======================================================================
template <class Traits>
QRelations<Traits>::QRelations( DataDirectory* pDataDir, const string& inName )
: CompatibleRelations(pDataDir,inName),
Old(pDataDir->strictGet<DataDirectory>("CompatibleHydroOld")),
linearQ(1.0),
quadQ(1.0)
{
VectorField& EdgeLength = DataDir.get<VectorField>( "EdgeLength", Mesh.getField<VectorField> ( AllEdge() ) );
for(int d=0;d<Dim;++d) {
Interval<1> CVert(CoarseVert[d]);
RDomain.push_back(CoarseVert);
RDomain[d][d] = Interval<1>(CVert.first(), CVert.last() - 2);
LDomain.push_back(CoarseVert);
LDomain[d][d] = Interval<1>(CVert.first() + 1,CVert.last() - 1);
Loc<Dim> offset(0); offset[d] = 1;
RightEdgeNgbr.push_back(RDomain[d]);
RightEdgeNgbr[d] += offset;
LeftEdgeNgbr.push_back(LDomain[d]);
LeftEdgeNgbr[d] -= offset;
CoarseEdges.push_back(Range<Dim>(EdgeLength[d].domain()));
LowerVert.push_back(CoarseEdges[d]);
UpperVert.push_back(CoarseEdges[d] + offset);
EdgeNgbr.push_back( Range<Dim>(EdgeLength[d].domain() + offset) );
}
} // end constructor
//======================================================================
// Destructor -- QRelations<Traits>::~QRelations
//======================================================================
template <class Traits>
QRelations<Traits>::~QRelations( ) { }
//======================================================================
// Function -- QRelations<Traits>::createRelations
//======================================================================
template <class Traits>
void QRelations<Traits>::createRelations( ) {
ENTER("QRelations::createRelations");
gammaConst = 2.6666666666666667 * 0.25 * quadQ;
// Dt is an Independent Constant Field used in Relationships
ScalarField& Dt = DataDir.get<ScalarField>( "Dt", Mesh.getField<ScalarField> ( Vert () ) );
// -------------------------------------------------------------------------------------------------------------------------
// Input Fields to QRelations
// -------------------------------------------------------------------------------------------------------------------------
ScalarField& OldDt = Old.get<ScalarField> ( "Dt", Mesh.getField<ScalarField> ( Vert () ) );
VectorField& OldVelocity = Old.get<VectorField> ( "Velocity" , Mesh.getField<VectorField> ( Vert () ) );
VectorField& OldEdgeLength = Old.get<VectorField> ( "EdgeLength" , FineMesh.getField<VectorField>( AllEdge() ) );
VectorField& OldEdgeVelocity = Old.get<VectorField> ( "EdgeVelocity", Mesh.getField<VectorField> ( AllEdge() ) );
ScalarField& Volume = DataDir.get<ScalarField>( "Volume", Mesh.getField<ScalarField> ( Cell () ) );
ScalarField& SubVolume = DataDir.get<ScalarField>( "SubVolume", FineMesh.getField<ScalarField>( Cell () ) );
VectorField& SubFaceAreas = DataDir.get<VectorField>( "SubFaceAreas" , FineMesh.getField<VectorField>( AllFace() ) );
VectorField& EdgeLength = DataDir.get<VectorField>( "EdgeLength" , FineMesh.getField<VectorField>( AllEdge() ) );
VectorField& Velocity = DataDir.get<VectorField>( "Velocity", Mesh.getField<VectorField> ( Vert () ) );
ScalarField& CellDensity = DataDir.get<ScalarField>( "CellDensity", Mesh.getField<ScalarField> ( Cell () ) );
ScalarField& CellSoundSpeedSq = DataDir.get<ScalarField>( "CellSoundSpeedSq", Mesh.getField<ScalarField> ( Cell () ) );
ScalarField& CellGammaConstant = DataDir.get<ScalarField>( "CellGammaConstant", Mesh.getField<ScalarField> ( Cell () ) );
CellGammaConstant = gammaConst; // Should go away when we use more sophisticated EOSs.
// -------------------------------------------------------------------------------------------------------------------------
// Store the sqrt because it is used to calculate the VertSoundSpeed
// -------------------------------------------------------------------------------------------------------------------------
ScalarField& CellSoundSpeed = DataDir.get<ScalarField>( "CellSoundSpeed", Mesh.getField<ScalarField> ( Cell () ) );
// -------------------------------------------------------------------------------------------------------------------------
// Vertex Fields -- (weighted) sums of SubCell Fields
// -------------------------------------------------------------------------------------------------------------------------
ScalarField& VertVolume = DataDir.get<ScalarField>( "VertVolume", Mesh.getField<ScalarField> ( Vert () ) );
// ScalarField& WeightedVertVol = DataDir.get<ScalarField>( "WeightedVertVol", Mesh.getField<ScalarField> ( Vert () ) );
ScalarField& VertDensity = DataDir.get<ScalarField>( "VertDensity", Mesh.getField<ScalarField> ( Vert () ) );
ScalarField& VertSoundSpeed = DataDir.get<ScalarField>( "VertSoundSpeed", Mesh.getField<ScalarField> ( Vert () ) );
ScalarField& VertGammaConstant = DataDir.get<ScalarField>( "VertGammaConstant", Mesh.getField<ScalarField> ( Vert () ) );
// -------------------------------------------------------------------------------------------------------------------------
// Edge Fields that hold Q-related quantities
// -------------------------------------------------------------------------------------------------------------------------
VectorField& EdgeVelocity = DataDir.get<VectorField>( "EdgeVelocity", Mesh.getField<VectorField> ( AllEdge() ) );
ScalarField& EdgeDensity = DataDir.get<ScalarField>( "EdgeDensity", Mesh.getField<ScalarField> ( AllEdge() ) );
ScalarField& EdgeQPressure = DataDir.get<ScalarField>( "EdgeQPressure", Mesh.getField<ScalarField> ( AllEdge() ) );
ScalarField& EdgeSoundSpeed = DataDir.get<ScalarField>( "EdgeSoundSpeed", Mesh.getField<ScalarField> ( AllEdge() ) );
ScalarField& EdgeGammaConstant = DataDir.get<ScalarField>( "EdgeGammaConstant", Mesh.getField<ScalarField> ( AllEdge() ) );
ScalarField& EdgeQTmpMax = DataDir.get<ScalarField>( "EdgeQTmpMax", Mesh.getField<ScalarField> ( AllEdge() ) );
ScalarField& SpokeDVolDt = DataDir.get<ScalarField>( "SpokeDVolDt", FineMesh.getField<ScalarField>( AllFace() ) );
ScalarField& SpokeQSwitch = DataDir.get<ScalarField>( "SpokeQSwitch", FineMesh.getField<ScalarField>( AllFace() ) );
ScalarField& RightLimiterRatio = DataDir.get<ScalarField>( "RightLimiterRatio", Mesh.getField<ScalarField> ( AllEdge() ) );
ScalarField& LeftLimiterRatio = DataDir.get<ScalarField>( "LeftLimiterRatio", Mesh.getField<ScalarField> ( AllEdge() ) );
ScalarField& EdgePsiLimiter = DataDir.get<ScalarField>( "EdgePsiLimiter", Mesh.getField<ScalarField> ( AllEdge() ) );
// -------------------------------------------------------------------------------------------------------------------------
// Output Fields from QRelations
// -------------------------------------------------------------------------------------------------------------------------
ScalarField& CellQ = DataDir.get<ScalarField>( "CellQ", Mesh.getField<ScalarField> ( Cell () ) );
VectorField& SubForceQ = DataDir.get<VectorField>( "SubForceQ", FineMesh.getField<VectorField>( Cell () ) );
ScalarField& SubPressureQMod = DataDir.get<ScalarField>( "SubPressureQMod", FineMesh.getField<ScalarField>( Cell () ) );
Let::NewRelation( *(pParent)this, &Parent::calcSqrt, CellSoundSpeed, CellSoundSpeedSq );
Let::NewRelation( *(pParent)this, &Parent::calcWeightedVertAvg, VertSoundSpeed, CellSoundSpeed, SubVolume, VertVolume );
Let::NewRelation( *(pParent)this, &Parent::calcWeightedVertAvg, VertGammaConstant, CellGammaConstant, SubVolume, VertVolume );
Let::NewRelation( *this, &QRelations<Traits>::calcEdgeVelocity, EdgeVelocity, OldVelocity );
Let::NewRelation( *this, &QRelations<Traits>::calcEdgeDensity, EdgeDensity, VertDensity );
Let::NewRelation( *this, &QRelations<Traits>::calcEdgeGammaConstant, EdgeGammaConstant, VertGammaConstant );
Let::NewRelation( *this, &QRelations<Traits>::calcEdgeQPressure, EdgeQPressure, EdgeGammaConstant, EdgeSoundSpeed, EdgeVelocity, EdgePsiLimiter );
Let::NewRelation( *this, &QRelations<Traits>::calcEdgeSoundSpeed, EdgeSoundSpeed, VertSoundSpeed );
Let::NewRelation( *this, &QRelations<Traits>::calcEdgeQTmpMax, EdgeQTmpMax, EdgeDensity, EdgeQPressure );
Let::NewRelation( *this, &QRelations<Traits>::calcSpokeDVolDt, SpokeDVolDt, EdgeVelocity, SubFaceAreas );
Let::NewRelation( *this, &QRelations<Traits>::calcSpokeQSwitch, SpokeQSwitch, SpokeDVolDt, Volume, OldDt );
Let::NewRelation( *this, &QRelations<Traits>::calcRightLimiterRatio, RightLimiterRatio, OldEdgeLength, OldEdgeVelocity );
Let::NewRelation( *this, &QRelations<Traits>::calcLeftLimiterRatio, LeftLimiterRatio, OldEdgeLength, OldEdgeVelocity );
Let::NewRelation( *this, &QRelations<Traits>::calcEdgePsiLimiter, EdgePsiLimiter, RightLimiterRatio, LeftLimiterRatio, OldEdgeLength, OldEdgeVelocity, OldDt );
Let::NewRelation( *this, &QRelations<Traits>::calcCellQ, CellQ, EdgeQTmpMax, SpokeQSwitch );
Let::NewRelation( *this, &QRelations<Traits>::calcSubForceQ, SubForceQ, EdgeQTmpMax, EdgeVelocity, SpokeDVolDt, SpokeQSwitch );
Let::NewRelation( *this, &QRelations<Traits>::calcSubPressureQMod, SubPressureQMod, EdgeQPressure, SpokeQSwitch );
}
//======================================================================
// Function -- QRelations<Traits>::calcEdgeVelocity
//======================================================================
template <class Traits>
void QRelations<Traits>::calcEdgeVelocity ( const VectorField& EdgeVelocity,
const VectorField& OldVelocity)
{
ENTER("QRelations::calcEdgeVelocity");
int d = getEdgeDirection(EdgeVelocity);
EdgeVelocity(CoarseEdges[d]) = OldVelocity(EdgeNgbr[d]) - OldVelocity(CoarseEdges[d]);
} // end function calcEdgeVelocity
//======================================================================
// Function -- QRelations<Traits>::calcEdgeDensity
//======================================================================
template <class Traits>
void QRelations<Traits>::calcEdgeDensity ( const ScalarField& EdgeDensity,
const ScalarField& VertDensity)
{
ENTER("QRelations::calcEdgeDensity");
int d = getEdgeDirection(EdgeDensity);
EdgeDensity(CoarseEdges[d]) = 2.0 * ( VertDensity(CoarseEdges[d]) * VertDensity(EdgeNgbr[d]) ) /
( VertDensity(CoarseEdges[d]) + VertDensity(EdgeNgbr[d]) );
ScalarField& VertVolume = DataDir.get<ScalarField>( "VertVolume", Mesh.getField<ScalarField> ( Vert () ) );
ScalarField& CellDensity = DataDir.strictGet<ScalarField>( "CellDensity" );
} // end function calcEdgeDensity
//======================================================================
// Function -- QRelations<Traits>::calcEdgeGammaConstant
//======================================================================
template <class Traits>
void QRelations<Traits>::calcEdgeGammaConstant ( const ScalarField& EdgeGammaConstant,
const ScalarField& VertGammaConstant)
{
ENTER("QRelations::calcEdgeGammaConstant");
int d = getEdgeDirection(EdgeGammaConstant);
EdgeGammaConstant(CoarseEdges[d]) = 0.5 * ( VertGammaConstant(CoarseEdges[d]) + VertGammaConstant(EdgeNgbr[d]) );
} // end function calcEdgeGammaConstant
//======================================================================
// Function -- QRelations<Traits>::calcEdgeQPressure
//======================================================================
template <class Traits>
void QRelations<Traits>::calcEdgeQPressure ( const ScalarField& EdgeQPressure,
const ScalarField& EdgeGammaConstant,
const ScalarField& EdgeSoundSpeed,
const VectorField& EdgeVelocity,
const ScalarField& EdgePsiLimiter )
{
ENTER("QRelations::calcEdgeQPressure");
int d = getEdgeDirection(EdgeQPressure);
EdgeGammaConstant[d].update();
EdgeSoundSpeed[d].update();
EdgeVelocity[d].update();
EdgePsiLimiter[d].update();
ScalarEdgeQPressure<Traits> sEdgeQPressure(linearQ);
ScalarCode<ScalarEdgeQPressure<Traits> > scEdgeQPressure(sEdgeQPressure);
scEdgeQPressure(EdgeQPressure,EdgeGammaConstant[d], EdgeSoundSpeed[d],EdgeVelocity[d],EdgePsiLimiter[d]);
} // end function EdgeQPressure
//======================================================================
// Function -- QRelations<Traits>::calcEdgeSoundSpeed
//======================================================================
template <class Traits>
void QRelations<Traits>::calcEdgeSoundSpeed ( const ScalarField& EdgeSoundSpeed,
const ScalarField& VertSoundSpeed )
{
ENTER("QRelations::calcEdgeSoundSpeed");
int d = getEdgeDirection(EdgeSoundSpeed);
VertSoundSpeed.update();
EdgeSoundSpeed(CoarseEdges[d]) = min( VertSoundSpeed(CoarseEdges[d]),VertSoundSpeed(EdgeNgbr[d]) );
} // end function calcEdgeSoundSpeed
//======================================================================
// Function -- QRelations<Traits>::calcEdgeQTmpMax
//======================================================================
template <class Traits>
void QRelations<Traits>::calcEdgeQTmpMax ( const ScalarField& EdgeQTmpMax,
const ScalarField& EdgeDensity,
const ScalarField& EdgeQPressure )
{
ENTER("QRelations::calcEdgeQTmpMax");
int d = getEdgeDirection(EdgeQTmpMax);
EdgeQTmpMax = EdgeDensity[d] * EdgeQPressure[d];
} // end function calcEdgeQTmpMax
//======================================================================
// Function -- QRelations<Traits>::calcSpokeDVolDt
//======================================================================
template <class Traits>
void QRelations<Traits>::calcSpokeDVolDt ( const ScalarField& SpokeDVolDt,
const VectorField& EdgeVelocity,
const VectorField& SubFaceAreas)
{
ENTER("QRelations::calcSpokeDVolDt");
int d = getFaceDirection(SpokeDVolDt);
EdgeVelocity[d].update();
SubFaceAreas[d].update();
for(int edg=0;edg<nEdgesPerDimension;++edg) {
SpokeDVolDt(Spoke[d][edg]) = dot(EdgeVelocity[d](CellEdge[d][edg]),SubFaceAreas[d](Spoke[d][edg]));
}
} // end function calcSpokeDVolDt
//======================================================================
// Function -- QRelations<Traits>::calcSpokeQSwitch
//======================================================================
template <class Traits>
void QRelations<Traits>::calcSpokeQSwitch ( const ScalarField& SpokeQSwitch,
const ScalarField& SpokeDVolDt,
const ScalarField& Volume,
const ScalarField& Dt )
{
ENTER("QRelations::calcSpokeQSwitch");
int d = getFaceDirection(SpokeQSwitch);
for(int edg=0;edg<nEdgesPerDimension;++edg) {
SpokeQSwitch(Spoke[d][edg]) = (-Dt(CoarseCell) * SpokeDVolDt[d](Spoke[d][edg]))/Volume(CoarseCell);
}
SpokeQSwitch = where(SpokeQSwitch > spokeCutoff, 1.0, 0.0);
} // end function calcSpokeQSwitch
//======================================================================
// Function -- QRelations<Traits>::calcRightLimiterRatio
//======================================================================
template <class Traits>
void QRelations<Traits>::calcRightLimiterRatio ( const ScalarField& RightLimiterRatio,
const VectorField& EdgeLength,
const VectorField& EdgeVelocity)
{
ENTER("QRelations::calcRightLimiterRatio");
int d = getEdgeDirection(RightLimiterRatio);
EdgeLength[d].update();
EdgeVelocity[d].update();
RightLimiterRatio = 1.0;
RightLimiterRatio(RDomain[d]) = dot(EdgeLength[d](RightEdgeNgbr[d]),
EdgeLength[d](RDomain[d])) *
dot(EdgeVelocity[d](RDomain[d]),
EdgeVelocity[d](RDomain[d]));
RightLimiterRatio(RDomain[d]) = where(RightLimiterRatio(RDomain[d]) > numeric_limits<Real>::epsilon(),
(dot(EdgeVelocity[d](RightEdgeNgbr[d]),
EdgeVelocity[d](RDomain[d])) *
dot(EdgeLength[d](RDomain[d]),
EdgeLength[d](RDomain[d])) ) /
RightLimiterRatio(RDomain[d]),
1.0 );
//-------------------------------------------------------------
// if ( d == 0 ) {
// } else if ( d == 1) {
// RightLimiterRatio -= 1.0;
// tecout << "RightLimiterRatio[1] = " << endl;
// tecout << RightLimiterRatio << endl;
// RightLimiterRatio += 1.0;
// }
//-------------------------------------------------------------
} // end function calcRightLimiterRatio
//======================================================================
// Function -- QRelations<Traits>::calcLeftLimiterRatio
//======================================================================
template <class Traits>
void QRelations<Traits>::calcLeftLimiterRatio ( const ScalarField& LeftLimiterRatio,
const VectorField& EdgeLength,
const VectorField& EdgeVelocity)
{
ENTER("QRelations::calcLeftLimiterRatio");
int d = getEdgeDirection(LeftLimiterRatio);
EdgeLength[d].update();
EdgeVelocity[d].update();
LeftLimiterRatio = 1.0;
LeftLimiterRatio(LDomain[d]) = dot(EdgeLength[d](LeftEdgeNgbr[d]),
EdgeLength[d](LDomain[d])) *
dot(EdgeVelocity[d](LDomain[d]),
EdgeVelocity[d](LDomain[d]));
LeftLimiterRatio(LDomain[d]) = where(LeftLimiterRatio(LDomain[d]) > numeric_limits<Real>::epsilon(),
(dot(EdgeVelocity[d](LeftEdgeNgbr[d]),
EdgeVelocity[d](LDomain[d])) *
dot(EdgeLength[d](LDomain[d]),
EdgeLength[d](LDomain[d])) ) /
LeftLimiterRatio(LDomain[d]),
1.0 );
//-------------------------------------------------------------
// if ( d == 0 ) {
// } else if ( d == 1) {
// LeftLimiterRatio -= 1.0;
// tecout << "LeftLimiterRatio[1] = " << endl;
// tecout << LeftLimiterRatio << endl;
// LeftLimiterRatio += 1.0;
// }
//-------------------------------------------------------------
} // end function calcLeftLimiterRatio
//======================================================================
// Function -- QRelations<Traits>::calcEdgePsiLimiter
//======================================================================
template <class Traits>
void QRelations<Traits>::calcEdgePsiLimiter ( const ScalarField& EdgePsiLimiter,
const ScalarField& RightLimiterRatio,
const ScalarField& LeftLimiterRatio,
const VectorField& EdgeLength,
const VectorField& EdgeVelocity,
const ScalarField& Dt)
{
ENTER("QRelations::calcEdgePsiLimiter");
int d = getEdgeDirection(EdgePsiLimiter);
//------------------------------------------------------------------------
// Typically the EdgePsiLimiter is limited between 0 and 1,
// For this algorithm, 0 means do not have any artificial viscosity
// along the edge.
//
// (In Randy Christensen's paper, EdgePsiLimiter is
// PsiLimiter = 1.0 - EdgePsiLimiter, so that when he calculates the
// Q, there is a factor of ( 1.0 - PsiLimiter) in the expression. This
// means his PsiLimiter has a different meaning than the one used here.
// Don M )
//------------------------------------------------------------------------
RightLimiterRatio[d].update();
LeftLimiterRatio[d].update();
EdgeLength[d].update();
EdgeVelocity[d].update();
EdgePsiLimiter = 1.0 - max( 0.0, min( min( 2.0 * RightLimiterRatio[d],
2.0 * LeftLimiterRatio[d]),
min( 0.5 * (RightLimiterRatio[d] +
LeftLimiterRatio[d]),
1.0 ) ) );
//------------------------------------------------------------------------
// The where statement turns off the limiter when:
// 1. The limiter is non zero because of round off
// 2. The edge is degenerate (the vertex points coincide)
// 3. The change in the velocity along the edge is too small for
// anything to happen during the time step.
//------------------------------------------------------------------------
EdgePsiLimiter = where( (EdgePsiLimiter > numeric_limits<Real>::epsilon() &&
dot(EdgeLength[d],EdgeLength[d]) >
numeric_limits<Real>::epsilon() &&
dot(EdgeVelocity[d],EdgeVelocity[d]) * Dt * Dt >
numeric_limits<Real>::epsilon() *
dot(EdgeLength[d],EdgeLength[d]) ),
EdgePsiLimiter, 0.0 );
//---------------------------------------------------------------------------------------------------------
if ( d == 0 ) {
// // tecout << "eps = " << eps <<endl;
// // tecout << "numeric_limits<Real>::epsilon() = " << numeric_limits<Real>::epsilon() << endl;
// tecout << " Dt = " << Dt << endl;
// tecout << "EdgePsiLimiter[0] = " << endl;
// tecout << EdgePsiLimiter << endl;
// tecout << "EdgeVelocity[0] = " << endl;
// tecout << EdgeVelocity[0] << endl;
} else if ( d == 1) {
// // tecout << "RightLimiterRatio[1] = " << endl;
// // tecout << RightLimiterRatio[1] << endl;
// // tecout << "LeftLimiterRatio[1] = " << endl;
// // tecout << LeftLimiterRatio[1] << endl;
// tecout << "EdgePsiLimiter[1] = " << endl;
// tecout << EdgePsiLimiter << endl;
// tecout << "EdgeVelocity[1] = " << endl;
// tecout << EdgeVelocity[1] << endl;
}
//---------------------------------------------------------------------------------------------------------
} // end function calcEdgePsiLimiter
//======================================================================
// Function -- QRelations<Traits>::calcCellQ
//======================================================================
template <class Traits>
void QRelations<Traits>::calcCellQ ( const ScalarField& CellQ,
const ScalarField& EdgeQTmpMax,
const ScalarField& SpokeQSwitch )
{
ENTER("QRelations::calcCellQ");
EdgeQTmpMax.update();
SpokeQSwitch.update();
CellQ = 0.0;
for ( int d = 0; d < Dim; ++d ) {
for ( int edg = 0; edg < nEdgesPerDimension; ++edg ) {
CellQ(CoarseCell) = max( CellQ(CoarseCell), EdgeQTmpMax[d](CellEdge[d][edg]) * SpokeQSwitch[d](Spoke[d][edg]) );
}
}
//----------------------------------------
// tecout << " in calcCellQ: CellQ = " << endl;
// tecout << CellQ << endl;
//----------------------------------------
} // end function calcCellQ
//======================================================================
// Function -- QRelations<Traits>::calcSubForceQ
//======================================================================
template <class Traits>
void QRelations<Traits>::calcSubForceQ ( const VectorField& SubForceQ,
const ScalarField& EdgeQTmpMax,
const VectorField& EdgeVelocity,
const ScalarField& SpokeDVolDt,
const ScalarField& SpokeQSwitch )
{
ENTER("QRelations::calcSubForceQ");
SubForceQ = 0.0;
for ( int d = 0; d < Dim; ++d ) {
Loc<Dim> offset(0); offset[d] = 1;
for ( int edg = 0; edg < nEdgesPerDimension; ++edg ) {
SubForceQ(Spoke[d][edg] - offset) -=(EdgeVelocity[d](CellEdge[d][edg]) *
SpokeDVolDt[d] (Spoke[d][edg]) *
EdgeQTmpMax[d] (CellEdge[d][edg]) *
SpokeQSwitch[d](Spoke[d][edg]) ) /
(dot(EdgeVelocity[d](CellEdge[d][edg]),EdgeVelocity[d](CellEdge[d][edg]))+eps);
SubForceQ(Spoke[d][edg]) +=(EdgeVelocity[d](CellEdge[d][edg]) *
SpokeDVolDt[d] (Spoke[d][edg]) *
EdgeQTmpMax[d] (CellEdge[d][edg]) *
SpokeQSwitch[d](Spoke[d][edg]) ) /
(dot(EdgeVelocity[d](CellEdge[d][edg]),EdgeVelocity[d](CellEdge[d][edg]))+eps);
}
}
//--------------------------------------------
// tecout << " in calcSubForceQ: SubForceQ = " << endl;
// tecout << SubForceQ << endl;
//--------------------------------------------
} // end function calcSubForceQ
//======================================================================
// Function -- QRelations<Traits>::calcSubPressureQMod
//======================================================================
template <class Traits>
void QRelations<Traits>::calcSubPressureQMod ( const ScalarField& SubPressureQMod,
const ScalarField& EdgeQPressure,
const ScalarField& SpokeQSwitch )
{
ENTER("QRelations::calcSubPressureQMod");
SubPressureQMod = 0.0;
for ( int d = 0; d < Dim; ++d ) {
Loc<Dim> offset(0); offset[d] = 1;
for ( int edg = 0; edg < nEdgesPerDimension; ++edg ) {
SubPressureQMod(Spoke[d][edg] - offset) += EdgeQPressure[d](CellEdge[d][edg]) *
SpokeQSwitch[d](Spoke[d][edg]);
SubPressureQMod(Spoke[d][edg]) += EdgeQPressure[d](CellEdge[d][edg]) *
SpokeQSwitch[d](Spoke[d][edg]);
}
}
} // end function calcSubPressureQMod
} // end namespace conejo
#endif // end shroud __conejo_QRelations_t_hh
-------------- next part --------------
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// CCSource File: QRelations
// Author: jcm
// Date: Sat - Nov 18, 2000
// Namespace: conejo
// Framework: Tecolote
// Copyright: Los Alamos National Laboratory
// Full Copyright=$(TECOLOTE_ROOT)/Doc/Copyright
// RCS_VERSION_ID: $Id:
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#if !defined (__MWERKS__)
#include "Demo_src/include/DemoPCH.hh"
#pragma hdrstop
#endif // !__MWERKS__
#include "Demo_src/Model/CompatibleHydro/QRelations.t.hh"
#include "TecFramework_src/MetaTypes/MetaTypes.hh"
#include "TecFramework_src/Foundation/LoadObject.hh"
namespace conejo
{
using namespace TecFramework;
using namespace poomalote;
using namespace Hydrodynamics;
using namespace PhysicsBaseClasses;
using namespace std;
static MetaClass<ITecoloteTraits<QRelations<OneDF<DefaultTraits> >, RelationPkg> >
QRelations1DMeta("QRelations1D", QRelations<OneDF<DefaultTraits> >::MakePersistents());
LoadObjectGroup QRelations1DBase_cc = { &QRelations1DMeta };
} // end namespace conejo
namespace TecFramework
{
using namespace conejo;
LoadObjectGroup QRelations1D_cc = { &QRelations1DBase_cc };
} // end namespace TecFramework
-------------- next part --------------
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// CCSource File: QRelations
// Author: jcm
// Date: Sat - Nov 18, 2000
// Namespace: conejo
// Framework: Tecolote
// Copyright: Los Alamos National Laboratory
// Full Copyright=$(TECOLOTE_ROOT)/Doc/Copyright
// RCS_VERSION_ID: $Id:
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#if !defined (__MWERKS__)
#include "Demo_src/include/DemoPCH.hh"
#pragma hdrstop
#endif // !__MWERKS__
#include "Demo_src/Model/CompatibleHydro/QRelations.t.hh"
#include "TecFramework_src/MetaTypes/MetaTypes.hh"
#include "TecFramework_src/Foundation/LoadObject.hh"
namespace conejo
{
using namespace TecFramework;
using namespace poomalote;
using namespace Hydrodynamics;
using namespace PhysicsBaseClasses;
using namespace std;
static MetaClass<ITecoloteTraits<QRelations<TwoDF<DefaultTraits> >, RelationPkg> >
QRelations2DMeta("QRelations2D", QRelations<TwoDF<DefaultTraits> >::MakePersistents());
LoadObjectGroup QRelations2DBase_cc = { &QRelations2DMeta };
} // end namespace conejo
namespace TecFramework
{
using namespace conejo;
LoadObjectGroup QRelations2D_cc = { &QRelations2DBase_cc };
} // end namespace TecFramework
-------------- next part --------------
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// Header File: QRelations
// Author: jcm
// Date: Sat - Nov 18, 2000
// Namespace: conejo
// Framework: Tecolote
// Copyright: Los Alamos National Laboratory
// Full Copyright=$(TECOLOTE_ROOT)/Doc/Copyright
// RCS_VERSION_ID: $Id:
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#ifndef __conejo_QRelations_hh
#define __conejo_QRelations_hh
#define DEBUG_UPDATERS 1
#include "Hydrodynamics_src/CompatibleHydro/CompatibleRelations.hh"
#include "Pooma2Integration_src/TecMesh/TecMesh.hh"
#include "Pooma/NewFields.h"
namespace conejo
{
using namespace TecFramework;
using namespace Hydrodynamics;
using namespace poomalote;
using namespace std;
//*******************************************************************
// Tecolote Class - QRelations
// ------------------------------------------------------------------
// Summary:
//
// This class holds functions that define Independent/Dependent
// Field Relationships. These Relationships cause Dependent
// Fields to be updated automatically when they are used and
// are out-of-date with respect to their Independent Fields.
//
//
// When declaring Relation functions, include in the parameter
// list the Fields that will participate in the relationships.
// The Dependent Field is listed first (the LField - Left
// of the equals sign) and the Independent Fields are listed
// next (the RFields - Right of the equals sign), followed
// by any scalar values. It is not necessary to include the
// parameter types if they default to the types declared in
// your Model class.
//
// An example of a relation-function declaration is:
//
// GammaLaw( Pressure, Density, IntEnergy);
//
// where Pressure is a function of Density and IntEnergy,
// and the GammaLaw function defines the expressions to update
// the Pressure.
//*******************************************************************
template <class Traits>
class QRelations : public CompatibleRelations<Traits>
{
public:
PERSISTENT_MEMBERS(QRelations)
// Typedefs
FIELD_TYPEDEFS(Traits)
typedef CompatibleRelations<Traits> Parent;
typedef CompatibleRelations<Traits> *pParent;
// Tecolote Constructor
QRelations( DataDirectory* pDataDir, const string& inName);
// Public Member Functions
void createRelations (void);
void calcWeightedVertVol ( const ScalarField& WeightedVertVol,
const ScalarField& SubVertVolume,
const ScalarField& VertVolume );
// void calcEdgeVolume ( const ScalarField& EdgeVolume,
// const ScalarField& WeightedVertVol );
void calcEdgeVelocity ( const VectorField& EdgeVelocity,
const VectorField& OldVelocity);
void calcSpokeDVolDt( const ScalarField& SpokeDVolDt,
const VectorField& EdgeVelocity,
const VectorField& SubFaceAreas);
void calcEdgeDensity ( const ScalarField& EdgeDensity,
const ScalarField& VertDensity);
void calcEdgeGammaConstant( const ScalarField& EdgeGammaConstant,
const ScalarField& VertGammaConstant);
void calcEdgeQPressure ( const ScalarField& EdgeQPressure,
const ScalarField& EdgeGammaConstant,
const ScalarField& EdgeSoundSpeed,
const VectorField& EdgeVelocity,
const ScalarField& EdgePsiLimiter );
void calcEdgeSoundSpeed ( const ScalarField& EdgeSoundSpeed,
const ScalarField& VertSoundSpeed );
void calcEdgeQTmpMax ( const ScalarField& EdgeQTmpMax,
const ScalarField& EdgeDensity,
const ScalarField& EdgeQPressure );
void calcRightLimiterRatio( const ScalarField& RightLimiterRatio,
const VectorField& EdgeDeltaLength,
const VectorField& EdgeVelocity);
void calcLeftLimiterRatio ( const ScalarField& LeftLimiterRatio,
const VectorField& EdgeDeltaLength,
const VectorField& EdgeVelocity);
void calcEdgePsiLimiter ( const ScalarField& EdgePsiLimiter,
const ScalarField& RightLimiterRatio,
const ScalarField& LeftLimiterRatio,
const VectorField& EdgeDeltaLength,
const VectorField& EdgeVelocity,
const ScalarField& Dt);
void calcSpokeQSwitch ( const ScalarField& SpokeQSwitch,
const ScalarField& SpokeDVolDt,
const ScalarField& Volume,
const ScalarField& Dt);
void calcCellQ ( const ScalarField& CellQ,
const ScalarField& EdgeQTmpMax,
const ScalarField& SpokeQSwitch );
void calcSubForceQ ( const VectorField& SubForceQ,
const ScalarField& EdgeQTmpMax,
const VectorField& EdgeVelocity,
const ScalarField& SpokeDVolDt,
const ScalarField& SpokeQSwitch );
void calcSubPressureQMod( const ScalarField& SubPressureQMod,
const ScalarField& EdgeQPressure,
const ScalarField& SpokeQSwitch );
//..............................................................
// The 6 C++ default methods
//
public:
virtual ~QRelations();
private:
QRelations();
// QRelations(const QRelations& c);
QRelations& operator=(const QRelations& c);
// QRelations* operator&();
// const QRelations* operator&() const;
//..............................................................
private:
// Private Member Functions
//----------------
// Member Data
//----------------
DataDirectory& Old;
Real linearQ;
Real quadQ;
Real gammaConst;
vector<Range<Dim> > RDomain;
vector<Range<Dim> > LDomain;
vector<Range<Dim> > RightEdgeNgbr;
vector<Range<Dim> > LeftEdgeNgbr;
vector<Range<Dim> > CoarseEdges;
vector<Range<Dim> > FineEdges;
vector<Range<Dim> > EdgeNgbr;
vector<Range<Dim> > UpperVert;
vector<Range<Dim> > LowerVert;
}; // end class QRelations
} // end namespace conejo
#ifndef TEC_INCLUDE_T_HH_FILE
#include "Demo_src/Model/CompatibleHydro/QRelations.t.hh"
#endif
#endif // end shroud __conejo_QRelations_hh
More information about the pooma-dev
mailing list