38 #include <boost/algorithm/string.hpp>
39 #include <boost/core/ignore_unused.hpp>
68 int shapedim =
m_fields[0]->GetShapeDimension();
70 for (
int j = 0; j < shapedim; ++j)
79 "No TESTTYPE defined in session.");
80 std::string TestTypeStr =
m_session->GetSolverInfo(
"TESTTYPE");
99 std::vector<std::string> vel;
149 std::string riemName;
150 m_session->LoadSolverInfo(
"AdvectionType", advName,
"WeakDG");
154 m_session->LoadSolverInfo(
"UpwindType", riemName,
"Upwind");
167 std::cout <<
"m_vellc is generated with mag = " <<
AvgInt(
m_vellc)
195 ASSERTL0(
false,
"Implicit unsteady Advection not set up.");
217 for (i = 0; i < nfields; ++i)
221 nvariables = nfields;
233 for (i = 0; i < nvariables; ++i)
246 "Only one of IO_CheckTime and IO_CheckSteps "
250 bool doCheckTime =
false;
264 for (
int i = 0; i < nvariables; ++i)
283 std::cout <<
"Steps: " << std::setw(8) << std::left << step + 1
285 <<
"Time: " << std::setw(12) << std::left <<
m_time;
287 std::stringstream ss;
288 ss << cpuTime <<
"s";
289 std::cout <<
" CPU Time: " << std::setw(8) << std::left << ss.str()
297 std::cout <<
"dMass = " << std::setw(8) << std::left << dMass[indx]
304 for (i = 0; i < nvariables; ++i)
324 std::cout <<
"dMass = ";
325 for (i = 0; i < Ntot; ++i)
327 std::cout << dMass[i] <<
" , ";
329 std::cout << std::endl << std::endl;
332 if (
m_session->GetComm()->GetRank() == 0)
338 <<
"CFL time-step : " <<
m_timestep << std::endl;
341 if (
m_session->GetSolverInfo(
"Driver") !=
"SteadyState")
343 std::cout <<
"Time-integration : " << intTime <<
"s" << std::endl;
347 for (i = 0; i < nvariables; ++i)
353 for (i = 0; i < nvariables; ++i)
397 boost::ignore_unused(time);
400 int nvariables = inarray.size();
407 int ncoeffs = inarray[0].size();
414 for (i = 1; i < nvariables; ++i)
416 WeakAdv[i] = WeakAdv[i - 1] + ncoeffs;
422 for (i = 0; i < nvariables; ++i)
424 m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
425 m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
439 for (i = 0; i < nvariables; ++i)
449 ASSERTL0(
false,
"Unknown projection scheme");
469 int nVariables = inarray.size();
484 for (i = 0; i < nVariables; ++i)
486 Vmath::Vcopy(nQuadraturePts, inarray[i], 1, outarray[i], 1);
496 for (i = 0; i < nVariables; ++i)
498 m_fields[i]->FwdTrans(inarray[i], coeffs);
499 m_fields[i]->BwdTrans(coeffs, outarray[i]);
505 ASSERTL0(
false,
"Unknown projection scheme");
522 "Dimension of flux array and velocity array do not match");
525 int nq = physfield[0].size();
527 for (i = 0; i < flux.size(); ++i)
529 for (j = 0; j < flux[0].size(); ++j)
549 for (i = 0; i < nvariables; ++i)
551 physfield[i] = InField[i];
555 for (i = 0; i < nvariables; ++i)
564 m_fields[0]->IProductWRTDirectionalDerivBase(
581 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
598 m_fields[i]->AddFwdBwdTraceIntegral(fluxFwd, fluxBwd, WeakDeriv[j]);
601 Vmath::Vadd(ncoeffs, &WeakDeriv[j][0], 1, &OutField[i][0], 1,
612 NekDouble vel_phi, vel_theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
622 for (
int j = 0; j < nq; j++)
636 -1.0 * vel_phi * sin_varphi - vel_theta * sin_theta * cos_varphi;
638 vel_phi * cos_varphi - vel_theta * sin_theta * sin_varphi;
639 velocity[2][j] = vel_theta * cos_theta;
674 Vmath::Vmul(nq, &tmp[0], 1, &tmp[0], 1, &tmp[0], 1);
679 <<
"Velocity vector is projected onto the tangent plane: Diff = "
868 NekDouble Tol = 0.0001, Maxiter = 1000, N = 100;
869 NekDouble newy, F = 0.0, dF = 1.0, y0, tmp;
880 intval =
sqrt(Rhs - zlevel * zlevel);
886 intval =
sqrt(0.5 * (Rhs - zlevel * zlevel * zlevel * zlevel -
894 (Rhs - zlevel * zlevel * zlevel * zlevel - zlevel * zlevel);
895 intval =
sqrt(0.5 * (1.0 +
sqrt(1.0 + 4.0 * tmp)));
910 for (
int j = 0; j < N + 1; ++j)
912 xp[j] = j * 2.0 * intval / N - intval;
915 for (
int i = 0; i < Maxiter; ++i)
923 F = xp[j] * xp[j] + y0 * y0 + zlevel * zlevel - Rhs;
930 F = 2.0 * xp[j] * xp[j] + y0 * y0 * y0 * y0 +
931 y0 * y0 + zlevel * zlevel * zlevel * zlevel +
932 zlevel * zlevel - Rhs;
933 dF = 4.0 * y0 * y0 * y0 + 2.0 * y0;
943 if (fabs(F / dF) < Tol)
955 "Advection Velocity convergence fails");
964 for (
int j = 0; j < N + 1; ++j)
966 xp[j] = j * 2.0 * intval / N - intval;
968 0.5 * (zlevel * zlevel * zlevel * zlevel +
970 (xp[j] * xp[j] * xp[j] * xp[j] - xp[j] * xp[j]);
987 for (
int j = 0; j < N; ++j)
990 arclength +
sqrt((yp[j + 1] - yp[j]) * (yp[j + 1] - yp[j]) +
991 (xp[j + 1] - xp[j]) * (xp[j + 1] - xp[j])) /
999 bool dumpInitialConditions,
1002 boost::ignore_unused(domain);
1004 int nq =
m_fields[0]->GetNpoints();
1018 for (
int i = 0; i <
m_fields.size(); ++i)
1035 for (
int i = 0; i <
m_fields.size(); ++i)
1050 std::cout <<
"m_Mass0 = " <<
m_Mass0 << std::endl;
1053 for (
int i = 0; i <
m_fields.size(); ++i)
1068 for (
int i = 0; i <
m_fields.size(); ++i)
1081 if (dumpInitialConditions)
1091 int nq =
m_fields[0]->GetNpoints();
1102 m_radius_limit = 0.5;
1106 for (
int j = 0; j < nq; ++j)
1111 dist =
sqrt(x0j * x0j + x1j * x1j);
1113 if (dist < m_radius_limit)
1115 outfield[j] = 0.5 * (1.0 + cos(
m_pi * dist / m_radius_limit));
1126 int nq =
m_fields[0]->GetNpoints();
1128 NekDouble dist, radius, cosdiff, sin_theta, cos_theta, sin_varphi,
1130 NekDouble m_theta_c, m_varphi_c, m_radius_limit, m_c0;
1140 m_varphi_c = 3.0 *
m_pi / 2.0;
1141 m_radius_limit = 7.0 *
m_pi / 64.0;
1146 for (
int j = 0; j < nq; ++j)
1152 radius =
sqrt(x0j * x0j + x1j * x1j + x2j * x2j);
1154 sin_varphi = x1j /
sqrt(x0j * x0j + x1j * x1j);
1155 cos_varphi = x0j /
sqrt(x0j * x0j + x1j * x1j);
1157 sin_theta = x2j / radius;
1158 cos_theta =
sqrt(x0j * x0j + x1j * x1j) / radius;
1160 cosdiff = cos_varphi * cos(m_varphi_c) + sin_varphi * sin(m_varphi_c);
1161 dist = radius * acos(sin(m_theta_c) * sin_theta +
1162 cos(m_theta_c) * cos_theta * cosdiff);
1164 if (dist < m_radius_limit)
1167 0.5 * (1.0 + cos(
m_pi * dist / m_radius_limit)) + m_c0;
1179 int nq =
m_fields[0]->GetNpoints();
1185 m_fields[0]->GetCoords(x0, x1, x2);
1189 for (
int i = 0; i < nq; ++i)
1201 int nq =
m_fields[0]->GetNpoints();
1207 m_fields[0]->GetCoords(x0, x1, x2);
1211 for (
int i = 0; i < nq; ++i)
1223 int nq =
m_fields[0]->GetNpoints();
1245 1, &velcoeff[0], 1, &velcoeff[0], 1);
1249 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
1259 1, &vellc[0], 1, &vellc[0], 1);
1266 1, &vellc[0], 1, &vellc[0], 1);
1273 1, &vellc[0], 1, &vellc[0], 1);
1284 int nq =
m_fields[0]->GetNpoints();
1296 1, &veldotMF[j][0], 1, &veldotMF[j][0], 1);
1305 boost::ignore_unused(field);
#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)
virtual void v_DoSolve()
Solves an unsteady problem.
virtual ~MMFAdvection()
Destructor.
Array< OneD, Array< OneD, NekDouble > > m_veldotMF
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)
Session reader.
void ComputeNablaCdotVelocity(Array< OneD, NekDouble > &vellc)
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print Summary.
void EvaluateAdvectionVelocity(Array< OneD, Array< OneD, NekDouble >> &velocity)
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the RHS.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the projection.
virtual void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
void WeakDGDirectionalAdvection(const Array< OneD, Array< OneD, NekDouble >> &InField, Array< OneD, 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.
void ComputeveldotMF(Array< OneD, Array< OneD, NekDouble >> &veldotMF)
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
Array< OneD, Array< OneD, NekDouble > > m_velocity
Advection velocity.
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity.
virtual void v_SetInitialConditions(const NekDouble initialtime, bool dumpInitialConditions, const int domain)
static std::string className
Name of class.
Array< OneD, NekDouble > m_vellc
void AdvectionBellPlane(Array< OneD, NekDouble > &outfield)
virtual void v_InitObject(bool DeclareFields=true)
Initialise the object.
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).
SOLVER_UTILS_EXPORT int GetTraceNpoints()
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
NekDouble m_fintime
Finish time of the simulation.
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.
SOLVER_UTILS_EXPORT int GetNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
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 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)
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 MMFInitObject(const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem=-1)
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
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
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.
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
std::vector< int > m_intVariables
int m_infosteps
Number of time steps between outputting status information.
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
The above copyright notice and this permission notice shall be included.
@ 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)