Commit c3923ac4 authored by Spencer Sherwin's avatar Spencer Sherwin Committed by Chris Cantwell
Browse files

Update to put in radiation outflow condition and also some playing around with matrix multiplies

Conflicts:
	solvers/IncNavierStokesSolver/EquationSystems/IncNavierStokes.cpp
parent 7955376b
......@@ -106,8 +106,8 @@ namespace Nektar
CanGetRawPtr<NekMatrix<LhsDataType, MatrixType> >
>::type* p=0)
{
int m = lhs.GetRows();
int n = lhs.GetColumns();
int m = lhs.GetRows();
int n = lhs.GetColumns();
int kl = lhs.GetNumberOfSubDiagonals();
int ku = lhs.GetNumberOfSuperDiagonals();
double alpha = lhs.Scale();
......
......@@ -137,6 +137,7 @@ namespace Nektar
NekMatrix<NekMatrix< DataType, InnerMatrixType>, ScaledMatrixTag>::CreateWrapper(const boost::shared_ptr<NekMatrix<NekMatrix< DataType, InnerMatrixType>, ScaledMatrixTag> >& rhs)
{
return boost::shared_ptr<ThisType>(new ThisType(*rhs));
//return boost::shared_ptr<ThisType>(new NekMatrix<NekMatrix< DataType, InnerMatrixType>, ScaledMatrixTag>(*rhs));
}
......
......@@ -469,7 +469,7 @@ namespace Nektar
{
return NekMatrix<DataType, StandardMatrixTag>(rhs, eWrapper);
}
template<typename DataType>
boost::shared_ptr<NekMatrix<DataType, StandardMatrixTag> > NekMatrix<DataType, StandardMatrixTag>::CreateWrapper(const boost::shared_ptr<NekMatrix<DataType, StandardMatrixTag> >& rhs)
{
......
......@@ -1036,7 +1036,8 @@ namespace Nektar
Array<OneD,int> sign;
StdRegions::VarCoeffMap varcoeffs;
varcoeffs[StdRegions::eVarCoeffPrimative] = primCoeffs;
//varcoeffs[StdRegions::eVarCoeffPrimative] = primCoeffs;
varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
LocalRegions::MatrixKey mkey(StdRegions::eMass,StdRegions::eSegment, *edgeExp, StdRegions::NullConstFactorMap, varcoeffs);
DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
......@@ -1158,7 +1159,8 @@ namespace Nektar
Array<OneD,int> sign;
StdRegions::VarCoeffMap varcoeffs;
varcoeffs[StdRegions::eVarCoeffPrimative] = primCoeffs;
//varcoeffs[StdRegions::eVarCoeffPrimative] = primCoeffs;
varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
LocalRegions::MatrixKey mkey(StdRegions::eMass,StdRegions::eSegment, *edgeExp, StdRegions::NullConstFactorMap, varcoeffs);
DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
......
......@@ -978,8 +978,9 @@ namespace Nektar
Array<OneD, int> sign;
StdRegions::VarCoeffMap varcoeffs;
varcoeffs[StdRegions::eVarCoeffPrimative] = primCoeffs;
//varcoeffs[StdRegions::eVarCoeffPrimative] = primCoeffs;
varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
StdRegions::ExpansionType expType =
faceExp->DetExpansionType();
......
......@@ -1297,7 +1297,8 @@ cout<<"deps/dx ="<<inarray_d0[i]<<" deps/dy="<<inarray_d1[i]<<endl;
{
case StdRegions::eMass:
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)||
(mkey.GetNVarCoeff()))
{
fac = 1.0;
goto UseLocRegionsMatrix;
......@@ -1311,7 +1312,8 @@ cout<<"deps/dx ="<<inarray_d0[i]<<" deps/dy="<<inarray_d1[i]<<endl;
break;
case StdRegions::eInvMass:
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)||
(mkey.GetNVarCoeff()))
{
NekDouble one = 1.0;
StdRegions::StdMatrixKey masskey(StdRegions::eMass,DetExpansionType(),
......
......@@ -1680,6 +1680,8 @@ namespace Nektar
int nblks = 2;
returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks, nblks, exp_size, exp_size); //Really need a constructor which takes Arrays
NekDouble factor = 1.0;
MatrixStorage AMatStorage = eFULL;
//MatrixStorage AMatStorage = eUPPER_TRIANGULAR;
switch(mkey.GetMatrixType())
{
......@@ -1726,7 +1728,7 @@ namespace Nektar
NekDouble invfactor = 1.0/factor;
NekDouble one = 1.0;
DNekScalMat &mat = *GetLocMatrix(mkey);
DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry,AMatStorage);
DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
......
......@@ -1097,7 +1097,7 @@ namespace Nektar
}
/**
* @brief Set up a list of elemeent IDs and edge IDs that link to the
* @brief Set up a list of element IDs and edge IDs that link to the
* boundary conditions.
*/
void DisContField2D::v_GetBoundaryToElmtMap(
......@@ -1690,7 +1690,7 @@ namespace Nektar
std::vector<SpatialDomains::
FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
std::vector<std::vector<NekDouble> > FieldData;
m_graph->Import(filebcs,FieldDef, FieldData);
......
......@@ -297,14 +297,14 @@ namespace Nektar
}
}
/**
/**
* Fills the list of local expansions with the segments in one
* subdomain specified in an inputfile by \a domain. This
* CompositeMap contains a list of Composites which define the
* subdomains.
* subdomain specified in an inputfile by \a domain. This
* CompositeMap contains a list of Composites which define the
* subdomains.
* @param domain A domain, comprising of one or more composite
* regions.
* @param i Index of currently processed subdomain
* @param i Index of currently processed subdomain
* @param graph1D A mesh, containing information about the
* domain and the spectral/hp element expansion.
* @param DeclareCoeffPhysArrays If true, create general segment expansions
......
......@@ -170,8 +170,7 @@ namespace Nektar
* Given a block matrix, construct a global matrix system according to
* a local to global mapping. #m_linSys is constructed by
* AssembleFullMatrix().
* @param mkey Associated linear system key.
* @param Mat Block matrix.
* @param pkey Associated linear system key.
* @param locToGloMap Local to global mapping.
*/
GlobalLinSys::GlobalLinSys(const GlobalLinSysKey &pKey,
......
......@@ -867,9 +867,13 @@ namespace Nektar
else
{
// Do matrix multiply locally
NekVector<NekDouble> loc(nLocal, m_wsp, eWrapper);
m_locToGloMap->GlobalToLocalBnd(pInput, m_wsp);
#if 1
NekVector<NekDouble> loc(nLocal, m_wsp, eWrapper);
loc = (*m_schurCompl)*loc;
#else
#endif
m_locToGloMap->AssembleBnd(m_wsp, pOutput);
}
}
......
......@@ -71,6 +71,7 @@ namespace Nektar
eSymmetry,
eRinglebFlow,
eTimeDependent,
eRadiation,
eIsentropicVortex,
eCalcBC,
eQinflow,
......@@ -105,6 +106,7 @@ namespace Nektar
known_type["RinglebFlow"] = eRinglebFlow;
known_type["Symmetry"] = eSymmetry;
known_type["TimeDependent"] = eTimeDependent;
known_type["Radiation"] = eRadiation;
known_type["IsentropicVortex"] = eIsentropicVortex;
known_type["InflowCFE"] = eInflowCFE;
known_type["OutflowCFE"] = eOutflowCFE;
......
......@@ -139,9 +139,13 @@ namespace Nektar
{
if (m_fields[0]->GetBndConditions()[n]->GetUserDefined() != SpatialDomains::eTimeDependent)
{
if(m_fields[0]->GetBndConditions()[n]->GetUserDefined() != SpatialDomains::eI)
{
ASSERTL0(false,"Unknown USERDEFINEDTYPE boundary condition");
if(m_fields[0]->GetBndConditions()[n]->GetUserDefined() != SpatialDomains::eRadiation)
{
if(m_fields[0]->GetBndConditions()[n]->GetUserDefined() != SpatialDomains::eI)
{
SpatialDomains::BndUserDefinedType btype = m_fields[0]->GetBndConditions()[n]->GetUserDefined();
ASSERTL0(false,"Unknown USERDEFINEDTYPE boundary condition");
}
}
}
}
......@@ -175,22 +179,71 @@ namespace Nektar
}
m_advObject = GetAdvectionTermFactory().CreateInstance(vConvectiveType, m_session, m_graph);
}
if(m_equationType == eUnsteadyStokes)
{
std::string vConvectiveType = "NoAdvection";
m_advObject = GetAdvectionTermFactory().CreateInstance(vConvectiveType, m_session, m_graph);
}
#if 0 // Not required if building on an UnsteadySystem rather than an EquationSystem
// Set up filters
LibUtilities::FilterMap::const_iterator x;
LibUtilities::FilterMap f = m_session->GetFilters();
for (x = f.begin(); x != f.end(); ++x)
if(m_equationType == eUnsteadyStokes)
{
m_filters.push_back(SolverUtils::GetFilterFactory().CreateInstance(x->first, m_session, x->second));
std::string vConvectiveType = "NoAdvection";
m_advObject = GetAdvectionTermFactory().CreateInstance(vConvectiveType, m_session, m_graph);
}
// check to see if any Robin boundary conditions and if so set
// up m_field to boundary condition maps;
m_fieldsBCToElmtID = Array<OneD, Array<OneD, int> >(m_fields.num_elements());
m_fieldsBCToTraceID = Array<OneD, Array<OneD, int> >(m_fields.num_elements());
m_fieldsRadiationFactor = Array<OneD, Array<OneD, NekDouble> > (m_fields.num_elements());
for (i = 0; i < m_fields.num_elements(); ++i)
{
bool Set = false;
Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
int radpts = 0;
BndConds = m_fields[i]->GetBndConditions();
BndExp = m_fields[i]->GetBndCondExpansions();
for(int n = 0; n < BndConds.num_elements(); ++n)
{
if(BndConds[n]->GetUserDefined() == SpatialDomains::eRadiation)
{
if(Set == false)
{
m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],m_fieldsBCToTraceID[i]);
Set = true;
}
radpts += BndExp[n]->GetTotPoints();
}
}
m_fieldsRadiationFactor[i] = Array<OneD, NekDouble>(radpts);
radpts = 0; // reset to use as a counter
for(int n = 0; n < BndConds.num_elements(); ++n)
{
if(BndConds[n]->GetUserDefined() == SpatialDomains::eRadiation)
{
int npoints = BndExp[n]->GetNpoints();
Array<OneD, NekDouble> x0(npoints,0.0);
Array<OneD, NekDouble> x1(npoints,0.0);
Array<OneD, NekDouble> x2(npoints,0.0);
Array<OneD, NekDouble> tmpArray;
BndExp[n]->GetCoords(x0,x1,x2);
LibUtilities::Equation coeff =
boost::static_pointer_cast<
SpatialDomains::RobinBoundaryCondition
>(BndConds[n])->m_robinPrimitiveCoeff;
coeff.Evaluate(x0,x1,x2,m_time,
tmpArray = m_fieldsRadiationFactor[i]+ radpts);
//Vmath::Neg(npoints,tmpArray = m_fieldsRadiationFactor[i]+ radpts,1);
radpts += npoints;
}
}
}
#endif
}
IncNavierStokes::~IncNavierStokes(void)
......@@ -820,19 +873,86 @@ namespace Nektar
void IncNavierStokes::SetBoundaryConditions(NekDouble time)
{
int nvariables = m_fields.num_elements();
SpatialDomains::BndUserDefinedType BndType;
for (int i = 0; i < nvariables; ++i)
{
for(int n = 0; n < m_fields[i]->GetBndConditions().num_elements(); ++n)
{
if(m_fields[i]->GetBndConditions()[n]->GetUserDefined() == SpatialDomains::eTimeDependent)
if(m_fields[i]->GetBndConditions()[n]->GetUserDefined() ==
SpatialDomains::eTimeDependent)
{
m_fields[i]->EvaluateBoundaryConditions(time);
}
}
SetRadiationBoundaryForcing(i); // Set Radiation conditions if required.
}
}
// Probably should be pushed back into ContField?
void IncNavierStokes::SetRadiationBoundaryForcing(int fieldid)
{
int i,n;
Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
BndConds = m_fields[fieldid]->GetBndConditions();
BndExp = m_fields[fieldid]->GetBndCondExpansions();
StdRegions::StdExpansionSharedPtr elmt;
StdRegions::StdExpansion1DSharedPtr Bc;
int cnt;
int elmtid,nq,offset, boundary;
Array<OneD, NekDouble> Bvals, U;
int cnt1 = 0;
for(cnt = n = 0; n < BndConds.num_elements(); ++n)
{
SpatialDomains::BndUserDefinedType type = BndConds[n]->GetUserDefined();
if((BndConds[n]->GetBoundaryConditionType() == SpatialDomains::eRobin)&&(type == SpatialDomains::eRadiation))
{
for(i = 0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
{
elmtid = m_fieldsBCToElmtID[fieldid][cnt];
elmt = m_fields[fieldid]->GetExp(elmtid);
offset = m_fields[fieldid]->GetPhys_Offset(elmtid);
U = m_fields[fieldid]->UpdatePhys() + offset;
Bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion1D> (BndExp[n]->GetExp(i));
boundary = m_fieldsBCToTraceID[fieldid][cnt];
// Get edge values and put into ubc
nq = Bc->GetTotPoints();
Array<OneD, NekDouble> ubc(nq);
elmt->GetEdgePhysVals(boundary,Bc,U,ubc);
Vmath::Vmul(nq,&m_fieldsRadiationFactor[fieldid][cnt1 + BndExp[n]->GetPhys_Offset(i)],1,&ubc[0],1,&ubc[0],1);
Bvals = BndExp[n]->UpdateCoeffs()+BndExp[n]->GetCoeff_Offset(i);
Bc->IProductWRTBase(ubc,Bvals);
}
cnt1 += BndExp[n]->GetTotPoints();
}
else if(type == SpatialDomains::eNoUserDefined || type == SpatialDomains::eTimeDependent || type == SpatialDomains::eHigh)
{
cnt += BndExp[n]->GetExpSize();
}
else
{
ASSERTL0(false,"Unknown USERDEFINEDTYPE in pressure boundary condition");
}
}
}
// Decide if at a steady state if the discrerte L2 sum of the
// coefficients is the same as the previous step to within the
// tolerance m_steadyStateTol;
......@@ -840,7 +960,7 @@ namespace Nektar
{
static NekDouble previousL2 = 0.0;
bool returnval = false;
NekDouble L2 = 0.0;
// calculate L2 discrete summation
......
......@@ -175,15 +175,16 @@ namespace Nektar
EquationType m_equationType; ///< equation type;
Array<OneD, Array<OneD, int> > m_fieldsBCToElmtID; // Mapping from BCs to Elmt IDs
Array<OneD, Array<OneD, int> > m_fieldsBCToTraceID; // Mapping from BCs to Elmt Edge IDs
Array<OneD, Array<OneD, NekDouble> > m_fieldsRadiationFactor; // RHS Factor for Radiation Condition
// Time integration classes
LibUtilities::TimeIntegrationSchemeOperators m_integrationOps;
Array<OneD, LibUtilities::TimeIntegrationSchemeSharedPtr> m_integrationScheme;
int m_intSteps; ///< Number of time integration steps AND Order of extrapolation for pressure boundary conditions.
// not required if building from an UnsteadySystem
// std::vector<SolverUtils::FilterSharedPtr> m_filters;
/**
* Constructor.
*/
......@@ -211,6 +212,9 @@ namespace Nektar
void SetBoundaryConditions(NekDouble time);
//Set Radiation forcing term
void SetRadiationBoundaryForcing(int fieldid);
// evaluate steady state
bool CalcSteadyState(void);
......
......@@ -409,6 +409,34 @@ namespace Nektar
int nqtot = m_fields[0]->GetTotPoints();
Timer timer;
#if 0
timer.Start();
for(int k = 0; k < 1000; ++k)
{
m_fields[0]->IProductWRTBase(m_fields[0]->GetPhys(),
m_fields[0]->UpdateCoeffs());
}
timer.Stop();
cout << "\t 1000 Iprods : "<< timer.TimePerTest(1) << endl;
#endif
#if 0
timer.Start();
Array<OneD, NekDouble> out (m_fields[0]->GetTotPoints());
Array<OneD, NekDouble> out1(m_fields[0]->GetTotPoints());
Array<OneD, NekDouble> out2(m_fields[0]->GetTotPoints());
for(int k = 0; k < 10000; ++k)
{
m_fields[0]->PhysDeriv(out,out1,out2);
}
timer.Stop();
cout << "\t 10000 Physderiv : "<< timer.TimePerTest(1) << endl;
exit(1);
#endif
timer.Start();
// evaluate convection terms
......@@ -521,6 +549,20 @@ namespace Nektar
// Solve Helmholtz system and put in Physical space
timer.Start();
#if 0
Array<OneD, NekDouble> Save(ncoeffs);
Vmath::Vcopy(ncoeffs,m_fields[0]->GetCoeffs(),1,Save,1);
for(int k = 0; k < 200; ++k)
{
SetBoundaryConditions(time);
m_fields[0]->HelmSolve(F[0], m_fields[0]->UpdateCoeffs(), NullFlagList, factors);
Vmath::Vcopy(ncoeffs,Save,1,m_fields[0]->UpdateCoeffs(),1);
}
timer.Stop();
cout << "\t 200 V_solve : "<< timer.TimePerTest(1) << endl;
exit(1);
#endif
for(i = 0; i < m_nConvectiveFields; ++i)
{
m_fields[i]->HelmSolve(F[i], m_fields[i]->UpdateCoeffs(), NullFlagList, factors);
......@@ -782,7 +824,7 @@ namespace Nektar
}
}
// setting if just standard BC no High order
// setting if just standard BC not High order
else if(type == SpatialDomains::eNoUserDefined || type == SpatialDomains::eTimeDependent)
{
cnt += PBndExp[n]->GetExpSize();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment