Commit 5d24f091 authored by mike turner's avatar mike turner
Browse files

clean rest

parent 0b799111
......@@ -227,7 +227,7 @@ int main(int argc, char *argv[])
}
}
char *vars[3] = { "u", "v", "w" };
string vars = "uvw";
for (int i = 0; i < dim; ++i)
{
......
......@@ -14,7 +14,7 @@ NekDouble Tri_sol(NekDouble x, NekDouble y, int order1, int order2);
int main(int argc, char *argv[])
{
int i, j;
int i;
int order, nq1, nq2;
LibUtilities::PointsType Qtype1, Qtype2;
LibUtilities::BasisType btype1, btype2;
......
///////////////////////////////////////////////////////////////////////////////
//
// File: SubSteppingExtrapolate.cpp
//
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
......@@ -76,130 +76,130 @@ namespace Nektar
SubSteppingExtrapolate::~SubSteppingExtrapolate()
{
}
void SubSteppingExtrapolate::v_EvaluatePressureBCs(const Array<OneD, const Array<OneD, NekDouble> > &fields, const Array<OneD, const Array<OneD, NekDouble> > &N, NekDouble kinvis)
{
ASSERTL0(false,"This method should not be called by Substepping routine");
}
void SubSteppingExtrapolate::v_SubSteppingTimeIntegration(
int intMethod,
const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)
{
int i;
// Set to 1 for first step and it will then be increased in
// time advance routines
switch(intMethod)
{
case LibUtilities::eBackwardEuler:
case LibUtilities::eBDFImplicitOrder1:
case LibUtilities::eBDFImplicitOrder1:
{
std::string vSubStepIntScheme = "ForwardEuler";
if(m_session->DefinesSolverInfo("SubStepIntScheme"))
{
vSubStepIntScheme = m_session->GetSolverInfo("SubStepIntScheme");
}
m_subStepIntegrationScheme = LibUtilities::GetTimeIntegrationWrapperFactory().CreateInstance(vSubStepIntScheme);
int nvel = m_velocity.num_elements();
// Fields for linear interpolation
m_previousVelFields = Array<OneD, Array<OneD, NekDouble> >(2*nvel);
m_previousVelFields = Array<OneD, Array<OneD, NekDouble> >(2*nvel);
int ntotpts = m_fields[0]->GetTotPoints();
m_previousVelFields[0] = Array<OneD, NekDouble>(2*nvel*ntotpts);
for(i = 1; i < 2*nvel; ++i)
{
m_previousVelFields[i] = m_previousVelFields[i-1] + ntotpts;
m_previousVelFields[i] = m_previousVelFields[i-1] + ntotpts;
}
}
break;
case LibUtilities::eBDFImplicitOrder2:
{
std::string vSubStepIntScheme = "RungeKutta2_ImprovedEuler";
if(m_session->DefinesSolverInfo("SubStepIntScheme"))
{
vSubStepIntScheme = m_session->GetSolverInfo("SubStepIntScheme");
}
m_subStepIntegrationScheme = LibUtilities::GetTimeIntegrationWrapperFactory().CreateInstance(vSubStepIntScheme);
int nvel = m_velocity.num_elements();
// Fields for quadratic interpolation
m_previousVelFields = Array<OneD, Array<OneD, NekDouble> >(3*nvel);
int ntotpts = m_fields[0]->GetTotPoints();
m_previousVelFields[0] = Array<OneD, NekDouble>(3*nvel*ntotpts);
for(i = 1; i < 3*nvel; ++i)
{
m_previousVelFields[i] = m_previousVelFields[i-1] + ntotpts;
m_previousVelFields[i] = m_previousVelFields[i-1] + ntotpts;
}
}
break;
default:
ASSERTL0(0,"Integration method not suitable: Options include BackwardEuler or BDFImplicitOrder1");
break;
}
m_intSteps = IntegrationScheme->GetIntegrationSteps();
// set explicit time-integration class operators
m_subStepIntegrationOps.DefineOdeRhs(&SubSteppingExtrapolate::SubStepAdvection, this);
m_subStepIntegrationOps.DefineProjection(&SubSteppingExtrapolate::SubStepProjection, this);
}
/**
/**
* Explicit Advection terms used by SubStepAdvance time integration
*/
void SubSteppingExtrapolate::SubStepAdvection(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble time)
{
int i;
int nVariables = inarray.num_elements();
int nQuadraturePts = inarray[0].num_elements();
/// Get the number of coefficients
int ncoeffs = m_fields[0]->GetNcoeffs();
/// Define an auxiliary variable to compute the RHS
int ncoeffs = m_fields[0]->GetNcoeffs();
/// Define an auxiliary variable to compute the RHS
Array<OneD, Array<OneD, NekDouble> > WeakAdv(nVariables);
WeakAdv[0] = Array<OneD, NekDouble> (ncoeffs*nVariables);
for(i = 1; i < nVariables; ++i)
{
WeakAdv[i] = WeakAdv[i-1] + ncoeffs;
}
Array<OneD, Array<OneD, NekDouble> > Velfields(m_velocity.num_elements());
Velfields[0] = Array<OneD, NekDouble> (nQuadraturePts*m_velocity.num_elements());
for(i = 1; i < m_velocity.num_elements(); ++i)
{
Velfields[i] = Velfields[i-1] + nQuadraturePts;
Velfields[i] = Velfields[i-1] + nQuadraturePts;
}
SubStepExtrapolateField(fmod(time,m_timestep), Velfields);
m_advObject->Advect(m_velocity.num_elements(), m_fields, Velfields, inarray, outarray, time);
for(i = 0; i < nVariables; ++i)
{
m_fields[i]->IProductWRTBase(outarray[i],WeakAdv[i]);
m_fields[i]->IProductWRTBase(outarray[i],WeakAdv[i]);
// negation requried due to sign of DoAdvection term to be consistent
Vmath::Neg(ncoeffs, WeakAdv[i], 1);
}
AddAdvectionPenaltyFlux(Velfields, inarray, WeakAdv);
/// Operations to compute the RHS
for(i = 0; i < nVariables; ++i)
{
......@@ -208,18 +208,18 @@ namespace Nektar
/// Multiply the flux by the inverse of the mass matrix
m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
/// Store in outarray the physical values of the RHS
m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
}
}
/**
/**
* Projection used by SubStepAdvance time integration
*/
void SubSteppingExtrapolate::SubStepProjection(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble time)
{
ASSERTL1(inarray.num_elements() == outarray.num_elements(),"Inarray and outarray of different sizes ");
......@@ -231,17 +231,17 @@ namespace Nektar
}
/**
*
/**
*
*/
void SubSteppingExtrapolate::v_SubStepSetPressureBCs(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const NekDouble Aii_Dt,
NekDouble kinvis)
{
int nConvectiveFields =m_fields.num_elements()-1;
//int nConvectiveFields =m_fields.num_elements()-1;
Array<OneD, Array<OneD, NekDouble> > nullvelfields;
m_pressureCalls++;
// Calculate non-linear and viscous BCs at current level and
......@@ -250,20 +250,20 @@ namespace Nektar
// Extrapolate to m_pressureHBCs to n+1
ExtrapolateArray(m_pressureHBCs);
// Add (phi,Du/Dt) term to m_presureHBC
// Add (phi,Du/Dt) term to m_presureHBC
AddDuDt();
// Copy m_pressureHBCs to m_PbndExp
CopyPressureHBCsToPbndExp();
CopyPressureHBCsToPbndExp();
// Evaluate High order outflow conditiosn if required.
// Evaluate High order outflow conditiosn if required.
CalcOutflowBCs(inarray, kinvis);
}
/**
*
/**
*
*/
void SubSteppingExtrapolate::v_SubStepSaveFields(const int nstep)
{
......@@ -271,61 +271,61 @@ namespace Nektar
int nvel = m_velocity.num_elements();
int npts = m_fields[0]->GetTotPoints();
// rotate fields
int nblocks = m_previousVelFields.num_elements()/nvel;
Array<OneD, NekDouble> save;
// rotate fields
int nblocks = m_previousVelFields.num_elements()/nvel;
Array<OneD, NekDouble> save;
// rotate storage space
for(n = 0; n < nvel; ++n)
{
save = m_previousVelFields[(nblocks-1)*nvel+n];
for(i = nblocks-1; i > 0; --i)
{
m_previousVelFields[i*nvel+n] = m_previousVelFields[(i-1)*nvel+n];
}
m_previousVelFields[n] = save;
m_previousVelFields[n] = save;
}
// Put previous field
// Put previous field
for(i = 0; i < nvel; ++i)
{
m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
m_fields[m_velocity[i]]->UpdatePhys());
Vmath::Vcopy(npts,m_fields[m_velocity[i]]->GetPhys(),1,
m_previousVelFields[i],1);
m_previousVelFields[i],1);
}
if(nstep == 0)// initialise all levels with first field
if(nstep == 0)// initialise all levels with first field
{
for(n = 0; n < nvel; ++n)
{
for(i = 1; i < nblocks; ++i)
{
Vmath::Vcopy(npts,m_fields[m_velocity[n]]->GetPhys(),1,
m_previousVelFields[i*nvel+n],1);
m_previousVelFields[i*nvel+n],1);
}
}
}
}
/**
*
/**
*
*/
void SubSteppingExtrapolate::v_SubStepAdvance(
const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln,
int nstep,
const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln,
int nstep,
NekDouble time)
{
int n;
int nsubsteps;
NekDouble dt;
NekDouble dt;
Array<OneD, Array<OneD, NekDouble> > fields;
static int ncalls = 1;
int nint = min(ncalls++, m_intSteps);
......@@ -335,17 +335,17 @@ namespace Nektar
// Get the proper time step with CFL control
dt = GetSubstepTimeStep();
nsubsteps = (m_timestep > dt)? ((int)(m_timestep/dt)+1):1;
nsubsteps = (m_timestep > dt)? ((int)(m_timestep/dt)+1):1;
nsubsteps = max(m_minsubsteps, nsubsteps);
ASSERTL0(nsubsteps < m_maxsubsteps,"Number of substeps has exceeded maximum");
dt = m_timestep/nsubsteps;
if (m_infosteps && !((nstep+1)%m_infosteps) && m_comm->GetRank() == 0)
{
cout << "Sub-integrating using "<< nsubsteps
<< " steps over Dt = " << m_timestep
cout << "Sub-integrating using "<< nsubsteps
<< " steps over Dt = " << m_timestep
<< " (SubStep CFL=" << m_cflSafetyFactor << ")"<< endl;
}
......@@ -353,60 +353,60 @@ namespace Nektar
{
// We need to update the fields held by the m_integrationSoln
fields = integrationSoln->UpdateSolutionVector()[m];
// Initialise NS solver which is set up to use a GLM method
// with calls to EvaluateAdvection_SetPressureBCs and
// SolveUnsteadyStokesSystem
LibUtilities::TimeIntegrationSolutionSharedPtr
LibUtilities::TimeIntegrationSolutionSharedPtr
SubIntegrationSoln = m_subStepIntegrationScheme->
InitializeScheme(dt, fields, time, m_subStepIntegrationOps);
for(n = 0; n < nsubsteps; ++n)
{
fields = m_subStepIntegrationScheme->TimeIntegrate(n, dt, SubIntegrationSoln,
m_subStepIntegrationOps);
}
// set up HBC m_acceleration field for Pressure BCs
IProductNormVelocityOnHBC(fields,m_iprodnormvel[m]);
// Reset time integrated solution in m_integrationSoln
// Reset time integrated solution in m_integrationSoln
integrationSoln->SetSolVector(m,fields);
}
}
/**
*
/**
*
*/
NekDouble SubSteppingExtrapolate::GetSubstepTimeStep()
{
int n_element = m_fields[0]->GetExpSize();
{
int n_element = m_fields[0]->GetExpSize();
const Array<OneD, int> ExpOrder=m_fields[0]->EvalBasisNumModesMaxPerExp();
const NekDouble cLambda = 0.2; // Spencer book pag. 317
Array<OneD, NekDouble> tstep (n_element, 0.0);
Array<OneD, NekDouble> stdVelocity(n_element, 0.0);
Array<OneD, Array<OneD, NekDouble> > velfields(m_velocity.num_elements());
for(int i = 0; i < m_velocity.num_elements(); ++i)
{
velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
}
}
stdVelocity = GetMaxStdVelocity(velfields);
for(int el = 0; el < n_element; ++el)
{
tstep[el] = m_cflSafetyFactor /
(stdVelocity[el] * cLambda *
tstep[el] = m_cflSafetyFactor /
(stdVelocity[el] * cLambda *
(ExpOrder[el]-1) * (ExpOrder[el]-1));
}
NekDouble TimeStep = Vmath::Vmin(n_element, tstep, 1);
m_comm->AllReduce(TimeStep,LibUtilities::ReduceMin);
m_comm->AllReduce(TimeStep,LibUtilities::ReduceMin);
return TimeStep;
}
......@@ -414,47 +414,47 @@ namespace Nektar
void SubSteppingExtrapolate::AddAdvectionPenaltyFlux(
const Array<OneD, const Array<OneD, NekDouble> > &velfield,
const Array<OneD, const Array<OneD, NekDouble> > &physfield,
const Array<OneD, const Array<OneD, NekDouble> > &velfield,
const Array<OneD, const Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &Outarray)
{
ASSERTL1(
physfield.num_elements() == Outarray.num_elements(),
"Physfield and outarray are of different dimensions");
int i;
/// Number of trace points
int nTracePts = m_fields[0]->GetTrace()->GetNpoints();
/// Number of spatial dimensions
int nDimensions = m_bnd_dim;
/// Forward state array
Array<OneD, NekDouble> Fwd(3*nTracePts);
/// Backward state array
Array<OneD, NekDouble> Bwd = Fwd + nTracePts;
/// upwind numerical flux state array
Array<OneD, NekDouble> numflux = Bwd + nTracePts;
/// Normal velocity array
Array<OneD, NekDouble> Vn (nTracePts, 0.0);
// Extract velocity field along the trace space and multiply by trace normals
for(i = 0; i < nDimensions; ++i)
{
m_fields[0]->ExtractTracePhys(velfield[i], Fwd);
Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Fwd, 1, Vn, 1, Vn, 1);
}
for(i = 0; i < physfield.num_elements(); ++i)
{
/// Extract forwards/backwards trace spaces
/// Note: Needs to have correct i value to get boundary conditions
m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
/// Upwind between elements
m_fields[0]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux);
......@@ -470,13 +470,13 @@ namespace Nektar
m_fields[0]->AddFwdBwdTraceIntegral(Fwd,Bwd,Outarray[i]);
}
}
/**
/**
* Extrapolate field using equally time spaced field un,un-1,un-2, (at
* dt intervals) to time n+t at order Ord
* dt intervals) to time n+t at order Ord
*/
void SubSteppingExtrapolate::SubStepExtrapolateField(
NekDouble toff,
NekDouble toff,
Array< OneD, Array<OneD, NekDouble> > &ExtVel)
{
int npts = m_fields[0]->GetTotPoints();
......@@ -485,10 +485,10 @@ namespace Nektar
Array<OneD, NekDouble> l(4);
int ord = m_intSteps;
// calculate Lagrange interpolants
Vmath::Fill(4,1.0,l,1);
for(i = 0; i <= ord; ++i)
{
for(j = 0; j <= ord; ++j)
......@@ -504,7 +504,7 @@ namespace Nektar
for(i = 0; i < nvel; ++i)
{
Vmath::Smul(npts,l[0],m_previousVelFields[i],1,ExtVel[i],1);
for(j = 1; j <= ord; ++j)
{
Blas::Daxpy(npts,l[j],m_previousVelFields[j*nvel+i],1,
......@@ -512,14 +512,14 @@ namespace Nektar
}
}
}
/**
*
/**
*
*/
void SubSteppingExtrapolate::v_MountHOPBCs(
int HBCdata,
NekDouble kinvis,
Array<OneD, NekDouble> &Q,
int HBCdata,
NekDouble kinvis,
Array<OneD, NekDouble> &Q,
Array<OneD, const NekDouble> &Advection)
{
Vmath::Smul(HBCdata,-kinvis,Q,1,Q,1);
......@@ -531,4 +531,3 @@ namespace Nektar
}
}
......@@ -61,14 +61,14 @@ namespace Nektar
SubSteppingExtrapolateWeakPressure::~SubSteppingExtrapolateWeakPressure()
{
}
void SubSteppingExtrapolateWeakPressure::v_SubStepSetPressureBCs(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const NekDouble Aii_Dt,