Commit 288508fa authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'feature/APE-merge' of /opt/gitlab/repositories/nektar

parents 0181fc0f aca8cfce
......@@ -169,12 +169,12 @@ namespace Nektar
Array<OneD, const NekDouble> &from,
Array<OneD, const NekDouble> &to,
NekDouble *mat);
void rotateToNormal (
SOLVER_UTILS_EXPORT void rotateToNormal (
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const Array<OneD, const Array<OneD, NekDouble> > &normals,
const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
Array<OneD, Array<OneD, NekDouble> > &outarray);
void rotateFromNormal(
SOLVER_UTILS_EXPORT void rotateFromNormal(
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
const Array<OneD, const Array<OneD, NekDouble> > &normals,
const Array<OneD, const Array<OneD, NekDouble> > &vecLocs,
......
......@@ -34,7 +34,7 @@
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/EquationSystem.h>
#include <SolverUtils/Driver.h>
#include <LibUtilities/BasicUtils/SessionReader.h>
using namespace Nektar;
......@@ -42,53 +42,33 @@ using namespace Nektar::SolverUtils;
int main(int argc, char *argv[])
{
// Create session reader.
LibUtilities::SessionReaderSharedPtr session;
session = LibUtilities::SessionReader::CreateInstance(argc, argv);
string vDriverModule;
DriverSharedPtr drv;
time_t starttime, endtime;
NekDouble CPUtime;
EquationSystemSharedPtr equ;
// Record start time.
time(&starttime);
// Create instance of module to solve the equation specified in the session.
try
{
equ = GetEquationSystemFactory().CreateInstance(
session->GetSolverInfo("EQTYPE"), session);
}
catch (int e)
{
ASSERTL0(e == -1, "No such solver class defined.");
}
// Create session reader.
session = LibUtilities::SessionReader::CreateInstance(argc, argv);
// Print a summary of solver and problem parameters and initialise the
// solver.
equ->PrintSummary(cout);
equ->DoInitialise();
// Create driver
session->LoadSolverInfo("Driver", vDriverModule, "Standard");
drv = GetDriverFactory().CreateInstance(vDriverModule, session);
// Solve the problem.
equ->DoSolve();
// Execute driver
drv->Execute();
// Record end time.
time(&endtime);
CPUtime = (1.0/60.0/60.0)*difftime(endtime,starttime);
// Write output to .fld file
equ->Output();
// Evaluate and output computation time and solution accuracy.
// The specific format of the error output is essential for the
// regression tests to work.
cout << "-------------------------------------------" << endl;
cout << "Total Computation Time = " << CPUtime << " hr." << endl;
for(int i = 0; i < equ->GetNvariables(); ++i)
// Finalise session
session->Finalise();
}
catch (const std::runtime_error& e)
{
return 1;
}
catch (const std::string& eStr)
{
cout << "L 2 error (variable " << equ->GetVariable(i) << "): " << equ->L2Error(i,true) << endl;
cout << "L inf error (variable " << equ->GetVariable(i) << "): " << equ->LinfError(i) << endl;
cout << "Error: " << eStr << endl;
}
session->Finalise();
return 0;
}
......@@ -5,12 +5,23 @@ CMAKE_DEPENDENT_OPTION(NEKTAR_SOLVER_APE
IF( NEKTAR_SOLVER_APE )
SET(APESolverSource
./APESolver.cpp
./EquationSystems/APESystem.cpp
./EquationSystems/APE.cpp)
./EquationSystems/APE.cpp
./RiemannSolvers/APEUpwindSolver.cpp)
ADD_SOLVER_EXECUTABLE(APESolver solvers
${APESolverSource})
ADD_NEKTAR_TEST(Channel)
ADD_NEKTAR_TEST(Pulse)
# TODO: We have a bug somewhere in the 1D boundary conditions. Therefore 1D
# problems are currently disabled. This should get fixed in the future. (See APE.cpp:64)
# ADD_NEKTAR_TEST(APE_1DPulseSource_FRDG_MODIFIED)
# ADD_NEKTAR_TEST(APE_1DPulseSource_Galerkin_MODIFIED)
# ADD_NEKTAR_TEST(APE_1DPulseWall_FRDG_MODIFIED)
ADD_NEKTAR_TEST(APE_2DChannel_FRDG_MODIFIED)
ADD_NEKTAR_TEST(APE_2DPulseAdv_FRDG_MODIFIED)
ADD_NEKTAR_TEST(APE_2DPulseAdv_Galerkin_MODIFIED)
ADD_NEKTAR_TEST(APE_2DPulseWall_FRDG_MODIFIED)
# ADD_NEKTAR_TEST(APE_2DPulseWall_Galerkin_MODIFIED)
ADD_NEKTAR_TEST(APE_3DPulse_FRDG_MODIFIED)
ADD_NEKTAR_TEST(APE_3DPulseWall_FRDG_MODIFIED)
ENDIF( NEKTAR_SOLVER_APE )
This diff is collapsed.
......@@ -36,142 +36,98 @@
#ifndef NEKTAR_SOLVERS_APESOLVER_EQUATIONSYSTEMS_APE_H
#define NEKTAR_SOLVERS_APESOLVER_EQUATIONSYSTEMS_APE_H
#include <SolverUtils/EquationSystem.h>
#include <APESolver/EquationSystems/APESystem.h>
#include <SolverUtils/UnsteadySystem.h>
#include <SolverUtils/RiemannSolvers/RiemannSolver.h>
using namespace Nektar::SolverUtils;
namespace Nektar
{
/**
*
*
**/
class APE : public APESystem
{
public:
friend class MemoryManager<APE>;
/// Creates an instance of this class
static EquationSystemSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession)
{
EquationSystemSharedPtr p = MemoryManager<APE>::AllocateSharedPtr(pSession);
p->InitObject();
return p;
}
/// Name of class
static std::string className;
virtual ~APE();
protected:
APE(const LibUtilities::SessionReaderSharedPtr& pSession);
virtual void v_InitObject();
Array<OneD, Array<OneD, NekDouble> > basefield;
void DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble time);
void DoOdeProjection(const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble time);
virtual void v_GetFluxVector(const int i, Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &flux)
{
switch(m_expdim)
{
case 1:
ASSERTL0(false,"1D not implemented for Acoustic perturbation equations");
break;
case 2:
GetFluxVector2D(i,physfield,flux);
break;
case 3:
ASSERTL0(false,"3D not implemented for Acoustic perturbation equations");
break;
default:
ASSERTL0(false,"Illegal dimension");
}
}
virtual void v_GenerateSummary(SolverUtils::SummaryList& s);
virtual void v_NumericalFlux(Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numfluxX,
Array<OneD, Array<OneD, NekDouble> > &numfluxY)
{
switch(m_expdim)
{
case 1:
ASSERTL0(false,"1D not implemented for Acoustic perturbation equations");
break;
case 2:
NumericalFlux2D(physfield,numfluxX,numfluxY);
break;
case 3:
ASSERTL0(false,"3D not implemented for Acoustic perturbation equations");
break;
default:
ASSERTL0(false,"Illegal dimension");
}
}
virtual void v_PrimitiveToConservative( );
virtual void v_ConservativeToPrimitive( );
void InitialiseBaseFlowAnalytical(Array<OneD, Array<OneD, NekDouble> > &base, const NekDouble time);
void GetSource(Array<OneD, NekDouble> &source, const NekDouble time);
void AddSource(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray);
class APE : public UnsteadySystem
{
public:
private:
friend class MemoryManager<APE>;
/// Creates an instance of this class
static EquationSystemSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession)
{
EquationSystemSharedPtr p = MemoryManager<APE>::AllocateSharedPtr(pSession);
p->InitObject();
return p;
}
/// Name of class
static std::string className;
/// Destructor
virtual ~APE();
protected:
SolverUtils::RiemannSolverSharedPtr m_riemannSolver;
Array<OneD, Array<OneD, NekDouble> > m_traceBasefield;
Array<OneD, Array<OneD, NekDouble> > m_vecLocs;
/// Constant incompressible density (APE)
NekDouble m_Rho0;
/// Isentropic coefficient, Ratio of specific heats (APE)
NekDouble m_gamma;
Array<OneD, Array<OneD, NekDouble> > m_basefield;
std::vector<std::string> m_basefield_names;
/// Initialises UnsteadySystem class members.
APE(const LibUtilities::SessionReaderSharedPtr& pSession);
virtual void v_InitObject();
void GetFluxVector1D(const int i, Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &flux);
void GetFluxVector2D(const int i, const Array<OneD, const Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &flux);
/// Sets up initial conditions.
virtual void v_DoInitialise();
void NumericalFlux1D(Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numfluxX);
void NumericalFlux2D(Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numfluxX,
Array<OneD, Array<OneD, NekDouble> > &numfluxY);
void DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble time);
void DoOdeProjection(const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble time);
void RiemannSolverUpwind(NekDouble hL,NekDouble uL,NekDouble vL,NekDouble hR,NekDouble uR, NekDouble vR, NekDouble P0, NekDouble U0, NekDouble V0,
NekDouble &hflux, NekDouble &huflux,NekDouble &hvflux );
virtual void v_GetFluxVector(const int i,
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &flux);
void SetBoundaryConditions(Array<OneD, Array<OneD, NekDouble> > &physarray, NekDouble time);
/// Evaulate flux = m_fields*ivel for i th component of Vu for direction j
virtual void v_GetFluxVector(const int i, const int j,
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &flux);
void WallBoundary1D(int bcRegion, Array<OneD, Array<OneD, NekDouble> > &physarray);
///
virtual void v_NumericalFlux(Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numflux);
void WallBoundary2D(int bcRegion, int cnt, Array<OneD, Array<OneD, NekDouble> > &physarray);
virtual void v_NumericalFlux(Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numfluxX,
Array<OneD, Array<OneD, NekDouble> > &numfluxY);
void AddSource(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray);
const Array<OneD, const Array<OneD, NekDouble> > &GetNormals();
const Array<OneD, const Array<OneD, NekDouble> > &GetVecLocs();
const Array<OneD, const Array<OneD, NekDouble> > &GetBasefield();
NekDouble GetGamma();
NekDouble GetRho();
private:
void ConservativeToPrimitive(const Array<OneD, const Array<OneD, NekDouble> >&physin,
Array<OneD, Array<OneD, NekDouble> >&physout);
void PrimitiveToConservative(const Array<OneD, const Array<OneD, NekDouble> >&physin,
Array<OneD, Array<OneD, NekDouble> >&physout);
void SetBoundaryConditions(Array<OneD, Array<OneD, NekDouble> > &physarray, NekDouble time);
};
void WallBC(int bcRegion, int cnt, Array<OneD, Array<OneD, NekDouble> > &physarray);
};
}
#endif
......
///////////////////////////////////////////////////////////////////////////////
//
// File APESystem.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Generic timestepping for APE solvers
//
///////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <LibUtilities/TimeIntegration/TimeIntegrationScheme.h>
#include <APESolver/EquationSystems/APESystem.h>
namespace Nektar
{
/**
* @class APESystem
*
* Provides the underlying timestepping framework for APE flow solvers
* including the general timestepping routines. This class is not intended
* to be directly instantiated, but rather is a base class on which to
* define APE solvers, e.g. SWE, Boussinesq, linear and nonlinear versions.
*
* For details on implementing unsteady solvers see
* \ref sectionADRSolverModuleImplementation
*/
/**
* Processes SolverInfo parameters from the session file and sets up
* timestepping-specific code.
* @param pSession Session object to read parameters from.
*/
APESystem::APESystem(
const LibUtilities::SessionReaderSharedPtr& pSession)
: EquationSystem(pSession)
{
}
void APESystem::v_InitObject()
{
EquationSystem::v_InitObject();
// Load SolverInfo parameters
m_session->MatchSolverInfo("DIFFUSIONADVANCEMENT","Explicit",
m_explicitDiffusion,true);
m_session->MatchSolverInfo("ADVECTIONADVANCEMENT","Explicit",
m_explicitAdvection,true);
// Determine TimeIntegrationMethod to use.
ASSERTL0(m_session->DefinesSolverInfo("TIMEINTEGRATIONMETHOD"),
"No TIMEINTEGRATIONMETHOD defined in session.");
int i;
for (i = 0; i < (int)LibUtilities::SIZE_TimeIntegrationMethod; ++i)
{
bool match;
m_session->MatchSolverInfo("TIMEINTEGRATIONMETHOD",
LibUtilities::TimeIntegrationMethodMap[i], match, false);
if (match)
{
m_timeIntMethod = (LibUtilities::TimeIntegrationMethod) i;
break;
}
}
ASSERTL0(i != (int) LibUtilities::SIZE_TimeIntegrationMethod,
"Invalid time integration type.");
// if discontinuous determine numerical flux to use
if (m_projectionType == MultiRegions::eDiscontinuous)
{
ASSERTL0(m_session->DefinesSolverInfo("UPWINDTYPE"),
"No UPWINDTYPE defined in session.");
int i;
for (i = 0; i < (int)SIZE_UpwindType; ++i)
{
bool match;
m_session->MatchSolverInfo("UPWINDTYPE",
UpwindTypeMap[i], match, false);
if (match)
{
m_upwindType = (UpwindType) i;
break;
}
}
ASSERTL0(i != (int) SIZE_UpwindType,
"Invalid upwind type.");
}
else
{
m_upwindType = (UpwindType) 0;
}
// Load generic input parameters
m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
// Load acceleration of gravity
m_session->LoadParameter("Gravity", m_g, 9.81);
// Load constant incompressible density (APE)
m_session->LoadParameter("Rho0", m_Rho0, 1.204);
// Load isentropic coefficient, Ratio of specific heats (APE)
m_session->LoadParameter("Gamma", m_gamma, 1.4);
// input/output in primitive variables
m_primitive = true;
}
/**
*
*/
APESystem::~APESystem()
{
}
/**
* Initialises the time integration scheme (as specified in the session
* file), and perform the time integration.
*/
void APESystem::v_DoSolve()
{
int i,n,nchk = 0;
int nq = m_fields[0]->GetTotPoints();
int nvariables = m_fields.num_elements();
// Set up wrapper to fields data storage.
Array<OneD, Array<OneD, NekDouble> > fields(nvariables);
Array<OneD, Array<OneD, NekDouble> > tmp(nvariables);
for(i = 0; i < nvariables; ++i)
{
m_fields[i]->SetPhysState(false);
fields[i] = m_fields[i]->UpdatePhys();
}
/* Declare an array of TimeIntegrationSchemes For multi-stage
methods, this array will have just one entry containing the
actual multi-stage method...
For multi-steps method, this can have multiple entries
- the first scheme will used for the first timestep (this
is an initialization scheme)
- the second scheme will used for the second timestep
(this is an initialization scheme)
- ...
- the last scheme will be used for all other time-steps
(this will be the actual scheme)*/
Array<OneD, LibUtilities::TimeIntegrationSchemeSharedPtr> IntScheme;
LibUtilities::TimeIntegrationSolutionSharedPtr u;
int numMultiSteps;
switch(m_timeIntMethod)
{
case LibUtilities::eForwardEuler:
case LibUtilities::eClassicalRungeKutta4:
{
numMultiSteps = 1;
IntScheme = Array<OneD, LibUtilities::TimeIntegrationSchemeSharedPtr>(numMultiSteps);
LibUtilities::TimeIntegrationSchemeKey IntKey(m_timeIntMethod);
IntScheme[0] = LibUtilities::TimeIntegrationSchemeManager()[IntKey];
u = IntScheme[0]->InitializeScheme(m_timestep,fields,m_time,m_ode);
break;
}
case LibUtilities::eAdamsBashforthOrder2:
{
numMultiSteps = 2;
IntScheme = Array<OneD, LibUtilities::TimeIntegrationSchemeSharedPtr>(numMultiSteps);
// Used in the first time step to initalize the scheme
LibUtilities::TimeIntegrationSchemeKey IntKey0(LibUtilities::eClassicalRungeKutta4);
// Used for all other time steps
LibUtilities::TimeIntegrationSchemeKey IntKey1(m_timeIntMethod);
IntScheme[0] = LibUtilities::TimeIntegrationSchemeManager()[IntKey0];
IntScheme[1] = LibUtilities::TimeIntegrationSchemeManager()[IntKey1];
// Initialise the scheme for the actual time integration scheme
u = IntScheme[1]->InitializeScheme(m_timestep,fields,m_time,m_ode);
break;
}
default:
{
ASSERTL0(false,"populate switch statement for integration scheme");
}
}
std::string outname = m_session->GetFilename() + ".his";
std::ofstream hisFile (outname.c_str());
// Perform integration in time.
for(n = 0; n < m_steps; ++n)
{
// Integrate over timestep.
if( n < numMultiSteps-1)
{
// Use initialisation schemes if time step is less than the
// number of steps in the scheme.
fields = IntScheme[n]->TimeIntegrate(m_timestep,u,m_ode);
}
else
{
fields = IntScheme[numMultiSteps-1]->TimeIntegrate(m_timestep,u,m_ode);
}
// Increment time.
m_time += m_timestep;