Commit 88080de7 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Fixed minor formatting issues and put back Aliasing utility.

parent 91daf310
......@@ -78,7 +78,7 @@ namespace Nektar
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
int nDim = fields[0]->GetCoordim(0);
......
/*
* AdvectionSystem.cpp
*
* Created on: 28 May 2014
* Author: cc
*/
///////////////////////////////////////////////////////////////////////////////
//
// File: AdvectionSystem.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: Base class for advection-based equation systems.
//
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/AdvectionSystem.h>
......
......@@ -43,39 +43,39 @@ namespace Nektar
{
string DriverSteadyState::className = GetDriverFactory().RegisterCreatorFunction("SteadyState", DriverSteadyState::create);
string DriverSteadyState::driverLookupId = LibUtilities::SessionReader::RegisterEnumValue("Driver","SteadyState",0);
/**
*
*
*/
DriverSteadyState::DriverSteadyState(const LibUtilities::SessionReaderSharedPtr pSession)
: Driver(pSession)
{
}
/**
*
*
*/
DriverSteadyState:: ~DriverSteadyState()
{
}
/**
*
*
*/
void DriverSteadyState::v_InitObject(ostream &out)
{
Driver::v_InitObject(out);
}
void DriverSteadyState::v_Execute(ostream &out)
{
//With a loop over "DoSolve", this Driver implements the "encaplulated" Selective Frequency Damping method
//to find the steady state of a flow above the critical Reynolds number.
m_equ[0]->PrintSummary(out);
m_equ[0]->DoInitialise();
......@@ -83,56 +83,56 @@ namespace Nektar
// - SFD Routine -
// Compressible case
NumVar_SFD = m_equ[0]->UpdateFields()[0]->GetCoordim(0);
if (m_session->GetSolverInfo("EqType") == "EulerCFE" ||
NumVar_SFD = m_equ[0]->UpdateFields()[0]->GetCoordim(0);
if (m_session->GetSolverInfo("EqType") == "EulerCFE" ||
m_session->GetSolverInfo("EqType") == "NavierStokesCFE")
{
NumVar_SFD += 2; //Number of variables for the compressible equations
}
Array<OneD, Array<OneD, NekDouble> > q0(NumVar_SFD);
Array<OneD, Array<OneD, NekDouble> > q1(NumVar_SFD);
Array<OneD, Array<OneD, NekDouble> > qBar0(NumVar_SFD);
Array<OneD, Array<OneD, NekDouble> > qBar1(NumVar_SFD);
NekDouble TOL(0);
m_n=0;
m_Check=0;
m_Check=0;
m_session->LoadParameter("FilterWidth", m_Delta0, 1);
m_session->LoadParameter("ControlCoeff",m_X0, 1);
m_session->LoadParameter("TOL", TOL, 1.0e-08);
m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1000);
m_session->LoadParameter("IO_CheckSteps", m_checksteps, 100000);
m_dt = m_equ[0]->GetTimeStep();
//The time-step is applied here to clarify the matrix F expression
m_X = m_X0*m_dt;
m_Delta = m_Delta0/m_dt;
//Exact solution of the Filters equation (ici on a X_tilde et Delta_tile!!!)
c1 = 1.0/(1.0 + m_X*m_Delta);
F11 = c1*(1.0 + m_X*m_Delta*exp(-(m_X + 1.0/m_Delta)));
F12 = c1*(m_X*m_Delta*(1.0 - exp(-(m_X + 1.0/m_Delta))));
F21 = c1*(1.0 - exp(-(m_X + 1.0/m_Delta)));
F22 = c1*(m_X*m_Delta + exp(-(m_X + 1.0/m_Delta)));
cout << "------------------ SFD Parameters ------------------" << endl;
cout << "\tX = " << m_X0 << endl;
cout << "\tDelta = " << m_Delta0 << endl;
cout << "----------------------------------------------------" << endl;
m_equ[0]->SetStepsToOne(); //m_steps is set to 1. Then "m_equ[0]->DoSolve()" will run for only one time step
m_equ[0]->SetStepsToOne(); //m_steps is set to 1. Then "m_equ[0]->DoSolve()" will run for only one time step
ofstream m_file("ConvergenceHistory.txt", ios::out | ios::trunc);
for(int i = 0; i < NumVar_SFD; ++i)
{
q0[i] = Array<OneD, NekDouble> (m_equ[0]->GetTotPoints(), 123.456); //q0 is initialised
q0[i] = Array<OneD, NekDouble> (m_equ[0]->GetTotPoints(), 0.0); //q0 is initialised
qBar0[i] = Array<OneD, NekDouble> (m_equ[0]->GetTotPoints(), 0.0);
m_equ[0]->CopyFromPhysField(i, qBar0[i]); //qBar0 is initially set at beiing equal to the initial conditions provided in the input file
}
MaxNormDiff_q_qBar = 1.0;
MaxNormDiff_q1_q0 = 1.0;
Min_MaxNormDiff_q_qBar = MaxNormDiff_q_qBar;
......@@ -143,42 +143,42 @@ namespace Nektar
{
//First order Splitting with exact resolution of the filters equation
m_equ[0]->DoSolve();
for(int i = 0; i < NumVar_SFD; ++i)
{
m_equ[0]->CopyFromPhysField(i, q0[i]);
//Apply the linear operator F to the outcome of the solver
ExactFilters(i, q0, qBar0, q1, qBar1);
qBar0[i] = qBar1[i];
m_equ[0]->CopyToPhysField(i, q1[i]);
m_equ[0]->CopyToPhysField(i, q1[i]);
}
if(m_infosteps && !((m_n+1)%m_infosteps))
{
ConvergenceHistory(qBar1, q0, MaxNormDiff_q_qBar, MaxNormDiff_q1_q0);
}
if(m_checksteps && m_n&&(!((m_n+1)%m_checksteps)))
{
m_Check++;
m_equ[0]->Checkpoint_Output(m_Check);
{
m_Check++;
m_equ[0]->Checkpoint_Output(m_Check);
}
m_n++;
}
cout << "\nFINAL Filter Width: Delta = " << m_Delta << "; FINAL Control Coeff: X = " << m_X << "\n" << endl;
m_Check++;
m_equ[0]->Checkpoint_Output(m_Check);
m_file.close();
// - End SFD Routine -
m_equ[0]->Output();
// Evaluate and output computation time and solution accuracy.
// The specific format of the error output is essential for the
// regression tests to work.
......@@ -194,8 +194,8 @@ namespace Nektar
}
}
}
void DriverSteadyState::ExactFilters(const int i,
const Array<OneD, const Array<OneD, NekDouble> > &q0,
const Array<OneD, const Array<OneD, NekDouble> > &qBar0,
......@@ -204,66 +204,66 @@ namespace Nektar
{
q1[i] = Array<OneD, NekDouble> (m_equ[0]->GetTotPoints(),0.0);
qBar1[i] = Array<OneD, NekDouble> (m_equ[0]->GetTotPoints(),0.0);
//Exact solution of the Filters equation
Vmath::Svtvp(q1[i].num_elements(), F11, q0[i], 1, q1[i], 1, q1[i], 1 );
Vmath::Svtvp(q1[i].num_elements(), F12, qBar0[i], 1, q1[i], 1, q1[i], 1 );
Vmath::Svtvp(qBar1[i].num_elements(), F21, q0[i], 1, qBar1[i], 1, qBar1[i], 1 );
Vmath::Svtvp(qBar1[i].num_elements(), F22, qBar0[i], 1, qBar1[i], 1, qBar1[i], 1 );
//Exact solution of the Filters equation
Vmath::Svtvp(q1[i].num_elements(), F11, q0[i], 1, q1[i], 1, q1[i], 1 );
Vmath::Svtvp(q1[i].num_elements(), F12, qBar0[i], 1, q1[i], 1, q1[i], 1 );
Vmath::Svtvp(qBar1[i].num_elements(), F21, q0[i], 1, qBar1[i], 1, qBar1[i], 1 );
Vmath::Svtvp(qBar1[i].num_elements(), F22, qBar0[i], 1, qBar1[i], 1, qBar1[i], 1 );
}
void DriverSteadyState::ConvergenceHistory(const Array<OneD, const Array<OneD, NekDouble> > &qBar1,
const Array<OneD, const Array<OneD, NekDouble> > &q0,
NekDouble &MaxNormDiff_q_qBar,
NekDouble &MaxNormDiff_q1_q0)
{
//This routine evaluates |q-qBar|_L2 and save the value in "ConvergenceHistory.txt"
Array<OneD, NekDouble > NormDiff_q_qBar(NumVar_SFD, 1.0);
Array<OneD, NekDouble > NormDiff_q1_q0(NumVar_SFD, 1.0);
MaxNormDiff_q_qBar=0.0;
MaxNormDiff_q1_q0=0.0;
//Norm Calculation
for(int i = 0; i < NumVar_SFD; ++i)
{
//To check convergence of SFD
//NormDiff_q_qBar[i] = m_equ[0]->L2Error(i, qBar1[i], false);
NormDiff_q_qBar[i] = m_equ[0]->LinfError(i, qBar1[i]);
//To check convergence of Navier-Stokes
NormDiff_q1_q0[i] = m_equ[0]->LinfError(i, q0[i]);
if (MaxNormDiff_q_qBar < NormDiff_q_qBar[i])
{
MaxNormDiff_q_qBar = NormDiff_q_qBar[i];
}
if (MaxNormDiff_q1_q0 < NormDiff_q1_q0[i])
{
MaxNormDiff_q1_q0 = NormDiff_q1_q0[i];
}
}
#if NEKTAR_USE_MPI
#if NEKTAR_USE_MPI
MPI_Comm_rank(MPI_COMM_WORLD,&MPIrank);
if (MPIrank==0)
{
cout << "SFD (MPI) - Step: " << m_n+1 << "; Time: " << m_equ[0]->GetFinalTime() << "; |q-qBar|inf = " << MaxNormDiff_q_qBar << "; |q1-q0|inf = " << MaxNormDiff_q1_q0 << ";\t for X = " << m_X0 <<" and Delta = " << m_Delta0 <<endl;
std::ofstream m_file( "ConvergenceHistory.txt", std::ios_base::app);
std::ofstream m_file( "ConvergenceHistory.txt", std::ios_base::app);
m_file << m_n+1 << "\t" << m_equ[0]->GetFinalTime() << "\t" << MaxNormDiff_q_qBar << "\t" << MaxNormDiff_q1_q0 << endl;
m_file.close();
}
#else
cout << "SFD - Step: " << m_n+1 << "; Time: " << m_equ[0]->GetFinalTime() << "; |q-qBar|inf = " << MaxNormDiff_q_qBar << "; |q1-q0|inf = " << MaxNormDiff_q1_q0 << ";\t for X = " << m_X0 <<" and Delta = " << m_Delta0 <<endl;
std::ofstream m_file( "ConvergenceHistory.txt", std::ios_base::app);
std::ofstream m_file( "ConvergenceHistory.txt", std::ios_base::app);
m_file << m_n+1 << "\t" << m_equ[0]->GetFinalTime() << "\t" << MaxNormDiff_q_qBar << "\t" << MaxNormDiff_q1_q0 << endl;
m_file.close();
#endif
#endif
}
}
}
......
......@@ -297,7 +297,6 @@ void AdjointAdvection::v_Advect(
//Evaluation of the base flow for periodic cases
//(it requires fld files)
if(m_slices>1)
{
if (m_session->GetFunctionType("BaseFlow", 0)
......@@ -317,34 +316,34 @@ void AdjointAdvection::v_Advect(
//Evaluate the Adjoint advection term
switch(ndim)
{
// 1D
case 1:
fields[0]->PhysDeriv(inarray[n],grad0);
fields[0]->PhysDeriv(m_baseflow[0],grad_base_u0);
//Evaluate U du'/dx
Vmath::Vmul(nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
//Evaluate U du'/dx+ u' dU/dx
Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
break;
//2D
case 2:
grad1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
fields[0]->PhysDeriv(inarray[n],grad0,grad1);
//Derivates of the base flow
fields[0]-> PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1);
fields[0]-> PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1);
//Since the components of the velocity are passed one by one, it is necessary to distinguish which
//term is consider
switch (n)
{
//x-equation
// 1D
case 1:
fields[0]->PhysDeriv(inarray[n],grad0);
fields[0]->PhysDeriv(m_baseflow[0],grad_base_u0);
//Evaluate U du'/dx
Vmath::Vmul(nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
//Evaluate U du'/dx+ u' dU/dx
Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
break;
//2D
case 2:
grad1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
fields[0]->PhysDeriv(inarray[n],grad0,grad1);
//Derivates of the base flow
fields[0]-> PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1);
fields[0]-> PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1);
//Since the components of the velocity are passed one by one, it is necessary to distinguish which
//term is consider
switch (n)
{
//x-equation
case 0:
// Evaluate U du'/dx
Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
......@@ -358,7 +357,7 @@ void AdjointAdvection::v_Advect(
Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[1],1,outarray[n],1,outarray[n],1);
break;
//y-equation
//y-equation
case 1:
// Evaluate U dv'/dx
Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
......@@ -372,43 +371,43 @@ void AdjointAdvection::v_Advect(
Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
break;
}
break;
break;
//3D
case 3:
//3D
case 3:
grad1 = Array<OneD, NekDouble> (nPointsTot);
grad2 = Array<OneD, NekDouble> (nPointsTot);
grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_w1 = Array<OneD, NekDouble> (nPointsTot);
grad1 = Array<OneD, NekDouble> (nPointsTot);
grad2 = Array<OneD, NekDouble> (nPointsTot);
grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_w1 = Array<OneD, NekDouble> (nPointsTot);
grad_base_u2 = Array<OneD, NekDouble> (nPointsTot);
grad_base_v2 = Array<OneD, NekDouble> (nPointsTot);
grad_base_w2 = Array<OneD, NekDouble> (nPointsTot);
grad_base_u2 = Array<OneD, NekDouble> (nPointsTot);
grad_base_v2 = Array<OneD, NekDouble> (nPointsTot);
grad_base_w2 = Array<OneD, NekDouble> (nPointsTot);
fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1,grad_base_u2);
fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1,grad_base_v2);
fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0, grad_base_w1, grad_base_w2);
fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1,grad_base_u2);
fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1,grad_base_v2);
fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0, grad_base_w1, grad_base_w2);
//HalfMode has W(x,y,t)=0
if(m_HalfMode || m_SingleMode)
//HalfMode has W(x,y,t)=0
if(m_HalfMode || m_SingleMode)
{
for(int i=0; i<grad_base_u2.num_elements();++i)
{
for(int i=0; i<grad_base_u2.num_elements();++i)
{
grad_base_u2[i]=0;
grad_base_v2[i]=0;
grad_base_w2[i]=0;
grad_base_u2[i]=0;
grad_base_v2[i]=0;
grad_base_w2[i]=0;
}
}
}
fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
switch (n)
{
//x-equation
//x-equation
case 0:
//Evaluate U du'/dx
Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
......@@ -425,7 +424,7 @@ void AdjointAdvection::v_Advect(
//Evaluate -(U du'/dx+ V du'/dy+W du'/dz)+u'dU/dx+ v' dV/dx+ w' dW/dz
Vmath::Vvtvp(nPointsTot,grad_base_w0,1,advVel[2],1,outarray[n],1,outarray[n],1);
break;
//y-equation
//y-equation
case 1:
//Evaluate U dv'/dx
Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
......@@ -443,7 +442,7 @@ void AdjointAdvection::v_Advect(
Vmath::Vvtvp(nPointsTot,grad_base_w1,1,advVel[2],1,outarray[n],1,outarray[n],1);
break;
//z-equation
//z-equation
case 2:
//Evaluate U dw'/dx
Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
......@@ -461,11 +460,9 @@ void AdjointAdvection::v_Advect(
Vmath::Vvtvp(nPointsTot,grad_base_w2,1,advVel[2],1,outarray[n],1,outarray[n],1);
break;
}
break;
default:
ASSERTL0(false,"dimension unknown");
break;
default:
ASSERTL0(false,"dimension unknown");
}
Vmath::Neg(nqtot,outarray[n],1);
......
......@@ -57,6 +57,11 @@ public:
static std::string className;
static std::string className2;
void SetSpecHPDealiasing(bool value)
{
m_specHP_dealiasing = value;
}
protected:
NavierStokesAdvection();
......
......@@ -36,7 +36,6 @@
#ifndef NEKTAR_SOLVERS_NOADVECTION_H
#define NEKTAR_SOLVERS_NOADVECTION_H
///#include <IncNavierStokesSolver/AdvectionTerms/AdvectionTerm.h>
#include <SolverUtils/Advection/Advection.h>
namespace Nektar
......
......@@ -357,8 +357,8 @@ namespace Nektar
Deriv = Array<OneD, NekDouble> (nqtot*VelDim);
}
m_advObject->Advect(m_nConvectiveFields,m_fields,
velocity,inarray,outarray,m_time);
m_advObject->Advect(m_nConvectiveFields, m_fields,
velocity, inarray, outarray, m_time);
}
/**
......
......@@ -119,12 +119,6 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, NekDouble> > &numflux);
// AdvectionTermSharedPtr GetAdvObject(void)
// {
// return m_advObject;
// }
int GetNConvectiveFields(void)
{
return m_nConvectiveFields;
......
......@@ -42,7 +42,6 @@
#include <MultiRegions/ExpList.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/TimeIntegration/TimeIntegrationWrapper.h>
///#include <IncNavierStokesSolver/AdvectionTerms/AdvectionTerm.h>
#include <SolverUtils/AdvectionSystem.h>
#include <IncNavierStokesSolver/EquationSystems/Extrapolate.h>
......
......@@ -42,7 +42,6 @@
#include <MultiRegions/ExpList.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/TimeIntegration/TimeIntegrationWrapper.h>
///#include <IncNavierStokesSolver/AdvectionTerms/AdvectionTerm.h>
#include <SolverUtils/AdvectionSystem.h>
#include <IncNavierStokesSolver/EquationSystems/Extrapolate.h>
......
......@@ -5,6 +5,7 @@
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <IncNavierStokesSolver/EquationSystems/IncNavierStokes.h>
#include <IncNavierStokesSolver/AdvectionTerms/NavierStokesAdvection.h>
using namespace Nektar;
using namespace Nektar::SolverUtils;
......@@ -53,17 +54,26 @@ int main(int argc, char *argv[])
NonLinearDealia