38#include <boost/algorithm/string.hpp>
67 int shapedim =
m_fields[0]->GetShapeDimension();
69 for (
int j = 0; j < shapedim; ++j)
78 "No TESTTYPE defined in session.");
79 std::string TestTypeStr =
m_session->GetSolverInfo(
"TESTTYPE");
98 std::vector<std::string> vel;
147 std::string riemName;
148 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
152 m_session->LoadSolverInfo(
"UpwindType", riemName,
"Upwind");
165 std::cout <<
"m_vellc is generated with mag = " <<
AvgInt(
m_vellc)
193 ASSERTL0(
false,
"Implicit unsteady Advection not set up.");
208 for (i = 0; i < nfields; ++i)
212 nvariables = nfields;
224 for (i = 0; i < nvariables; ++i)
237 "Only one of IO_CheckTime and IO_CheckSteps "
241 bool doCheckTime =
false;
255 for (
int i = 0; i < nvariables; ++i)
275 std::cout <<
"Steps: " << std::setw(8) << std::left << step + 1
277 <<
"Time: " << std::setw(12) << std::left <<
m_time;
279 std::stringstream ss;
280 ss << cpuTime <<
"s";
281 std::cout <<
" CPU Time: " << std::setw(8) << std::left << ss.str()
289 std::cout <<
"dMass = " << std::setw(8) << std::left << dMass[indx]
296 for (i = 0; i < nvariables; ++i)
316 std::cout <<
"dMass = ";
317 for (i = 0; i < Ntot; ++i)
319 std::cout << dMass[i] <<
" , ";
321 std::cout << std::endl << std::endl;
324 if (
m_session->GetComm()->GetRank() == 0)
330 <<
"CFL time-step : " <<
m_timestep << std::endl;
333 if (
m_session->GetSolverInfo(
"Driver") !=
"SteadyState")
335 std::cout <<
"Time-integration : " << intTime <<
"s" << std::endl;
339 for (i = 0; i < nvariables; ++i)
345 for (i = 0; i < nvariables; ++i)
391 int nvariables = inarray.size();
398 int ncoeffs = inarray[0].size();
405 for (i = 1; i < nvariables; ++i)
407 WeakAdv[i] = WeakAdv[i - 1] + ncoeffs;
413 for (i = 0; i < nvariables; ++i)
415 m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
416 m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
425 for (i = 0; i < nvariables; ++i)
435 ASSERTL0(
false,
"Unknown projection scheme");
456 int nVariables = inarray.size();
471 if (inarray != outarray)
473 for (i = 0; i < nVariables; ++i)
475 Vmath::Vcopy(nQuadraturePts, inarray[i], 1, outarray[i], 1);
486 for (i = 0; i < nVariables; ++i)
488 m_fields[i]->FwdTrans(inarray[i], coeffs);
489 m_fields[i]->BwdTrans(coeffs, outarray[i]);
495 ASSERTL0(
false,
"Unknown projection scheme");
512 "Dimension of flux array and velocity array do not match");
515 int nq = physfield[0].size();
517 for (i = 0; i < flux.size(); ++i)
519 for (j = 0; j < flux[0].size(); ++j)
539 for (i = 0; i < nvariables; ++i)
541 physfield[i] = InField[i];
545 for (i = 0; i < nvariables; ++i)
554 m_fields[0]->IProductWRTDirectionalDerivBase(
571 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
588 m_fields[i]->AddFwdBwdTraceIntegral(fluxFwd, fluxBwd, WeakDeriv[j]);
591 Vmath::Vadd(ncoeffs, &WeakDeriv[j][0], 1, &OutField[i][0], 1,
602 NekDouble vel_phi, vel_theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
612 for (
int j = 0; j < nq; j++)
626 -1.0 * vel_phi * sin_varphi - vel_theta * sin_theta * cos_varphi;
628 vel_phi * cos_varphi - vel_theta * sin_theta * sin_varphi;
629 velocity[2][j] = vel_theta * cos_theta;
664 Vmath::Vmul(nq, &tmp[0], 1, &tmp[0], 1, &tmp[0], 1);
669 <<
"Velocity vector is projected onto the tangent plane: Diff = "
683 NekDouble Tol = 0.0001, Maxiter = 1000, N = 100;
684 NekDouble newy, F = 0.0, dF = 1.0, y0, tmp;
695 intval =
sqrt(Rhs - zlevel * zlevel);
701 intval =
sqrt(0.5 * (Rhs - zlevel * zlevel * zlevel * zlevel -
709 (Rhs - zlevel * zlevel * zlevel * zlevel - zlevel * zlevel);
710 intval =
sqrt(0.5 * (1.0 +
sqrt(1.0 + 4.0 * tmp)));
725 for (
int j = 0; j < N + 1; ++j)
727 xp[j] = j * 2.0 * intval / N - intval;
730 for (
int i = 0; i < Maxiter; ++i)
738 F = xp[j] * xp[j] + y0 * y0 + zlevel * zlevel - Rhs;
745 F = 2.0 * xp[j] * xp[j] + y0 * y0 * y0 * y0 +
746 y0 * y0 + zlevel * zlevel * zlevel * zlevel +
747 zlevel * zlevel - Rhs;
748 dF = 4.0 * y0 * y0 * y0 + 2.0 * y0;
758 if (fabs(F / dF) < Tol)
770 "Advection Velocity convergence fails");
779 for (
int j = 0; j < N + 1; ++j)
781 xp[j] = j * 2.0 * intval / N - intval;
783 0.5 * (zlevel * zlevel * zlevel * zlevel +
785 (xp[j] * xp[j] * xp[j] * xp[j] - xp[j] * xp[j]);
802 for (
int j = 0; j < N; ++j)
805 arclength +
sqrt((yp[j + 1] - yp[j]) * (yp[j + 1] - yp[j]) +
806 (xp[j + 1] - xp[j]) * (xp[j + 1] - xp[j])) /
814 bool dumpInitialConditions,
815 [[maybe_unused]]
const int domain)
831 for (
int i = 0; i <
m_fields.size(); ++i)
848 for (
int i = 0; i <
m_fields.size(); ++i)
863 std::cout <<
"m_Mass0 = " <<
m_Mass0 << std::endl;
866 for (
int i = 0; i <
m_fields.size(); ++i)
881 for (
int i = 0; i <
m_fields.size(); ++i)
894 if (dumpInitialConditions)
915 m_radius_limit = 0.5;
919 for (
int j = 0; j < nq; ++j)
924 dist =
sqrt(x0j * x0j + x1j * x1j);
926 if (dist < m_radius_limit)
928 outfield[j] = 0.5 * (1.0 + cos(
m_pi * dist / m_radius_limit));
941 NekDouble dist, radius, cosdiff, sin_theta, cos_theta, sin_varphi,
943 NekDouble m_theta_c, m_varphi_c, m_radius_limit, m_c0;
953 m_varphi_c = 3.0 *
m_pi / 2.0;
954 m_radius_limit = 7.0 *
m_pi / 64.0;
959 for (
int j = 0; j < nq; ++j)
965 radius =
sqrt(x0j * x0j + x1j * x1j + x2j * x2j);
967 sin_varphi = x1j /
sqrt(x0j * x0j + x1j * x1j);
968 cos_varphi = x0j /
sqrt(x0j * x0j + x1j * x1j);
970 sin_theta = x2j / radius;
971 cos_theta =
sqrt(x0j * x0j + x1j * x1j) / radius;
973 cosdiff = cos_varphi * cos(m_varphi_c) + sin_varphi * sin(m_varphi_c);
974 dist = radius * acos(sin(m_theta_c) * sin_theta +
975 cos(m_theta_c) * cos_theta * cosdiff);
977 if (dist < m_radius_limit)
980 0.5 * (1.0 + cos(
m_pi * dist / m_radius_limit)) + m_c0;
1002 for (
int i = 0; i < nq; ++i)
1014 int nq =
m_fields[0]->GetNpoints();
1020 m_fields[0]->GetCoords(x0, x1, x2);
1024 for (
int i = 0; i < nq; ++i)
1036 int nq =
m_fields[0]->GetNpoints();
1058 1, &velcoeff[0], 1, &velcoeff[0], 1);
1062 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
1072 1, &vellc[0], 1, &vellc[0], 1);
1079 1, &vellc[0], 1, &vellc[0], 1);
1086 1, &vellc[0], 1, &vellc[0], 1);
1097 int nq =
m_fields[0]->GetNpoints();
1109 1, &veldotMF[j][0], 1, &veldotMF[j][0], 1);
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
void Test2Dproblem(const NekDouble time, Array< OneD, NekDouble > &outfield)
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
Evaluate the flux at each solution point.
Array< OneD, Array< OneD, NekDouble > > m_veldotMF
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
void Test3Dproblem(const NekDouble time, Array< OneD, NekDouble > &outfield)
void AdvectionBellSphere(Array< OneD, NekDouble > &outfield)
NekDouble ComputeCirculatingArclength(const NekDouble zlevel, const NekDouble Rhs)
MMFAdvection(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void WeakDGDirectionalAdvection(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField)
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
void ComputeNablaCdotVelocity(Array< OneD, NekDouble > &vellc)
void ComputeveldotMF(Array< OneD, Array< OneD, NekDouble > > &veldotMF)
void v_SetInitialConditions(const NekDouble initialtime, bool dumpInitialConditions, const int domain) override
void EvaluateAdvectionVelocity(Array< OneD, Array< OneD, NekDouble > > &velocity)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print Summary.
void v_InitObject(bool DeclareFields=true) override
Initialise the object.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
Array< OneD, Array< OneD, NekDouble > > m_velocity
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity.
void v_DoSolve() override
Solves an unsteady problem.
static std::string className
Name of class.
Array< OneD, NekDouble > m_vellc
void AdvectionBellPlane(Array< OneD, NekDouble > &outfield)
Array< OneD, NekDouble > m_traceVn
A base class for PDEs which include an advection component.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
int m_spacedim
Spatial dimension (>= expansion dim).
NekDouble m_timestep
Time step size.
int m_infosteps
Number of time steps between outputting status information.
NekDouble m_time
Current time of simulation.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
NekDouble m_fintime
Finish time of the simulation.
SOLVER_UTILS_EXPORT int GetNcoeffs()
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
NekDouble m_checktime
Time between checkpoints.
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
std::string m_sessionName
Name of the session.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
int m_steps
Number of steps to take.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
int m_checksteps
Number of steps between checkpoints.
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
A base class for PDEs which include an advection component.
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Virtual function for generating summary information.
SOLVER_UTILS_EXPORT void VectorCrossProd(const Array< OneD, const Array< OneD, NekDouble > > &v1, const Array< OneD, const Array< OneD, NekDouble > > &v2, Array< OneD, Array< OneD, NekDouble > > &v3)
SOLVER_UTILS_EXPORT NekDouble AvgInt(const Array< OneD, const NekDouble > &inarray)
Array< OneD, Array< OneD, NekDouble > > m_movingframes
SOLVER_UTILS_EXPORT void GramSchumitz(const Array< OneD, const Array< OneD, NekDouble > > &v1, const Array< OneD, const Array< OneD, NekDouble > > &v2, Array< OneD, Array< OneD, NekDouble > > &outarray, bool KeepTheMagnitude=true)
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
SOLVER_UTILS_EXPORT void ComputencdotMF()
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
SOLVER_UTILS_EXPORT NekDouble RootMeanSquare(const Array< OneD, const NekDouble > &inarray)
SurfaceType m_surfaceType
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
std::vector< int > m_intVariables
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ eMixed_CG_Discontinuous
static const NekDouble kNekZeroTol
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
std::vector< std::pair< std::string, std::string > > SummaryList
EquationSystemFactory & GetEquationSystemFactory()
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::vector< double > z(NPUPPER)
@ SIZE_TestType
Length of enum list.
const char *const TestTypeMap[]
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Neg(int n, T *x, const int incx)
Negate x = -x.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
scalarT< T > sqrt(scalarT< T > in)