Commit 764acd12 authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'feature/diffusion' of localhost:nektar

parents ae7f46b9 2253fab3
......@@ -33,7 +33,7 @@
//
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/Advection.h>
#include <SolverUtils/Advection/Advection.h>
namespace Nektar
{
......@@ -53,28 +53,7 @@ namespace Nektar
{
v_InitObject(pSession, pFields);
}
void Advection::SetupMetrics(
const LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
v_SetupMetrics(pSession, pFields);
}
void Advection::SetupCFunctions(
const LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
v_SetupCFunctions(pSession, pFields);
}
void Advection::SetupInterpolationMatrices(
const LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
v_SetupInterpolationMatrices(pSession, pFields);
}
void Advection::Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
......@@ -84,38 +63,5 @@ namespace Nektar
{
v_Advect(nConvectiveFields, fields, advVel, inarray, outarray);
}
void Advection::DivCFlux_1D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, const NekDouble> &fluxX1,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux)
{
v_DivCFlux_1D(nConvectiveFields, fields, fluxX1, numericalFlux, divCFlux);
}
void Advection::DivCFlux_2D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, const NekDouble> &fluxX1,
const Array<OneD, const NekDouble> &fluxX2,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux)
{
v_DivCFlux_2D(nConvectiveFields, fields, fluxX1, fluxX2, numericalFlux, divCFlux);
}
void Advection::DivCFlux_3D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, const NekDouble> &fluxX1,
const Array<OneD, const NekDouble> &fluxX2,
const Array<OneD, const NekDouble> &fluxX3,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux)
{
v_DivCFlux_3D(nConvectiveFields, fields, fluxX1, fluxX2, fluxX3, numericalFlux, divCFlux);
}
}
}
......@@ -60,50 +60,14 @@ namespace Nektar
SOLVER_UTILS_EXPORT void InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
SOLVER_UTILS_EXPORT void SetupMetrics(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
SOLVER_UTILS_EXPORT void SetupCFunctions(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
SOLVER_UTILS_EXPORT void SetupInterpolationMatrices(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields);
SOLVER_UTILS_EXPORT void Advect(
const int nConvectiveFields,
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);
SOLVER_UTILS_EXPORT void DivCFlux_1D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, const NekDouble> &fluxX1,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux);
SOLVER_UTILS_EXPORT void DivCFlux_2D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, const NekDouble> &fluxX1,
const Array<OneD, const NekDouble> &fluxX2,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux);
SOLVER_UTILS_EXPORT void DivCFlux_3D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, const NekDouble> &fluxX1,
const Array<OneD, const NekDouble> &fluxX2,
const Array<OneD, const NekDouble> &fluxX3,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux);
template<typename FuncPointerT, typename ObjectPointerT>
void SetFluxVector(FuncPointerT func, ObjectPointerT obj)
{
......@@ -122,28 +86,7 @@ namespace Nektar
{
};
virtual void v_SetupMetrics(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
};
virtual void v_SetupCFunctions(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
};
virtual void v_SetupInterpolationMatrices(
LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
};
virtual void v_Advect(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
......@@ -151,39 +94,6 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)=0;
virtual void v_DivCFlux_1D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, const NekDouble> &fluxX1,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux)
{
};
virtual void v_DivCFlux_2D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, const NekDouble> &fluxX1,
const Array<OneD, const NekDouble> &fluxX2,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux)
{
};
virtual void v_DivCFlux_3D(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, const NekDouble> &fluxX1,
const Array<OneD, const NekDouble> &fluxX2,
const Array<OneD, const NekDouble> &fluxX3,
const Array<OneD, const NekDouble> &numericalFlux,
Array<OneD, NekDouble> &divCFlux)
{
};
AdvectionFluxVecCB m_fluxVector;
RiemannSolverSharedPtr m_riemann;
};
......
......@@ -33,7 +33,7 @@
//
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/AdvectionFR.h>
#include <SolverUtils/Advection/AdvectionFR.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/Foundations/Basis.h>
#include <LibUtilities/Foundations/Points.h>
......@@ -895,9 +895,6 @@ namespace Nektar
// 2D-Problems
case 2:
{
// =========================================================
// IMPLEMENTATION FOR SYSTEMS OF EQUATIONS
// =========================================================
Array<OneD, Array<OneD, NekDouble> > fluxvector(
nDimensions);
Array<OneD, Array<OneD, NekDouble> > fluxvectorX1(
......@@ -947,10 +944,7 @@ namespace Nektar
for (n = 0; n < nElements; n++)
{
// -------------------------------------------------
// Proper implementation discontinuous flux:
// first attempt
// -------------------------------------------------
// Discontinuous flux
nLocalSolutionPts = fields[0]->GetExp(n)->
GetTotPoints();
phys_offset = fields[0]->GetPhys_Offset(n);
......@@ -1007,23 +1001,6 @@ namespace Nektar
fields[0]->GetExp(n)->StdPhysDeriv(1,
auxArray1 = g_hat,
auxArray2 = DfluxvectorX2[i] + phys_offset);
// -------------------------------------------------
/*
// -------------------------------------------------
// Slightly incorrect implementation
// -------------------------------------------------
phys_offset = fields[0]->GetPhys_Offset(n);
fields[i]->GetExp(n)->PhysDeriv(
0, auxArray1 = fluxvectorX1[i] + phys_offset,
auxArray2 = DfluxvectorX1[i] + phys_offset);
fields[i]->GetExp(n)->PhysDeriv(
1, auxArray1 = fluxvectorX2[i] + phys_offset,
auxArray2 = DfluxvectorX2[i] + phys_offset);
// -------------------------------------------------
*/
}
// Divergence of the discontinuous flux
......@@ -1040,20 +1017,13 @@ namespace Nektar
numflux[i],
divFC[i]);
// -----------------------------------------------------
// Divergence of the final flux for
// correct implementation
// -----------------------------------------------------
// Divergence of the final flux
Vmath::Vadd(nSolutionPts,
auxArray1 = divFD[i], 1,
auxArray2 = divFC[i], 1,
auxArray3 = outarray[i], 1);
// -----------------------------------------------------
// -----------------------------------------------------
// Multiply by the metric terms for
// correct implementation
// -----------------------------------------------------
// Multiply by the metric terms
for (n = 0; n < nElements; ++n)
{
nLocalSolutionPts = fields[0]->
......@@ -1081,195 +1051,9 @@ namespace Nektar
auxArray2 = outarray[i] + phys_offset, 1);
}
}
// -----------------------------------------------------
/*
// -----------------------------------------------------
// Multiply by the metric terms for sligthly
// incorrect implementation
// -----------------------------------------------------
for (n = 0; n < nElements; ++n)
{
nLocalSolutionPts = fields[0]->
GetExp(n)->GetTotPoints();
phys_offset = fields[0]->GetPhys_Offset(n);
jac = fields[0]->GetExp(n)->GetGeom2D()->GetJac();
gmat = fields[0]->GetExp(n)->GetGeom2D()->GetGmat();
if (fields[0]->GetExp(n)->GetGeom2D()->GetGtype()
== SpatialDomains::eDeformed)
{
for (j = 0; j < nLocalSolutionPts; j++)
{
divFC[i][phys_offset + j] =
divFC[i][phys_offset + j]/jac[j];
}
}
else
{
Vmath::Smul(
nLocalSolutionPts,
1/jac[0],
auxArray1 = divFC[i] + phys_offset, 1,
auxArray2 = divFC[i] + phys_offset, 1);
}
}
// -----------------------------------------------------
// -----------------------------------------------------
// Divergence of the final flux for sligthly
// -----------------------------------------------------
Vmath::Vadd(nSolutionPts,
auxArray1 = divFD[i], 1,
auxArray2 = divFC[i], 1,
auxArray3 = outarray[i], 1);
// -----------------------------------------------------
*/
} // close nConvectiveFields loop
// =========================================================
// =========================================================
// =========================================================
/*
// =========================================================
// IMPLEMENTATION FOR SINGLE EQUATION
// =========================================================
Array<OneD, Array<OneD, NekDouble> > fluxvector(
nDimensions);
Array<OneD, Array<OneD, NekDouble> > Dfluxvector(
nDimensions);
for(i = 0; i < nDimensions; ++i)
{
fluxvector[i] = Array<OneD, NekDouble>(nSolutionPts);
Dfluxvector[i] = Array<OneD, NekDouble>(nSolutionPts);
}
// Get the discontinuous flux FD ("i" is used by inarray)
for(i = 0; i < nConvectiveFields; ++i)
{
// Get the ith component of the flux vector
m_fluxVector(i, inarray, fluxvector);
}
// Divergence of the discontinuous flux
for (n = 0; n < nElements; ++n)
{
nLocalSolutionPts = fields[0]->
GetExp(n)->GetTotPoints();
phys_offset = fields[0]->GetPhys_Offset(n);
jac = fields[0]->GetExp(n)->GetGeom2D()->GetJac();
gmat = fields[0]->GetExp(n)->GetGeom2D()->GetGmat();
Array<OneD, NekDouble> f_hat(nLocalSolutionPts, 0.0);
Array<OneD, NekDouble> g_hat(nLocalSolutionPts, 0.0);
if (fields[0]->GetExp(n)->GetGeom2D()->GetGtype()
== SpatialDomains::eDeformed)
{
for (j = 0; j < nLocalSolutionPts; j++)
{
f_hat[j] =
(fluxvector[0][j+phys_offset]
* gmat[0][j] +
fluxvector[1][j+phys_offset]
* gmat[2][j]) * jac[j];
g_hat[j] =
(fluxvector[0][j+phys_offset]
* gmat[1][j] +
fluxvector[1][j+phys_offset]
* gmat[3][j]) * jac[j];
}
}
else
{
for (j = 0; j < nLocalSolutionPts; j++)
{
f_hat[j] =
(fluxvector[0][j+phys_offset]
* gmat[0][0] +
fluxvector[1][j+phys_offset]
* gmat[2][0]) * jac[0];
g_hat[j] =
(fluxvector[0][j+phys_offset]
* gmat[1][0] +
fluxvector[1][j+phys_offset]
* gmat[3][0])*jac[0];
}
}
fields[0]->GetExp(n)->StdPhysDeriv(0,
auxArray1 = f_hat,
auxArray2 = Dfluxvector[0] + phys_offset);
fields[0]->GetExp(n)->StdPhysDeriv(1,
auxArray1 = g_hat,
auxArray2 = Dfluxvector[1] + phys_offset);
}
// Computation of the divergence of the discontinuous flux
Array<OneD, NekDouble> divFD(nSolutionPts, 0.0);
Vmath::Vadd(nSolutionPts,
auxArray1 = Dfluxvector[0], 1,
auxArray2 = Dfluxvector[1], 1,
auxArray3 = divFD, 1);
// Computation of the divergence of the correction flux
Array<OneD, NekDouble> divFC(nSolutionPts, 0.0);
for (j = 0; j < nConvectiveFields; ++j)
{
v_DivCFlux_2D(nConvectiveFields,
fields,
fluxvector[0],
fluxvector[1],
numflux[j],
divFC);
}
// Computation of the divergence of the final flux
Vmath::Vadd(nSolutionPts,
auxArray1 = divFD, 1,
auxArray2 = divFC, 1,
auxArray3 = outarray[0], 1);
// Multiplication by the inverse of the jacobian
for (n = 0; n < nElements; ++n)
{
nLocalSolutionPts = fields[0]->GetExp(n)->GetTotPoints();
phys_offset = fields[0]->GetPhys_Offset(n);
jac = fields[0]->GetExp(n)->GetGeom2D()->GetJac();
gmat = fields[0]->GetExp(n)->GetGeom2D()->GetGmat();
if (fields[0]->GetExp(n)->GetGeom2D()->GetGtype()
== SpatialDomains::eDeformed)
{
for (j = 0; j < nLocalSolutionPts; j++)
{
outarray[0][phys_offset + j] =
outarray[0][phys_offset + j]/jac[j];
}
}
else
{
Vmath::Smul(
nLocalSolutionPts,
1/jac[0],
auxArray1 = outarray[0] + phys_offset, 1,
auxArray2 = outarray[0] + phys_offset, 1);
}
}
// =========================================================
// =========================================================
// =========================================================
*/
} // close nConvectiveFields loop
break;
}
......
......@@ -36,7 +36,7 @@
#ifndef NEKTAR_SOLVERUTILS_ADVECTIONFR
#define NEKTAR_SOLVERUTILS_ADVECTIONFR
#include <SolverUtils/Advection.h>
#include <SolverUtils/Advection/Advection.h>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
......
......@@ -33,7 +33,7 @@
//
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/AdvectionNonConservative.h>
#include <SolverUtils/Advection/AdvectionNonConservative.h>
namespace Nektar
{
......
......@@ -36,7 +36,7 @@
#ifndef NEKTAR_SOLVERUTILS_ADVECTIONNONCONSERVATIVE
#define NEKTAR_SOLVERUTILS_ADVECTIONNONCONSERVATIVE
#include <SolverUtils/Advection.h>
#include <SolverUtils/Advection/Advection.h>
namespace Nektar
{
......
......@@ -33,7 +33,7 @@
//
///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/AdvectionWeakDG.h>
#include <SolverUtils/Advection/AdvectionWeakDG.h>
namespace Nektar
{
......@@ -75,7 +75,7 @@ namespace Nektar
{
tmp[i] = Array<OneD,NekDouble>(nCoeffs, 0.0);
// Get the ith component of the flux vector in (physical space)
// Get the ith component of the flux vector in physical space
m_fluxVector(i, inarray, fluxvector);
for (j = 0; j < nVelDim; ++j)
......
......@@ -29,14 +29,14 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Upwind Riemann solver.
// Description: Weak DG advection class.
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERUTILS_ADVECTIONWEAKDG
#define NEKTAR_SOLVERUTILS_ADVECTIONWEAKDG
#include <SolverUtils/Advection.h>
#include <SolverUtils/Advection/Advection.h>
namespace Nektar
{
......
SET(SOLVER_UTILS_SOURCES
Advection.cpp
AdvectionFR.cpp
AdvectionNonConservative.cpp
AdvectionWeakDG.cpp
Advection/Advection.cpp
Advection/AdvectionFR.cpp
Advection/AdvectionNonConservative.cpp
Advection/AdvectionWeakDG.cpp
Diffusion/Diffusion.cpp
Diffusion/DiffusionLDG.cpp
Diffusion/DiffusionLDGNS.cpp
Driver.cpp
DriverArnoldi.cpp
DriverModifiedArnoldi.cpp
......@@ -19,10 +22,13 @@ SET(SOLVER_UTILS_SOURCES
)
SET(SOLVER_UTILS_HEADERS
Advection.h
AdvectionFR.h
AdvectionNonConservative.h
AdvectionWeakDG.h
Advection/Advection.h
Advection/AdvectionFR.h
Advection/AdvectionNonConservative.h
Advection/AdvectionWeakDG.h
Diffusion/Diffusion.h
Diffusion/DiffusionLDG.h
Diffusion/DiffusionLDGNS.h
Driver.h
DriverArnoldi.h
DriverModifiedArnoldi.h
......
///////////////////////////////////////////////////////////////////////////////
//
// File: Diffusion.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: Abstract base class for diffusion.
//
///////////////////////////////////////////////////////////////////////////////