Commit 370a310d authored by Dave Moxey's avatar Dave Moxey
Browse files

Adding initial 3D support to the compressible flow solver.

Fixed some minor bugs with parallelisation in compressible flow solver.
Fixed minor formatting issues DisContField2D.


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4022 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 6f017c63
......@@ -302,102 +302,107 @@ namespace Nektar
*/
void DisContField2D::SetUpDG()
{
if(m_trace == NullExpListSharedPtr) // check for multiple calls
// Check for multiple calls
if (m_trace != NullExpListSharedPtr)
{
ExpList1DSharedPtr trace;
SpatialDomains::MeshGraph2DSharedPtr graph2D =
boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(
m_graph);
// Set up matrix map
m_globalBndMat = MemoryManager<GlobalLinSysMap>::
AllocateSharedPtr();
return;
}
ExpList1DSharedPtr trace;
SpatialDomains::MeshGraph2DSharedPtr graph2D =
boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(
m_graph);
// Set up Trace space
trace = MemoryManager<ExpList1D>::AllocateSharedPtr(m_bndCondExpansions, m_bndConditions, *m_exp, graph2D, m_periodicEdges);
// Set up matrix map
m_globalBndMat = MemoryManager<GlobalLinSysMap>::
AllocateSharedPtr();
m_trace = boost::dynamic_pointer_cast<ExpList>(trace);
m_traceMap = MemoryManager<AssemblyMapDG>::
AllocateSharedPtr(m_session, graph2D, trace, *this,
m_bndCondExpansions, m_bndConditions,
m_periodicEdges);
// Set up trace space
trace = MemoryManager<ExpList1D>::AllocateSharedPtr(
m_bndCondExpansions, m_bndConditions, *m_exp,
graph2D, m_periodicEdges);
m_trace = boost::dynamic_pointer_cast<ExpList>(trace);
m_traceMap = MemoryManager<AssemblyMapDG>::
AllocateSharedPtr(m_session, graph2D, trace, *this,
m_bndCondExpansions, m_bndConditions,
m_periodicEdges);
Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
// Scatter trace segments to 2D elements. For each element, we find
// the trace segment associated to each edge. The element then
// retains a pointer to the trace space segments, to ensure
// uniqueness of normals when retrieving from two adjoining elements
// which do not lie in a plane.
SpatialDomains::Geometry1DSharedPtr ElmtSegGeom;
SpatialDomains::Geometry1DSharedPtr TraceSegGeom;
for (int i = 0; i < m_exp->size(); ++i)
// Scatter trace segments to 2D elements. For each element, we find
// the trace segment associated to each edge. The element then
// retains a pointer to the trace space segments, to ensure
// uniqueness of normals when retrieving from two adjoining elements
// which do not lie in a plane.
SpatialDomains::Geometry1DSharedPtr ElmtSegGeom;
SpatialDomains::Geometry1DSharedPtr TraceSegGeom;
for (int i = 0; i < m_exp->size(); ++i)
{
for (int j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
{
for (int j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
{
LocalRegions::Expansion2DSharedPtr exp2d =
boost::dynamic_pointer_cast<
LocalRegions::Expansion2D>((*m_exp)[i]);
LocalRegions::Expansion1DSharedPtr exp1d =
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(elmtToTrace[i][j]);
LocalRegions::ExpansionSharedPtr exp =
boost::dynamic_pointer_cast<
LocalRegions::Expansion> (elmtToTrace[i][j]);
exp2d->SetEdgeExp (j, exp );
exp1d->SetAdjacentElementExp(j, exp2d);
}
LocalRegions::Expansion2DSharedPtr exp2d =
boost::dynamic_pointer_cast<
LocalRegions::Expansion2D>((*m_exp)[i]);
LocalRegions::Expansion1DSharedPtr exp1d =
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(elmtToTrace[i][j]);
LocalRegions::ExpansionSharedPtr exp =
boost::dynamic_pointer_cast<
LocalRegions::Expansion> (elmtToTrace[i][j]);
exp2d->SetEdgeExp (j, exp );
exp1d->SetAdjacentElementExp(j, exp2d);
}
}
// Set up physical normals
SetUpPhysNormals();
// Set up physical normals
SetUpPhysNormals();
// Set up information for parallel jobs.
for (int i = 0; i < m_trace->GetExpSize(); ++i)
{
LocalRegions::Expansion1DSharedPtr traceEl =
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(m_trace->GetExp(i));
// Set up information for parallel jobs.
for (int i = 0; i < m_trace->GetExpSize(); ++i)
{
LocalRegions::Expansion1DSharedPtr traceEl =
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(m_trace->GetExp(i));
int offset = m_trace->GetPhys_Offset(i);
int offset = m_trace->GetPhys_Offset(i);
if (m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
{
traceEl->GetLeftAdjacentElementExp()->NegateEdgeNormal(
traceEl->GetLeftAdjacentElementEdge());
}
if (m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
{
traceEl->GetLeftAdjacentElementExp()->NegateEdgeNormal(
traceEl->GetLeftAdjacentElementEdge());
}
}
int cnt, n, e;
int cnt, n, e;
// Identify boundary edges
for(cnt = 0, n = 0; n < m_bndCondExpansions.num_elements(); ++n)
// Identify boundary edges
for(cnt = 0, n = 0; n < m_bndCondExpansions.num_elements(); ++n)
{
if (m_bndConditions[n]->GetBoundaryConditionType() !=
SpatialDomains::ePeriodic)
{
if (m_bndConditions[n]->GetBoundaryConditionType() !=
SpatialDomains::ePeriodic)
for(e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
{
for(e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
{
m_boundaryEdges.insert(m_trace->GetOffset_Elmt_Id(
m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e)));
}
m_boundaryEdges.insert(m_trace->GetOffset_Elmt_Id(
m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e)));
}
cnt += m_bndCondExpansions[n]->GetExpSize();
}
cnt += m_bndCondExpansions[n]->GetExpSize();
}
// Set up information for periodic boundary conditions.
for (n = 0; n < m_exp->size(); ++n)
// Set up information for periodic boundary conditions.
for (n = 0; n < m_exp->size(); ++n)
{
for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
{
for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
{
map<int,int>::iterator it = m_periodicEdges.find(
(*m_exp)[n]->GetGeom2D()->GetEid(e));
map<int,int>::iterator it = m_periodicEdges.find(
(*m_exp)[n]->GetGeom2D()->GetEid(e));
if (it != m_periodicEdges.end())
{
m_perEdgeToExpMap[it->first] = make_pair(n, e);
}
if (it != m_periodicEdges.end())
{
m_perEdgeToExpMap[it->first] = make_pair(n, e);
}
}
}
......@@ -1426,10 +1431,10 @@ namespace Nektar
* solving the system:
*
* \f[
* \begin{align*}
* \begin{aligned}
* (\nabla w, \nabla u^*) &= (\nabla w, u), \\
* \langle \nabla u^*, 1 \rangle &= \langle \grad u, 1 \rangle
* \end{align*}
* \langle \nabla u^*, 1 \rangle &= \langle \nabla u, 1 \rangle
* \end{aligned}
* \f]
*
* where \f$ u \f$ corresponds with the current solution as stored
......
......@@ -921,7 +921,7 @@ namespace Nektar
<< ntot << "\" NumberOfCells=\""
<< ntotminus << "\">" << endl;
outfile << " <Points>" << endl;
outfile << " <DataArray type=\"Float32\" "
outfile << " <DataArray type=\"Float64\" "
<< "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
outfile << " ";
for (i = 0; i < ntot; ++i)
......
......@@ -58,7 +58,7 @@ namespace Nektar
int nPointsTot = fields[0]->GetTotPoints();
int nCoeffs = fields[0]->GetNcoeffs();
int nTracePointsTot = fields[0]->GetTrace()->GetTotPoints();
Array<OneD, Array<OneD, NekDouble> > fluxvector(nVelDim);
Array<OneD, Array<OneD, NekDouble> > tmp (nConvectiveFields);
......
......@@ -58,7 +58,7 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
virtual void v_Advect(
const int nConvectiveFields,
const int nConvective,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &advVel,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
......
......@@ -574,7 +574,14 @@ namespace Nektar
}
break;
case 3:
ASSERTL0(false, "Force function not implemented for 3D.");
for (int i = 0; i < m_forces.num_elements(); i++)
{
m_forces[i] = MemoryManager<MultiRegions::ExpList3D>
::AllocateSharedPtr(*boost::static_pointer_cast<
MultiRegions::ExpList3D>(m_fields[i]));
Vmath::Zero(nq, m_forces[i]->UpdatePhys(), 1);
}
break;
}
// Check for file
......
......@@ -36,6 +36,22 @@
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <SolverUtils/RiemannSolvers/RiemannSolver.h>
#define EPSILON 0.000001
#define CROSS(dest, v1, v2){ \
dest[0] = v1[1] * v2[2] - v1[2] * v2[1]; \
dest[1] = v1[2] * v2[0] - v1[0] * v2[2]; \
dest[2] = v1[0] * v2[1] - v1[1] * v2[0];}
#define DOT(v1, v2) (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2])
#define SUB(dest, v1, v2){ \
dest[0] = v1[0] - v2[0]; \
dest[1] = v1[1] - v2[1]; \
dest[2] = v1[2] - v2[2];}
using namespace std;
namespace Nektar
{
namespace SolverUtils
......@@ -131,7 +147,7 @@ namespace Nektar
switch(normals.num_elements())
{
case 1:
ASSERTL0(false, "1D and 3D not implemented yet.");
ASSERTL0(false, "1D not implemented yet.");
break;
case 2:
......@@ -142,8 +158,10 @@ namespace Nektar
for (i = 0; i < normals[0].num_elements(); ++i)
{
double tmp1 = inarray[vx][i]*normals[0][i] + inarray[vy][i]*normals[1][i];
double tmp2 = -inarray[vx][i]*normals[1][i] + inarray[vy][i]*normals[0][i];
double tmp1 = inarray[vx][i]*normals[0][i] +
inarray[vy][i]*normals[1][i];
double tmp2 = -inarray[vx][i]*normals[1][i] +
inarray[vy][i]*normals[0][i];
outarray[vx][i] = tmp1;
outarray[vy][i] = tmp2;
}
......@@ -155,15 +173,52 @@ namespace Nektar
continue;
}
Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1, outarray[i], 1);
Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
outarray[i], 1);
}
break;
}
case 3:
ASSERTL0(false, "1D and 3D not implemented yet.");
break;
{
int vx = (int)velLoc[0];
int vy = (int)velLoc[1];
int vz = (int)velLoc[2];
// Generate matrices if they don't already exist.
if (m_rotMatrices.size() == 0)
{
GenerateRotationMatrices();
}
// Apply rotation matrices.
for (int i = 0; i < normals[0].num_elements(); ++i)
{
DNekMatSharedPtr m = m_rotMatrices[i];
double tmp1 = inarray[vx][i]*(*m)(0,0) +
inarray[vy][i]*(*m)(0,1) + inarray[vz][i]*(*m)(0,2);
double tmp2 = inarray[vx][i]*(*m)(1,0) +
inarray[vy][i]*(*m)(1,1) + inarray[vz][i]*(*m)(1,2);
double tmp3 = inarray[vx][i]*(*m)(2,0) +
inarray[vy][i]*(*m)(2,1) + inarray[vz][i]*(*m)(2,2);
outarray[vx][i] = tmp1;
outarray[vy][i] = tmp2;
outarray[vz][i] = tmp3;
}
for (int i = 0; i < inarray.num_elements(); ++i)
{
if (i == vx || i == vy || i == vz)
{
continue;
}
Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
outarray[i], 1);
}
break;
}
default:
ASSERTL0(false, "Invalid space dimension.");
break;
......@@ -182,13 +237,14 @@ namespace Nektar
const Array<OneD, const Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
const Array<OneD, const Array<OneD, NekDouble> > &normals = m_vectors["N"]();
const Array<OneD, const Array<OneD, NekDouble> > &normals =
m_vectors["N"]();
const Array<OneD, NekDouble> &velLoc = m_scalars["velLoc"]();
switch(normals.num_elements())
{
case 1:
ASSERTL0(false, "1D and 3D not implemented yet.");
ASSERTL0(false, "1D not implemented yet.");
break;
case 2:
......@@ -199,8 +255,10 @@ namespace Nektar
for (i = 0; i < normals[0].num_elements(); ++i)
{
double tmp1 = inarray[vx][i]*normals[0][i] - inarray[vy][i]*normals[1][i];
double tmp2 = inarray[vx][i]*normals[1][i] + inarray[vy][i]*normals[0][i];
double tmp1 = inarray[vx][i]*normals[0][i] -
inarray[vy][i]*normals[1][i];
double tmp2 = inarray[vx][i]*normals[1][i] +
inarray[vy][i]*normals[0][i];
outarray[vx][i] = tmp1;
outarray[vy][i] = tmp2;
}
......@@ -212,15 +270,52 @@ namespace Nektar
continue;
}
Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1, outarray[i], 1);
Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
outarray[i], 1);
}
break;
}
case 3:
ASSERTL0(false, "1D and 3D not implemented yet.");
break;
{
int vx = (int)velLoc[0];
int vy = (int)velLoc[1];
int vz = (int)velLoc[2];
// Generate matrices if they don't already exist.
if (m_rotMatrices.size() == 0)
{
GenerateRotationMatrices();
}
// Apply rotation matrices.
for (int i = 0; i < normals[0].num_elements(); ++i)
{
DNekMatSharedPtr m = m_invRotMatrices[i];
double tmp1 = inarray[vx][i]*(*m)(0,0) +
inarray[vy][i]*(*m)(0,1) + inarray[vz][i]*(*m)(0,2);
double tmp2 = inarray[vx][i]*(*m)(1,0) +
inarray[vy][i]*(*m)(1,1) + inarray[vz][i]*(*m)(1,2);
double tmp3 = inarray[vx][i]*(*m)(2,0) +
inarray[vy][i]*(*m)(2,1) + inarray[vz][i]*(*m)(2,2);
outarray[vx][i] = tmp1;
outarray[vy][i] = tmp2;
outarray[vz][i] = tmp3;
}
for (int i = 0; i < inarray.num_elements(); ++i)
{
if (i == vx || i == vy || i == vz)
{
continue;
}
Vmath::Vcopy(inarray[i].num_elements(), inarray[i], 1,
outarray[i], 1);
}
break;
}
default:
ASSERTL0(false, "Invalid space dimension.");
break;
......@@ -265,5 +360,139 @@ namespace Nektar
return it != m_params.end();
}
/**
* @brief Generate rotation matrices for 3D expansions.
*/
void RiemannSolver::GenerateRotationMatrices()
{
const Array<OneD, const Array<OneD, NekDouble> > &normals =
m_vectors["N"]();
Array<OneD, NekDouble> xdir(3,0.0);
Array<OneD, NekDouble> tn (3);
xdir[0] = 1.0;
for (int i = 0; i < normals[0].num_elements(); ++i)
{
DNekMatSharedPtr mat = MemoryManager<DNekMat>::
AllocateSharedPtr(3,3);
// Generate matrix which takes us from (1,0,0) vector to trace
// normal.
tn[0] = normals[0][i];
tn[1] = normals[1][i];
tn[2] = normals[2][i];
FromToRotation(tn, xdir, mat);
// Generate inverse rotation matrix.
DNekMatSharedPtr mat2 = MemoryManager<DNekMat>::
AllocateSharedPtr(*mat);
mat2->Invert();
m_rotMatrices. push_back(mat);
m_invRotMatrices.push_back(mat2);
}
}
/**
* @brief A function for creating a rotation matrix that rotates a
* vector @a from into another vector @a to.
*
* Authors: Tomas Möller, John Hughes
* "Efficiently Building a Matrix to Rotate One Vector to
* Another" Journal of Graphics Tools, 4(4):1-4, 1999
*
* @param from Normalised 3-vector to rotate from.
* @param to Normalised 3-vector to rotate to.
* @param out Resulting 3x3 rotation matrix (row-major order).
*/
void RiemannSolver::FromToRotation(
Array<OneD, const NekDouble> &from,
Array<OneD, const NekDouble> &to,
DNekMatSharedPtr mat)
{
NekDouble v[3];
NekDouble e, h, f;
CROSS(v, from, to);
e = DOT(from, to);
f = (e < 0)? -e:e;
if (f > 1.0 - EPSILON)
{
NekDouble u[3], v[3];
NekDouble x[3];
NekDouble c1, c2, c3;
int i, j;
x[0] = (from[0] > 0.0)? from[0] : -from[0];
x[1] = (from[1] > 0.0)? from[1] : -from[1];
x[2] = (from[2] > 0.0)? from[2] : -from[2];
if (x[0] < x[1])
{
if (x[0] < x[2])
{
x[0] = 1.0; x[1] = x[2] = 0.0;
}
else
{
x[2] = 1.0; x[0] = x[1] = 0.0;
}
}
else
{
if (x[1] < x[2])
{
x[1] = 1.0; x[0] = x[2] = 0.0;
}
else
{
x[2] = 1.0; x[0] = x[1] = 0.0;
}
}
u[0] = x[0] - from[0];
u[1] = x[1] - from[1];
u[2] = x[2] - from[2];
v[0] = x[0] - to [0];
v[1] = x[1] - to [1];
v[2] = x[2] - to [2];
c1 = 2.0 / DOT(u, u);
c2 = 2.0 / DOT(v, v);
c3 = c1 * c2 * DOT(u, v);
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
(*mat)(i,j) = - c1 * u[i] * u[j]
- c2 * v[i] * v[j]
+ c3 * v[i] * u[j];
}
(*mat)(i,i) += 1.0;
}
}
else
{
NekDouble hvx, hvz, hvxy, hvxz, hvyz;
h = 1.0/(1.0 + e);
hvx = h * v[0];
hvz = h * v[2];
hvxy = hvx * v[1];
hvxz = hvx * v[2];
hvyz = hvz * v[1];
(*mat)(0,0) = e + hvx * v[0];
(*mat)(0,1) = hvxy - v[2];
(*mat)(0,2) = hvxz + v[1];
(*mat)(1,0) = hvxy + v[2];
(*mat)(1,1) = e + h * v[1] * v[1];
(*mat)(1,2) = hvyz - v[0];
(*mat)(2,0) = hvxz - v[1];
(*mat)(2,1) = hvyz + v[0];
(*mat)(2,2) = e + hvz * v[2];
}
}
}
}
......@@ -40,6 +40,7 @@
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
#include <SolverUtils/SolverUtilsDeclspec.h>
#include <boost/function.hpp>
......@@ -99,6 +100,10 @@ namespace Nektar
std::map<std::string, RSVecFuncType> m_vectors;