#include <MMFSWE.h>
Inheritance diagram for Nektar::MMFSWE:
Static Public Member Functions | |
| static SolverUtils::EquationSystemSharedPtr | create (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph) |
| Creates an instance of this class. More... | |
Public Attributes | |
| TestType | m_TestType |
 Public Attributes inherited from Nektar::SolverUtils::MMFSystem | |
| NekDouble | m_pi |
| int | m_shapedim |
| SurfaceType | m_surfaceType |
| UpwindType | m_upwindType |
| TestMaxwellType | m_TestMaxwellType |
| PolType | m_PolType |
| IncType | m_IncType |
| Array< OneD, NekDouble > | m_MMFfactors |
Static Public Attributes | |
| static std::string | className |
| Name of class. More... | |
| Static Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem | |
| static std::string | cmdSetStartTime |
| static std::string | cmdSetStartChkNum |
Protected Attributes | |
| Array< OneD, NekDouble > | m_depth |
| Still water depth. More... | |
| Array< OneD, Array< OneD, NekDouble > > | m_Derivdepth |
| Array< OneD, NekDouble > | m_coriolis |
| Coriolis force. More... | |
| int | m_AddCoriolis |
| int | m_AddRotation |
| NekDouble | m_g |
| NekDouble | m_alpha |
| NekDouble | m_u0 |
| NekDouble | m_Omega |
| NekDouble | m_H0 |
| NekDouble | m_Hvar |
| NekDouble | m_k2 |
| NekDouble | m_hs0 |
| NekDouble | m_uthetamax |
| NekDouble | m_theta0 |
| NekDouble | m_theta1 |
| NekDouble | m_en |
| NekDouble | m_hbar |
| NekDouble | m_angfreq |
| NekDouble | m_K |
| NekDouble | m_Vorticity0 |
| NekDouble | m_Mass0 |
| NekDouble | m_Energy0 |
| NekDouble | m_Enstrophy0 |
| bool | m_primitive |
| Indicates if variables are primitive or conservative. More... | |
| Array< OneD, Array< OneD, NekDouble > > | m_velocity |
| Advection velocity. More... | |
| Array< OneD, NekDouble > | m_traceVn |
| Array< OneD, Array< OneD, NekDouble > > | m_veldotMF |
| Array< OneD, NekDouble > | m_vellc |
| int | m_planeNumber |
| Protected Attributes inherited from Nektar::SolverUtils::MMFSystem | |
| NekDouble | m_alpha |
| NekDouble | m_Incfreq |
| int | m_SmoothFactor |
| NekDouble | m_SFinit |
| Array< OneD, Array< OneD, NekDouble > > | m_movingframes |
| Array< OneD, Array< OneD, NekDouble > > | m_surfaceNormal |
| Array< OneD, Array< OneD, NekDouble > > | m_ncdotMFFwd |
| Array< OneD, Array< OneD, NekDouble > > | m_ncdotMFBwd |
| Array< OneD, Array< OneD, NekDouble > > | m_nperpcdotMFFwd |
| Array< OneD, Array< OneD, NekDouble > > | m_nperpcdotMFBwd |
| Array< OneD, Array< OneD, NekDouble > > | m_DivMF |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_CurlMF |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_MFtraceFwd |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_MFtraceBwd |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_ntimesMFFwd |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_ntimesMFBwd |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_ntimes_ntimesMFFwd |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_ntimes_ntimesMFBwd |
| Array< OneD, Array< OneD, NekDouble > > | m_ZimFwd |
| Array< OneD, Array< OneD, NekDouble > > | m_ZimBwd |
| Array< OneD, Array< OneD, NekDouble > > | m_YimFwd |
| Array< OneD, Array< OneD, NekDouble > > | m_YimBwd |
| Array< OneD, Array< OneD, NekDouble > > | m_epsvec |
| Array< OneD, Array< OneD, NekDouble > > | m_muvec |
| Array< OneD, Array< OneD, NekDouble > > | m_negepsvecminus1 |
| Array< OneD, Array< OneD, NekDouble > > | m_negmuvecminus1 |
| Array< OneD, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > > | m_dedxi_cdot_e |
| SpatialDomains::GeomMMF | m_MMFdir |
| Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem | |
| LibUtilities::TimeIntegrationSchemeSharedPtr | m_intScheme |
| Wrapper to the time integration scheme. More... | |
| LibUtilities::TimeIntegrationSchemeOperators | m_ode |
| The time integration scheme operators to use. More... | |
| Array< OneD, Array< OneD, NekDouble > > | m_previousSolution |
| Storage for previous solution for steady-state check. More... | |
| std::vector< int > | m_intVariables |
| NekDouble | m_cflSafetyFactor |
| CFL safety factor (comprise between 0 to 1). More... | |
| NekDouble | m_CFLGrowth |
| CFL growth rate. More... | |
| NekDouble | m_CFLEnd |
| Maximun cfl in cfl growth. More... | |
| int | m_abortSteps |
| Number of steps between checks for abort conditions. More... | |
| bool | m_explicitDiffusion |
| Indicates if explicit or implicit treatment of diffusion is used. More... | |
| bool | m_explicitAdvection |
| Indicates if explicit or implicit treatment of advection is used. More... | |
| bool | m_explicitReaction |
| Indicates if explicit or implicit treatment of reaction is used. More... | |
| int | m_steadyStateSteps |
| Check for steady state at step interval. More... | |
| NekDouble | m_steadyStateTol |
| Tolerance to which steady state should be evaluated at. More... | |
| int | m_filtersInfosteps |
| Number of time steps between outputting filters information. More... | |
| std::vector< std::pair< std::string, FilterSharedPtr > > | m_filters |
| bool | m_homoInitialFwd |
| Flag to determine if simulation should start in homogeneous forward transformed state. More... | |
| std::ofstream | m_errFile |
| NekDouble | m_epsilon |
| Diffusion coefficient. More... | |
| Protected Attributes inherited from Nektar::SolverUtils::EquationSystem | |
| LibUtilities::CommSharedPtr | m_comm |
| Communicator. More... | |
| bool | m_verbose |
| LibUtilities::SessionReaderSharedPtr | m_session |
| The session reader. More... | |
| std::map< std::string, SolverUtils::SessionFunctionSharedPtr > | m_sessionFunctions |
| Map of known SessionFunctions. More... | |
| LibUtilities::FieldIOSharedPtr | m_fld |
| Field input/output. More... | |
| Array< OneD, MultiRegions::ExpListSharedPtr > | m_fields |
| Array holding all dependent variables. More... | |
| SpatialDomains::BoundaryConditionsSharedPtr | m_boundaryConditions |
| Pointer to boundary conditions object. More... | |
| SpatialDomains::MeshGraphSharedPtr | m_graph |
| Pointer to graph defining mesh. More... | |
| std::string | m_sessionName |
| Name of the session. More... | |
| NekDouble | m_time |
| Current time of simulation. More... | |
| int | m_initialStep |
| Number of the step where the simulation should begin. More... | |
| NekDouble | m_fintime |
| Finish time of the simulation. More... | |
| NekDouble | m_timestep |
| Time step size. More... | |
| NekDouble | m_lambda |
| Lambda constant in real system if one required. More... | |
| NekDouble | m_checktime |
| Time between checkpoints. More... | |
| NekDouble | m_lastCheckTime |
| NekDouble | m_TimeIncrementFactor |
| int | m_nchk |
| Number of checkpoints written so far. More... | |
| int | m_steps |
| Number of steps to take. More... | |
| int | m_checksteps |
| Number of steps between checkpoints. More... | |
| int | m_infosteps |
| Number of time steps between outputting status information. More... | |
| int | m_iterPIT = 0 |
| Number of parallel-in-time time iteration. More... | |
| int | m_windowPIT = 0 |
| Index of windows for parallel-in-time time iteration. More... | |
| int | m_spacedim |
| Spatial dimension (>= expansion dim). More... | |
| int | m_expdim |
| Expansion dimension. More... | |
| bool | m_singleMode |
| Flag to determine if single homogeneous mode is used. More... | |
| bool | m_halfMode |
| Flag to determine if half homogeneous mode is used. More... | |
| bool | m_multipleModes |
| Flag to determine if use multiple homogenenous modes are used. More... | |
| bool | m_useFFT |
| Flag to determine if FFT is used for homogeneous transform. More... | |
| bool | m_homogen_dealiasing |
| Flag to determine if dealiasing is used for homogeneous simulations. More... | |
| bool | m_specHP_dealiasing |
| Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More... | |
| enum MultiRegions::ProjectionType | m_projectionType |
| Type of projection; e.g continuous or discontinuous. More... | |
| Array< OneD, Array< OneD, NekDouble > > | m_traceNormals |
| Array holding trace normals for DG simulations in the forwards direction. More... | |
| Array< OneD, bool > | m_checkIfSystemSingular |
| Flag to indicate if the fields should be checked for singularity. More... | |
| LibUtilities::FieldMetaDataMap | m_fieldMetaDataMap |
| Map to identify relevant solver info to dump in output fields. More... | |
| Array< OneD, NekDouble > | m_movingFrameVelsxyz |
| Moving frame of reference velocities. More... | |
| Array< OneD, NekDouble > | m_movingFrameTheta |
| Moving frame of reference angles with respect to the. More... | |
| boost::numeric::ublas::matrix< NekDouble > | m_movingFrameProjMat |
| Projection matrix for transformation between inertial and moving. More... | |
| int | m_NumQuadPointsError |
| Number of Quadrature points used to work out the error. More... | |
| enum HomogeneousType | m_HomogeneousType |
| NekDouble | m_LhomX |
| physical length in X direction (if homogeneous) More... | |
| NekDouble | m_LhomY |
| physical length in Y direction (if homogeneous) More... | |
| NekDouble | m_LhomZ |
| physical length in Z direction (if homogeneous) More... | |
| int | m_npointsX |
| number of points in X direction (if homogeneous) More... | |
| int | m_npointsY |
| number of points in Y direction (if homogeneous) More... | |
| int | m_npointsZ |
| number of points in Z direction (if homogeneous) More... | |
| int | m_HomoDirec |
| number of homogenous directions More... | |
Private Member Functions | |
| void | TestSWE2Dproblem (const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield) |
| void | Checkpoint_Output_Cartesian (std::string outname) |
Private Attributes | |
| int | m_RossbyDisturbance |
| int | m_PurturbedJet |
Friends | |
| class | MemoryManager< MMFSWE > |
Additional Inherited Members | |
| Protected Types inherited from Nektar::SolverUtils::EquationSystem | |
| enum | HomogeneousType { eHomogeneous1D , eHomogeneous2D , eHomogeneous3D , eNotHomogeneous } |
| Parameter for homogeneous expansions. More... | |
| Static Protected Attributes inherited from Nektar::SolverUtils::EquationSystem | |
| static std::string | equationSystemTypeLookupIds [] |
| static std::string | projectionTypeLookupIds [] |
Detailed Description
Constructor & Destructor Documentation
◆ ~MMFSWE()
|
virtual |
Destructor.
Unsteady linear advection equation destructor.
Definition at line 245 of file MMFSWE.cpp.
246{
247}
◆ MMFSWE()
|
protected |
Session reader.
Definition at line 52 of file MMFSWE.cpp.
55{
56 m_planeNumber = 0;
57}
SOLVER_UTILS_EXPORT MMFSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition: MMFSystem.cpp:43
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.
Definition: UnsteadySystem.cpp:76
References m_planeNumber.
Member Function Documentation
◆ AddCoriolis()
|
protected |
Definition at line 1201 of file MMFSWE.cpp.
1203{
1204 int ncoeffs = outarray[0].size();
1205 int nq = physarray[0].size();
1206
1207 Array<OneD, NekDouble> h(nq);
1208 Array<OneD, NekDouble> tmp(nq);
1209 Array<OneD, NekDouble> tmpc(ncoeffs);
1210
1211 // physarray is primitive
1212 // conservative formulation compute h
1213 // h = \eta + d
1215
1216 int indx = 0;
1218 {
1219 if (j == 0)
1220 {
1221 indx = 2;
1222 }
1223
1224 else if (j == 1)
1225 {
1226 indx = 1;
1227 }
1228
1229 // add to hu equation
1231 Vmath::Vmul(nq, h, 1, tmp, 1, tmp, 1);
1232
1233 if (j == 1)
1234 {
1235 Vmath::Neg(nq, tmp, 1);
1236 }
1237
1238 // N \cdot (e^1 \times e^2 )
1239 // Vmath::Vmul(nq, &m_MF1crossMF2dotSN[0], 1, &tmp[0], 1, &tmp[0], 1);
1240 m_fields[0]->IProductWRTBase(tmp, tmpc);
1241 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1242 }
1243}
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:363
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.
Definition: Vmath.cpp:207
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.
Definition: Vmath.cpp:354
References m_coriolis, m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Neg(), Vmath::Vadd(), and Vmath::Vmul().
Referenced by DoOdeRhs().
◆ AddDivForGradient()
|
protected |
Definition at line 552 of file MMFSWE.cpp.
554{
555 // routine works for both primitive and conservative formulations
556 int ncoeffs = outarray[0].size();
557 int nq = physarray[0].size();
558
559 Array<OneD, NekDouble> h(nq);
560 Array<OneD, NekDouble> tmp(nq);
561 Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
563 {
564 fluxvector[i] = Array<OneD, NekDouble>(nq);
565 }
566
567 // Get 0.5 g H*H / || e^m ||^2
568 GetSWEFluxVector(3, physarray, fluxvector);
569
570 Array<OneD, NekDouble> tmpc(ncoeffs);
572 {
574
575 Vmath::Neg(nq, &tmp[0], 1);
576 m_fields[0]->IProductWRTBase(tmp, tmpc);
577
578 Vmath::Vadd(ncoeffs, outarray[j + 1], 1, tmpc, 1, outarray[j + 1], 1);
579 }
580}
void GetSWEFluxVector(const int i, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: MMFSWE.cpp:582
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Definition: MMFSystem.h:194
References GetSWEFluxVector(), Nektar::SolverUtils::MMFSystem::m_DivMF, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Neg(), Vmath::Vadd(), and Vmath::Vmul().
Referenced by DoOdeRhs().
◆ AddElevationEffect()
|
protected |
Definition at line 1246 of file MMFSWE.cpp.
1248{
1249 int ncoeffs = outarray[0].size();
1250 int nq = physarray[0].size();
1251
1252 Array<OneD, NekDouble> h(nq);
1253 Array<OneD, NekDouble> tmp(nq);
1254 Array<OneD, NekDouble> tmpc(ncoeffs);
1255
1256 // physarray is primitive
1257 // conservative formulation compute h
1258 // h = \eta + d
1260
1262 {
1265
1266 m_fields[0]->IProductWRTBase(tmp, tmpc);
1267
1268 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1269 }
1270}
Array< OneD, Array< OneD, NekDouble > > m_Derivdepth
Definition: MMFSWE.h:87
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:245
References m_depth, m_Derivdepth, Nektar::SolverUtils::EquationSystem::m_fields, m_g, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Smul(), Vmath::Vadd(), and Vmath::Vmul().
Referenced by DoOdeRhs().
◆ AddRotation()
|
protected |
Definition at line 1275 of file MMFSWE.cpp.
1277{
1278 // routine works for both primitive and conservative formulations
1279 int ncoeffs = outarray[0].size();
1280 int nq = physarray[0].size();
1281
1282 // Compute h
1283 Array<OneD, NekDouble> h(nq);
1285
1286 Array<OneD, NekDouble> de0dt_cdot_e0;
1287 Array<OneD, NekDouble> de0dt_cdot_e1;
1288 Array<OneD, NekDouble> de1dt_cdot_e0;
1289 Array<OneD, NekDouble> de1dt_cdot_e1;
1290 Compute_demdt_cdot_ek(0, 0, physarray, de0dt_cdot_e0);
1291 Compute_demdt_cdot_ek(1, 0, physarray, de1dt_cdot_e0);
1292 Compute_demdt_cdot_ek(0, 1, physarray, de0dt_cdot_e1);
1293 Compute_demdt_cdot_ek(1, 1, physarray, de1dt_cdot_e1);
1294
1295 Array<OneD, NekDouble> Rott1(nq);
1296 Array<OneD, NekDouble> Rott2(nq);
1297 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e0, 1, Rott1, 1);
1298 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e1, 1, Rott2, 1);
1299 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e0, 1, Rott1, 1, Rott1, 1);
1300 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e1, 1, Rott2, 1, Rott2, 1);
1301
1302 // Multiply H and \partial \phi / \partial t which is assumed to be u_{\phi}
1303 Vmath::Vmul(nq, &h[0], 1, &Rott1[0], 1, &Rott1[0], 1);
1304 Vmath::Vmul(nq, &h[0], 1, &Rott2[0], 1, &Rott2[0], 1);
1305
1306 Vmath::Neg(nq, Rott1, 1);
1307 Vmath::Neg(nq, Rott2, 1);
1308
1309 Array<OneD, NekDouble> tmpc1(ncoeffs);
1310 Array<OneD, NekDouble> tmpc2(ncoeffs);
1311 m_fields[0]->IProductWRTBase(Rott1, tmpc1);
1312 m_fields[0]->IProductWRTBase(Rott2, tmpc2);
1313
1314 Vmath::Vadd(ncoeffs, tmpc1, 1, outarray[1], 1, outarray[1], 1);
1315 Vmath::Vadd(ncoeffs, tmpc2, 1, outarray[2], 1, outarray[2], 1);
1316}
void Compute_demdt_cdot_ek(const int indm, const int indk, const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1322
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
Definition: Vmath.cpp:569
References Compute_demdt_cdot_ek(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Neg(), Vmath::Vadd(), Vmath::Vmul(), and Vmath::Vvtvp().
Referenced by DoOdeRhs().
◆ AverageFlux()
|
protected |
Definition at line 960 of file MMFSWE.cpp.
964{
965 NekDouble MageF1, MageF2, MageB1, MageB2;
966 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
967 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
968
970 NekDouble uRF, vRF, uLB, vLB;
971
972 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
973 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
974
975 // uRF = uR component in moving frames e^{Fwd}
976 // vRF = vR component in moving frames e^{Fwd}
977 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
978 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
979
980 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
981 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
982 numfluxF[2] = 0.0;
983
984 numfluxF[3] =
985 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
986 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
987 numfluxF[5] = 0.0;
988
989 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
990 numfluxF[7] =
991 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
992 numfluxF[8] = 0.0;
993
994 // uLB = uL component in moving frames e^{Bwd}
995 // vLB = vL component in moving frames e^{Bwd}
996 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
997 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
998
999 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1000 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1001 numfluxB[2] = 0.0;
1002
1003 numfluxB[3] =
1004 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1005 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1006 numfluxB[5] = 0.0;
1007
1008 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1009 numfluxB[7] =
1010 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1011 numfluxB[8] = 0.0;
1012}
void ComputeMagAndDot(const int index, NekDouble &MageF1, NekDouble &MageF2, NekDouble &MageB1, NekDouble &MageB2, NekDouble &eF1_cdot_eB1, NekDouble &eF1_cdot_eB2, NekDouble &eF2_cdot_eB1, NekDouble &eF2_cdot_eB2)
Definition: MMFSWE.cpp:1163
References ComputeMagAndDot(), and m_g.
Referenced by NumericalSWEFlux().
◆ Checkpoint_Output_Cartesian()
|
private |
Definition at line 2767 of file MMFSWE.cpp.
2768{
2771
2772 NekDouble rad_earth = 6.37122 * 1000000;
2773
2774 int nvariables = 7;
2775
2776 // vector u in Cartesian coordinates
2777 std::vector<std::string> variables(nvariables);
2778
2779 variables[0] = "eta";
2780 variables[1] = "hstar";
2781 variables[2] = "vorticity";
2782 variables[3] = "ux";
2783 variables[4] = "uy";
2784 variables[5] = "uz";
2785 variables[6] = "null";
2786
2787 // Obtain \vec{u} in cartesian coordinate
2788 Array<OneD, Array<OneD, NekDouble>> fieldphys(nvariables);
2789 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvariables);
2790 for (int i = 0; i < nvariables; ++i)
2791 {
2792 fieldphys[i] = Array<OneD, NekDouble>(nq, 0.0);
2793 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2794 }
2795
2797 &fieldphys[0][0], 1);
2798 // Vmath::Vadd(nq, &(m_fields[0]->GetPhys())[0], 1, &m_depth[0], 1,
2799 // &fieldphys[1][0], 1);
2800
2802 Vmath::Neg(nq, &fieldphys[1][0], 1);
2804 Vmath::Smul(nq, rad_earth, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2805
2806 Array<OneD, NekDouble> utmp(nq);
2807 Array<OneD, NekDouble> vtmp(nq);
2808
2811
2812 ComputeVorticity(utmp, vtmp, fieldphys[2]);
2813 // u_x = u e^1_x + v e^2_x
2814 int indx = 3;
2816 {
2818 &fieldphys[k + indx][0], 1);
2820 &fieldphys[k + indx][0], 1, &fieldphys[k + indx][0], 1);
2821 }
2822
2823 for (int i = 0; i < nvariables; ++i)
2824 {
2825 m_fields[0]->FwdTrans(fieldphys[i], fieldcoeffs[i]);
2826 }
2827
2829}
void ComputeVorticity(const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
Definition: MMFSWE.cpp:2211
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1249
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:185
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add scalar y = alpha + x.
Definition: Vmath.cpp:379
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1191
References ComputeVorticity(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, m_H0, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Neg(), Vmath::Sadd(), Vmath::Smul(), Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vvtvp(), and Nektar::SolverUtils::EquationSystem::WriteFld().
Referenced by v_SetInitialConditions().
◆ Compute_demdt_cdot_ek()
|
protected |
Definition at line 1322 of file MMFSWE.cpp.
1326{
1327 int j, k;
1329
1330 Array<OneD, NekDouble> tmp(nq, 0.0);
1331 Array<OneD, NekDouble> tmpderiv(nq);
1332
1333 outarray = Array<OneD, NekDouble>(nq, 0.0);
1335 {
1337 {
1338 // Compute d e^m / d \xi_1 and d e^m / d \xi_2
1341
1342 Vmath::Vmul(nq, &physarray[j + 1][0], 1, &tmpderiv[0], 1,
1343 &tmpderiv[0], 1);
1344
1346 &outarray[0], 1, &outarray[0], 1);
1347 }
1348 }
1349}
References Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_shapedim, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Vcopy(), Vmath::Vmul(), and Vmath::Vvtvp().
Referenced by AddRotation().
◆ ComputeEnergy()
|
protected |
Definition at line 2151 of file MMFSWE.cpp.
2154{
2156
2157 Array<OneD, NekDouble> tmp(nq);
2158 Array<OneD, NekDouble> htmp(nq);
2159 Array<OneD, NekDouble> hstmp(nq);
2160
2161 Vmath::Vmul(nq, u, 1, u, 1, tmp, 1);
2162 Vmath::Vvtvp(nq, v, 1, v, 1, tmp, 1, tmp, 1);
2163 Vmath::Vmul(nq, eta, 1, tmp, 1, tmp, 1);
2164
2166 Vmath::Vmul(nq, htmp, 1, htmp, 1, htmp, 1);
2167
2169 Vmath::Vmul(nq, hstmp, 1, hstmp, 1, hstmp, 1);
2170
2171 Vmath::Vsub(nq, htmp, 1, hstmp, 1, htmp, 1);
2173
2174 Vmath::Vadd(nq, htmp, 1, tmp, 1, tmp, 1);
2175 Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2176
2178}
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.
Definition: Vmath.cpp:414
References m_depth, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, Vmath::Sadd(), Vmath::Smul(), Vmath::Vadd(), Vmath::Vmul(), Vmath::Vsub(), and Vmath::Vvtvp().
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ ComputeEnstrophy()
|
protected |
Definition at line 2180 of file MMFSWE.cpp.
2183{
2185
2186 Array<OneD, NekDouble> hstartmp(nq);
2187 Array<OneD, NekDouble> tmp(nq);
2188
2190
2191 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2192 Array<OneD, Array<OneD, NekDouble>> Curlu(m_spacedim);
2194 {
2195 uvec[i] = Array<OneD, NekDouble>(nq);
2196 Curlu[i] = Array<OneD, NekDouble>(nq);
2197 }
2198
2199 ComputeVorticity(u, v, tmp);
2200
2202 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
2203 Vmath::Vdiv(nq, tmp, 1, hstartmp, 1, tmp, 1);
2204 Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2205
2207}
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:280
References ComputeVorticity(), m_coriolis, m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, Vmath::Smul(), Vmath::Vadd(), Vmath::Vdiv(), and Vmath::Vmul().
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ Computehhuhvflux()
|
protected |
Definition at line 894 of file MMFSWE.cpp.
898{
900
901 NekDouble hC, huC, hvC, SL, SR, Sstar;
904
905 // Compute SL
906 if (hstar > hL)
907 SL = uL - cL * sqrt(0.5 * ((hstar * hstar + hstar * hL) / (hL * hL)));
908 else
909 SL = uL - cL;
910
911 // Compute SR
912 if (hstar > hR)
913 SR = uR + cR * sqrt(0.5 * ((hstar * hstar + hstar * hR) / (hR * hR)));
914 else
915 SR = uR + cR;
916
917 if (fabs(hR * (uR - SR) - hL * (uL - SL)) <= 1.0e-15)
918 Sstar = 0.0;
919 else
920 Sstar = (SL * hR * (uR - SR) - SR * hL * (uL - SL)) /
921 (hR * (uR - SR) - hL * (uL - SL));
922
923 if (SL >= 0)
924 {
925 hflux = hL * uL;
926 huflux = uL * uL * hL + 0.5 * g * hL * hL;
927 hvflux = hL * uL * vL;
928 }
929 else if (SR <= 0)
930 {
931 hflux = hR * uR;
932 huflux = uR * uR * hR + 0.5 * g * hR * hR;
933 hvflux = hR * uR * vR;
934 }
935 else
936 {
937 if ((SL < 0) && (Sstar >= 0))
938 {
939 hC = hL * ((SL - uL) / (SL - Sstar));
940 huC = hC * Sstar;
941 hvC = hC * vL;
942
943 hflux = hL * uL + SL * (hC - hL);
944 huflux = (uL * uL * hL + 0.5 * g * hL * hL) + SL * (huC - hL * uL);
945 hvflux = (uL * vL * hL) + SL * (hvC - hL * vL);
946 }
947 else
948 {
949 hC = hR * ((SR - uR) / (SR - Sstar));
950 huC = hC * Sstar;
951 hvC = hC * vR;
952
953 hflux = hR * uR + SR * (hC - hR);
954 huflux = (uR * uR * hR + 0.5 * g * hR * hR) + SR * (huC - hR * uR);
955 hvflux = (uR * vR * hR) + SR * (hvC - hR * vR);
956 }
957 }
958}
References m_g, and tinysimd::sqrt().
Referenced by RiemannSolverHLLC().
◆ ComputeMagAndDot()
|
protected |
Definition at line 1163 of file MMFSWE.cpp.
1168{
1169 NekDouble MF1x, MF1y, MF1z, MF2x, MF2y, MF2z;
1170 NekDouble MB1x, MB1y, MB1z, MB2x, MB2y, MB2z;
1171
1172 MF1x = m_MFtraceFwd[0][0][index];
1173 MF1y = m_MFtraceFwd[0][1][index];
1174 MF1z = m_MFtraceFwd[0][2][index];
1175
1176 MF2x = m_MFtraceFwd[1][0][index];
1177 MF2y = m_MFtraceFwd[1][1][index];
1178 MF2z = m_MFtraceFwd[1][2][index];
1179
1180 MB1x = m_MFtraceBwd[0][0][index];
1181 MB1y = m_MFtraceBwd[0][1][index];
1182 MB1z = m_MFtraceBwd[0][2][index];
1183
1184 MB2x = m_MFtraceBwd[1][0][index];
1185 MB2y = m_MFtraceBwd[1][1][index];
1186 MB2z = m_MFtraceBwd[1][2][index];
1187
1188 // MFtrace = MFtrace [ j*spacedim + k ], j = shape, k = sapce
1189 MageF1 = MF1x * MF1x + MF1y * MF1y + MF1z * MF1z;
1190 MageF2 = MF2x * MF2x + MF2y * MF2y + MF2z * MF2z;
1191 MageB1 = MB1x * MB1x + MB1y * MB1y + MB1z * MB1z;
1192 MageB2 = MB2x * MB2x + MB2y * MB2y + MB2z * MB2z;
1193
1194 eF1_cdot_eB1 = MF1x * MB1x + MF1y * MB1y + MF1z * MB1z;
1195 eF1_cdot_eB2 = MF1x * MB2x + MF1y * MB2y + MF1z * MB2z;
1196 eF2_cdot_eB1 = MF2x * MB1x + MF2y * MB1y + MF2z * MB1z;
1197 eF2_cdot_eB2 = MF2x * MB2x + MF2y * MB2y + MF2z * MB2z;
1198}
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
Definition: MMFSystem.h:198
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
Definition: MMFSystem.h:199
References Nektar::SolverUtils::MMFSystem::m_MFtraceBwd, and Nektar::SolverUtils::MMFSystem::m_MFtraceFwd.
Referenced by AverageFlux(), LaxFriedrichFlux(), and RusanovFlux().
◆ ComputeMass()
Definition at line 2141 of file MMFSWE.cpp.
2142{
2144
2145 Array<OneD, NekDouble> tmp(nq);
2147
2149}
References m_depth, Nektar::SolverUtils::EquationSystem::m_fields, and Vmath::Vadd().
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ ComputeNablaCdotVelocity()
Definition at line 2232 of file MMFSWE.cpp.
2233{
2235
2236 Array<OneD, NekDouble> velcoeff(nq, 0.0);
2237
2238 Array<OneD, NekDouble> Dtmp0(nq);
2239 Array<OneD, NekDouble> Dtmp1(nq);
2240 Array<OneD, NekDouble> Dtmp2(nq);
2241 Array<OneD, NekDouble> Drv(nq);
2242
2243 vellc = Array<OneD, NekDouble>(nq, 0.0);
2244
2245 // m_vellc = \nabla m_vel \cdot tan_i
2246 Array<OneD, NekDouble> tmp(nq);
2247 Array<OneD, NekDouble> vessel(nq);
2248
2250 {
2251 Vmath::Zero(nq, velcoeff, 1);
2253 {
2254 // a_j = tan_j cdot m_vel
2256 1, &velcoeff[0], 1, &velcoeff[0], 1);
2257 }
2258
2259 // d a_j / d x^k
2260 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
2261
2263 {
2264 // tan_j^k ( d a_j / d x^k )
2265 switch (k)
2266 {
2267 case (0):
2268 {
2270 1, &vellc[0], 1, &vellc[0], 1);
2271 }
2272 break;
2273
2274 case (1):
2275 {
2277 1, &vellc[0], 1, &vellc[0], 1);
2278 }
2279 break;
2280
2281 case (2):
2282 {
2284 1, &vellc[0], 1, &vellc[0], 1);
2285 }
2286 break;
2287 }
2288 }
2289 }
2290}
References Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_shapedim, Nektar::SolverUtils::EquationSystem::m_spacedim, m_velocity, Vmath::Vvtvp(), and Vmath::Zero().
◆ ComputeUnstableJetEta()
Definition at line 2737 of file MMFSWE.cpp.
2738{
2739 NekDouble uphi, f, dh;
2740
2741 uphi = ComputeUnstableJetuphi(theta);
2742 f = 2.0 * m_Omega * sin(theta);
2743
2744 dh = f * uphi + tan(theta) * uphi * uphi;
2745
2746 return dh;
2747}
NekDouble ComputeUnstableJetuphi(const NekDouble theta)
Definition: MMFSWE.cpp:2749
References ComputeUnstableJetuphi(), and m_Omega.
Referenced by UnstableJetFlow().
◆ ComputeUnstableJetuphi()
Definition at line 2749 of file MMFSWE.cpp.
2750{
2751 NekDouble uphi;
2752
2754 {
2757 }
2758
2759 else
2760 {
2761 uphi = 0.0;
2762 }
2763
2764 return uphi;
2765}
References m_en, m_theta0, m_theta1, and m_uthetamax.
Referenced by ComputeUnstableJetEta(), and UnstableJetFlow().
◆ ComputeVorticity()
|
protected |
Definition at line 2211 of file MMFSWE.cpp.
2214{
2216
2217 Array<OneD, NekDouble> tmp(nq);
2218
2219 Vorticity = Array<OneD, NekDouble>(nq, 0.0);
2220
2223 &Vorticity[0], 1);
2224
2226 Vmath::Neg(nq, tmp, 1);
2228
2229 Vmath::Vadd(nq, tmp, 1, Vorticity, 1, Vorticity, 1);
2230}
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Definition: MMFSystem.h:195
References Nektar::SolverUtils::MMFSystem::m_CurlMF, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Vmath::Neg(), Vmath::Vadd(), and Vmath::Vvtvp().
Referenced by Checkpoint_Output_Cartesian(), ComputeEnstrophy(), TestVorticityComputation(), and v_DoSolve().
◆ ConservativeToPrimitive() [1/2]
Definition at line 2912 of file MMFSWE.cpp.
2913{
2915
2916 // u = hu/h
2918 m_fields[1]->UpdatePhys(), 1);
2919
2920 // v = hv/ v
2922 m_fields[2]->UpdatePhys(), 1);
2923
2924 // \eta = h - d
2926 m_fields[0]->UpdatePhys(), 1);
2927}
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:765
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Vdiv(), and Vmath::Vsub().
Referenced by DoOdeProjection(), DoOdeRhs(), and v_DoSolve().
◆ ConservativeToPrimitive() [2/2]
|
protected |
Definition at line 2832 of file MMFSWE.cpp.
2835{
2837
2838 if (physin[0].get() == physout[0].get())
2839 {
2840 // copy indata and work with tmp array
2841 Array<OneD, Array<OneD, NekDouble>> tmp(3);
2842 for (int i = 0; i < 3; ++i)
2843 {
2844 // deep copy
2845 tmp[i] = Array<OneD, NekDouble>(nq);
2846 Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2847 }
2848
2849 // \eta = h - d
2851
2852 // u = hu/h
2853 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
2854
2855 // v = hv/h
2856 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
2857 }
2858 else
2859 {
2860 // \eta = h - d
2862
2863 // u = hu/h
2864 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
2865
2866 // v = hv/h
2867 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
2868 }
2869}
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vsub().
◆ create()
|
inlinestatic |
Creates an instance of this class.
Definition at line 67 of file MMFSWE.h.
70 {
72 MemoryManager<MMFSWE>::AllocateSharedPtr(pSession, pGraph);
73 p->InitObject();
75 }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:168
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition: EquationSystem.h:66
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.
◆ DoOdeProjection()
|
protected |
Compute the projection.
Compute the projection for the linear advection equation.
- Parameters
-
inarray Given fields. outarray Calculated solution. time Time.
Definition at line 1358 of file MMFSWE.cpp.
1361{
1363 {
1365 {
1366 ConservativeToPrimitive(inarray, outarray);
1367 SetBoundaryConditions(outarray, time);
1368 PrimitiveToConservative(outarray, outarray);
1369 }
1370 break;
1371 default:
1373 break;
1374 }
1375}
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &inarray, NekDouble time)
Definition: MMFSWE.cpp:1377
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:422
References ASSERTL0, ConservativeToPrimitive(), Nektar::MultiRegions::eDiscontinuous, Nektar::SolverUtils::EquationSystem::m_projectionType, PrimitiveToConservative(), and SetBoundaryConditions().
Referenced by v_InitObject().
◆ DoOdeRhs()
|
protected |
Compute the RHS.
Definition at line 426 of file MMFSWE.cpp.
429{
430 boost::ignore_unused(time);
431
432 int i;
433 int nvariables = inarray.size();
436
437 // inarray in physical space
438 Array<OneD, Array<OneD, NekDouble>> physarray(nvariables);
439 Array<OneD, Array<OneD, NekDouble>> modarray(nvariables);
440
441 for (i = 0; i < nvariables; ++i)
442 {
443 physarray[i] = Array<OneD, NekDouble>(nq);
444 modarray[i] = Array<OneD, NekDouble>(ncoeffs);
445 }
446
447 // (h, hu, hv) -> (\eta, u, v)
448 ConservativeToPrimitive(inarray, physarray);
449
450 // Weak Directional Derivative
451 WeakDGSWEDirDeriv(physarray, modarray);
452 AddDivForGradient(physarray, modarray);
453
454 for (i = 0; i < nvariables; ++i)
455 {
456 Vmath::Neg(ncoeffs, modarray[i], 1);
457 }
458
459 // coriolis forcing
461 {
462 AddCoriolis(physarray, modarray);
463 }
464
465 // Bottom Elevation Effect
466 // Add g H \nabla m_depth
467 AddElevationEffect(physarray, modarray);
468
469 // Add terms concerning the rotation of the moving frame
471 {
472 AddRotation(physarray, modarray);
473 }
474
475 for (i = 0; i < nvariables; ++i)
476 {
477 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
478 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
479 }
480}
void AddDivForGradient(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:552
void AddCoriolis(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1201
void WeakDGSWEDirDeriv(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField)
Definition: MMFSWE.cpp:482
void AddRotation(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1275
void AddElevationEffect(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1246
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:707
References AddCoriolis(), AddDivForGradient(), AddElevationEffect(), AddRotation(), ConservativeToPrimitive(), Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_AddCoriolis, m_AddRotation, Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Neg(), and WeakDGSWEDirDeriv().
Referenced by v_InitObject().
◆ EvaluateCoriolis()
Definition at line 1697 of file MMFSWE.cpp.
1698{
1700 {
1702 {
1704 }
1705 break;
1706
1708 {
1710 }
1711 break;
1712
1717 {
1719 }
1720 break;
1721
1722 default:
1723 break;
1724 }
1725}
void EvaluateStandardCoriolis(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1759
void EvaluateCoriolisForZonalFlow(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1727
SOLVER_UTILS_EXPORT SessionFunctionSharedPtr GetFunction(std::string name, const MultiRegions::ExpListSharedPtr &field=MultiRegions::NullExpListSharedPtr, bool cache=false)
Get a SessionFunction by name.
Definition: EquationSystem.cpp:737
References Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, EvaluateCoriolisForZonalFlow(), EvaluateStandardCoriolis(), Nektar::SolverUtils::EquationSystem::GetFunction(), m_coriolis, and m_TestType.
Referenced by v_DoInitialise().
◆ EvaluateCoriolisForZonalFlow()
Definition at line 1727 of file MMFSWE.cpp.
1728{
1730 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1731
1732 Array<OneD, NekDouble> x(nq);
1733 Array<OneD, NekDouble> y(nq);
1734 Array<OneD, NekDouble> z(nq);
1735
1737
1738 NekDouble x0j, x1j, x2j;
1739 NekDouble tmp;
1740
1741 outarray = Array<OneD, NekDouble>(nq, 0.0);
1742 for (int j = 0; j < nq; ++j)
1743 {
1744 x0j = x[j];
1745 x1j = y[j];
1746 x2j = z[j];
1747
1748 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1749 cos_theta);
1750
1751 // H = 2 \Omega *(- \cos \phi \cos \theta \sin \alpha + \sin \theta \cos
1752 // \alpha )
1753 tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
1754 sin_theta * cos(m_alpha);
1755 outarray[j] = 2.0 * m_Omega * tmp;
1756 }
1757}
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)
Definition: MMFSystem.cpp:778
std::vector< double > z(NPUPPER)
References Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_Omega, and Nektar::UnitTests::z().
Referenced by EvaluateCoriolis().
◆ EvaluateStandardCoriolis()
Definition at line 1759 of file MMFSWE.cpp.
1760{
1762
1763 NekDouble x0j, x1j, x2j;
1764 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1765
1766 Array<OneD, NekDouble> x(nq);
1767 Array<OneD, NekDouble> y(nq);
1768 Array<OneD, NekDouble> z(nq);
1769
1771
1772 outarray = Array<OneD, NekDouble>(nq, 0.0);
1773 for (int j = 0; j < nq; ++j)
1774 {
1775 x0j = x[j];
1776 x1j = y[j];
1777 x2j = z[j];
1778
1779 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1780 cos_theta);
1781
1782 outarray[j] = 2.0 * m_Omega * sin_theta;
1783 }
1784}
References Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_Omega, and Nektar::UnitTests::z().
Referenced by EvaluateCoriolis().
◆ EvaluateWaterDepth()
Definition at line 1557 of file MMFSWE.cpp.
1558{
1560
1561 m_depth = Array<OneD, NekDouble>(nq, 0.0);
1562
1564 {
1566 {
1568 }
1569 break;
1570
1572 {
1573 // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1574 Array<OneD, NekDouble> x(nq);
1575 Array<OneD, NekDouble> y(nq);
1576 Array<OneD, NekDouble> z(nq);
1577
1579
1580 NekDouble x0j, x1j, x2j;
1581 NekDouble sin_varphi, cos_varphi, sin_theta, cos_theta;
1582
1583 for (int j = 0; j < nq; ++j)
1584 {
1585 x0j = x[j];
1586 x1j = y[j];
1587 x2j = z[j];
1588
1589 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1590 sin_theta, cos_theta);
1591 m_depth[j] =
1594 }
1595 }
1596 break;
1597
1599 {
1600 // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1601 Array<OneD, NekDouble> x(nq);
1602 Array<OneD, NekDouble> y(nq);
1603 Array<OneD, NekDouble> z(nq);
1604
1606
1607 int indx = 0;
1608 NekDouble x0j, x1j, x2j, dist2;
1609 NekDouble phi, theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
1610 NekDouble hRad, phic, thetac;
1611 NekDouble Tol = 0.000001;
1612
1613 hRad = m_pi / 9.0;
1614 phic = -m_pi / 2.0;
1615 thetac = m_pi / 6.0;
1616
1617 for (int j = 0; j < nq; ++j)
1618 {
1619 x0j = x[j];
1620 x1j = y[j];
1621 x2j = z[j];
1622
1623 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1624 sin_theta, cos_theta);
1625
1627 std::abs(sin(thetac) - sin_theta)) < Tol)
1628 {
1629 std::cout << "A point " << j
1630 << " is coincient with the singularity "
1631 << std::endl;
1632 indx = 1;
1633 }
1634
1635 phi = atan2(sin_varphi, cos_varphi);
1636 theta = atan2(sin_theta, cos_theta);
1637
1638 // Compute r
1639 dist2 = (phi - phic) * (phi - phic) +
1640 (theta - thetac) * (theta - thetac);
1641
1642 if (dist2 > hRad * hRad)
1643 {
1644 dist2 = hRad * hRad;
1645 }
1646
1648 }
1649
1650 if (!indx)
1651 {
1652 std::cout << "No point is coincident with the singularity point"
1653 << std::endl;
1654 }
1655 }
1656 break;
1657
1659 {
1660 for (int j = 0; j < nq; ++j)
1661 {
1663 }
1664 }
1665 break;
1666
1669 {
1671 }
1672 break;
1673
1674 default:
1675 {
1677 }
1678 break;
1679 }
1680
1681 // Comptue \nabla m_depth \cdot e^i
1683
1685 {
1686 m_Derivdepth[j] = Array<OneD, NekDouble>(nq);
1688 m_Derivdepth[j]);
1689 }
1690
1691 std::cout << "Water Depth (m_depth) was generated with mag = "
1695}
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2296
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:43
T Vamax(int n, const T *x, const int incx)
Return the maximum absolute element in x called vamax to avoid conflict with max.
Definition: Vmath.cpp:994
References tinysimd::abs(), Nektar::SolverUtils::MMFSystem::AvgAbsInt(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, Vmath::Fill(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, m_Derivdepth, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, m_hs0, m_k2, Nektar::SolverUtils::MMFSystem::m_movingframes, m_Omega, Nektar::SolverUtils::MMFSystem::m_pi, Nektar::SolverUtils::MMFSystem::m_shapedim, m_TestType, tinysimd::sqrt(), Vmath::Vamax(), Nektar::UnitTests::z(), and Vmath::Zero().
Referenced by v_DoInitialise().
◆ GetNormalVelocity()
Get the normal velocity.
◆ GetSWEFluxVector()
|
protected |
Definition at line 582 of file MMFSWE.cpp.
585{
587
588 switch (i)
589 {
590 // flux function for the h equation = [(\eta + d ) u, (\eta + d) v ]
591 case 0:
592 {
593 // h in flux 1
595
596 // hu in flux 0
597 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
598
599 // hv in flux 1
600 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
601 }
602 break;
603
604 // flux function for the hu equation = [Hu^2 + 0.5 g H^2, Huv]
605 case 1:
606 {
607 Array<OneD, NekDouble> tmp(nq);
608
609 // h in tmp
611
612 // hu in flux 1
613 Vmath::Vmul(nq, tmp, 1, physfield[1], 1, flux[1], 1);
614
615 // huu in flux 0
616 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
617
618 // hh in tmp
619 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
620
621 // huu + 0.5 g hh in flux 0
622 // Daxpy overwrites flux[0] on exit
624
625 // huv in flux 1
626 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
627 }
628 break;
629
630 // flux function for the hv equation = [Huv, Hv^2]
631 case 2:
632 {
633 Array<OneD, NekDouble> tmp(nq);
634
635 // h in tmp
637
638 // hv in flux 0
639 Vmath::Vmul(nq, tmp, 1, physfield[2], 1, flux[0], 1);
640
641 // hvv in flux 1
642 Vmath::Vmul(nq, flux[0], 1, physfield[2], 1, flux[1], 1);
643
644 // huv in flux 0
645 Vmath::Vmul(nq, flux[0], 1, physfield[1], 1, flux[0], 1);
646
647 // hh in tmp
648 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
649
650 // hvv + 0.5 g hh in flux 1
652 }
653 break;
654
655 // flux function 0.5 g h * h
656 case 3:
657 {
658 Array<OneD, NekDouble> h(nq);
659 flux[0] = Array<OneD, NekDouble>(nq, 0.0);
660
661 // h in tmp
663
664 // hh in tmp
665 Vmath::Vmul(nq, h, 1, h, 1, h, 1);
666
667 // 0.5 g hh in flux 0
669 }
670 break;
671
672 default:
673 break;
674 }
675}
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:137
References Blas::Daxpy(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, m_g, Vmath::Vadd(), and Vmath::Vmul().
Referenced by AddDivForGradient(), and WeakDGSWEDirDeriv().
◆ IsolatedMountainFlow()
|
protected |
Definition at line 2393 of file MMFSWE.cpp.
2395{
2396 boost::ignore_unused(time);
2397
2399
2400 NekDouble uhat, vhat;
2401 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2402 NekDouble x0j, x1j, x2j;
2403
2404 Array<OneD, NekDouble> eta(nq, 0.0);
2405 Array<OneD, NekDouble> u(nq, 0.0);
2406 Array<OneD, NekDouble> v(nq, 0.0);
2407
2408 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2410 {
2411 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2412 }
2413
2414 Array<OneD, NekDouble> x(nq);
2415 Array<OneD, NekDouble> y(nq);
2416 Array<OneD, NekDouble> z(nq);
2417
2419
2420 for (int j = 0; j < nq; ++j)
2421 {
2422 x0j = x[j];
2423 x1j = y[j];
2424 x2j = z[j];
2425
2426 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2427 cos_theta);
2428
2429 // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2431 sin_theta * sin_theta;
2432
2433 uhat = m_u0 * cos_theta;
2434 vhat = 0.0;
2435
2436 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2437 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2438 uvec[2][j] = vhat * cos_theta;
2439 }
2440
2441 // Projection of u onto the tangent plane with conserving the mag. of the
2442 // velocity.
2443 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2444
2446 {
2447 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2448 }
2449
2450 // u is projected on the tangent plane with preserving its length
2451 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2452
2453 // Change it to the coordinate of moving frames
2454 // CartesianToMovingframes(0,uvecproj,u);
2455 // CartesianToMovingframes(1,uvecproj,v);
2456
2457 u = CartesianToMovingframes(uvec, 0);
2458 v = CartesianToMovingframes(uvec, 1);
2459
2460 switch (field)
2461 {
2462 case (0):
2463 {
2464 outfield = eta;
2465 }
2466 break;
2467
2468 case (1):
2469 {
2470 outfield = u;
2471 }
2472 break;
2473
2474 case (2):
2475 {
2476 outfield = v;
2477 }
2478 break;
2479 }
2480}
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble > > &uvec, unsigned int field)
Definition: MMFSystem.cpp:757
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_Omega, Nektar::SolverUtils::EquationSystem::m_spacedim, m_u0, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ LaxFriedrichFlux()
|
protected |
Definition at line 1014 of file MMFSWE.cpp.
1018{
1019 int nvariables = 3;
1020 NekDouble MageF1, MageF2, MageB1, MageB2;
1021 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1022 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1023
1025 NekDouble uRF, vRF, uLB, vLB;
1026 NekDouble velL, velR, lambdaF, lambdaB;
1027
1028 Array<OneD, NekDouble> EigF(nvariables);
1029 Array<OneD, NekDouble> EigB(nvariables);
1030
1031 // Compute Magnitude and Dot product of moving frames for the index
1032 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1033 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1034
1035 // Get the velocity in the normal to the edge
1038
1039 EigF[0] = velL - sqrt(g * hL);
1040 EigF[1] = velL;
1041 EigF[2] = velL + sqrt(g * hL);
1042
1043 EigB[0] = velR - sqrt(g * hR);
1044 EigB[1] = velR;
1045 EigB[2] = velR + sqrt(g * hR);
1046
1047 lambdaF = Vmath::Vamax(nvariables, EigF, 1);
1048 lambdaB = Vmath::Vamax(nvariables, EigB, 1);
1049
1050 // uRF = uR component in moving frames e^{Fwd}
1051 // vRF = vR component in moving frames e^{Fwd}
1052 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1053 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1054
1055 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1056 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1057 numfluxF[2] = 0.5 * lambdaF * (hL - hR);
1058
1059 numfluxF[3] = 0.5 * (hL * uL * uL * MageF1 + hR * uRF * uRF * MageB1 +
1060 0.5 * g * (hL * hL + hR * hR));
1061 numfluxF[4] = 0.5 * (hL * uL * vL * MageF1 + hR * uRF * vRF * MageB1);
1062 numfluxF[5] = 0.5 * lambdaF * (uL * hL - uRF * hR);
1063
1064 numfluxF[6] = 0.5 * (hL * uL * vL * MageF2 + hR * uRF * vRF * MageB2);
1065 numfluxF[7] = 0.5 * (hL * vL * vL * MageF2 + hR * vRF * vRF * MageB2 +
1066 0.5 * g * (hL * hL + hR * hR));
1067 numfluxF[8] = 0.5 * lambdaF * (vL * hL - vRF * hR);
1068
1069 // uLB = uL component in moving frames e^{Bwd}
1070 // vLB = vL component in moving frames e^{Bwd}
1071 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1072 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1073
1074 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1075 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1076 numfluxB[2] = 0.5 * lambdaB * (hL - hR);
1077
1078 numfluxB[3] = 0.5 * (hR * uR * uR * MageB1 + hR * uLB * uLB * MageF1 +
1079 0.5 * g * (hR * hR + hL * hL));
1080 numfluxB[4] = 0.5 * (hR * uR * vR * MageB1 + hR * uLB * vLB * MageF1);
1081 numfluxB[5] = 0.5 * lambdaB * (uLB * hL - uR * hR);
1082
1083 numfluxB[6] = 0.5 * (hR * uR * vR * MageB2 + hR * uLB * vLB * MageF2);
1084 numfluxB[7] = 0.5 * (hR * vR * vR * MageB2 + hR * vLB * vLB * MageF2 +
1085 0.5 * g * (hR * hR + hL * hL));
1086 numfluxB[8] = 0.5 * lambdaB * (vLB * hL - vR * hR);
1087}
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:189
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:188
References ComputeMagAndDot(), m_g, Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd, Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd, tinysimd::sqrt(), and Vmath::Vamax().
Referenced by NumericalSWEFlux().
◆ NumericalSWEFlux()
|
protected |
Definition at line 677 of file MMFSWE.cpp.
680{
681
682 int i, k;
684 int nvariables = 3; // only the dependent variables
685
686 // get temporary arrays
687 Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
688 Array<OneD, Array<OneD, NekDouble>> Bwd(nvariables);
689 Array<OneD, NekDouble> DepthFwd(nTraceNumPoints);
690 Array<OneD, NekDouble> DepthBwd(nTraceNumPoints);
691
692 for (i = 0; i < nvariables; ++i)
693 {
694 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
695 Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
696 }
697
698 // get the physical values at the trace
699 for (i = 0; i < nvariables; ++i)
700 {
701 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
702
703 // Copy Fwd to Bwd at boundaries
705 }
708
709 // note that we are using the same depth - i.e. the depth is assumed
710 // continuous...
712 {
714 {
715 // rotate the values to the normal direction
716 NekDouble tmpX, tmpY;
717
718 // Fwd[1] = hu^+ = ( h^1^+ (e^1 \cdot n) + h^2^+ ( e^2 \cdot n) )
719 // Fwd[2] = hv^+ = ( h^1^+ (e^1 \cdot n^{\perp}) + h^2^+ ( e^2 \cdot
720 // n^{\perp}) )
721 // Bwd[1] = hu^- = ( h^1^- (e^1 \cdot n) + h^2^- ( e^2 \cdot n) )
722 // Bwd[2] = hv^- = ( h^1^- (e^1 \cdot n^{\perp}) + h^2^- ( e^2 \cdot
723 // n^{\perp}) )
724
725 for (k = 0; k < nTraceNumPoints; ++k)
726 {
727 tmpX = Fwd[1][k] * m_ncdotMFFwd[0][k] +
728 Fwd[2][k] * m_ncdotMFFwd[1][k];
729 tmpY = Fwd[1][k] * m_nperpcdotMFFwd[0][k] +
730 Fwd[2][k] * m_nperpcdotMFFwd[1][k];
731 Fwd[1][k] = tmpX;
732 Fwd[2][k] = tmpY;
733
734 tmpX = Bwd[1][k] * m_ncdotMFBwd[0][k] +
735 Bwd[2][k] * m_ncdotMFBwd[1][k];
736 tmpY = Bwd[1][k] * m_nperpcdotMFBwd[0][k] +
737 Bwd[2][k] * m_nperpcdotMFBwd[1][k];
738 Bwd[1][k] = tmpX;
739 Bwd[2][k] = tmpY;
740 }
741
742 // Solve the Riemann problem
743 NekDouble denomFwd, denomBwd;
744 Array<OneD, NekDouble> numfluxF(nvariables);
745 Array<OneD, NekDouble> numfluxB(nvariables);
746
747 NekDouble eF1n, eF2n, eB1n, eB2n;
748 NekDouble eF1t, eF2t, eB1t, eB2t;
749 for (k = 0; k < nTraceNumPoints; ++k)
750 {
751 RiemannSolverHLLC(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
752 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
753 Bwd[2][k], numfluxF, numfluxB);
754
755 // uflux = 1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
756 // \cdot n)(e^1 \cdot n^{\perp} ) ] ( e^2 \cdot n^{\perp} huflux
757 // - e^2 \cdot n hvflux
758 // vflux = -1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
759 // \cdot n)(e^1 \cdot n^{\perp} ) ] ( -e^1 \cdot n^{\perp}
760 // huflux + e^1 \cdot n hvflux
761 eF1n = m_ncdotMFFwd[0][k];
762 eF2n = m_ncdotMFFwd[1][k];
763 eB1n = m_ncdotMFBwd[0][k];
764 eB2n = m_ncdotMFBwd[1][k];
765
766 eF1t = m_nperpcdotMFFwd[0][k];
767 eF2t = m_nperpcdotMFFwd[1][k];
768 eB1t = m_nperpcdotMFBwd[0][k];
769 eB2t = m_nperpcdotMFBwd[1][k];
770
771 denomFwd = eF1n * eF2t - eF2n * eF1t;
772 denomBwd = eB1n * eB2t - eB2n * eB1t;
773
774 numfluxFwd[0][k] = numfluxF[0];
775 numfluxFwd[1][k] = (1.0 / denomFwd) *
776 (eF2t * numfluxF[1] - eF2n * numfluxF[2]);
777 numfluxFwd[2][k] =
778 (1.0 / denomFwd) *
779 (-1.0 * eF1t * numfluxF[1] + eF1n * numfluxF[2]);
780
781 numfluxBwd[0][k] = 1.0 * numfluxB[0];
782 numfluxBwd[1][k] = (1.0 / denomBwd) *
783 (eB2t * numfluxB[1] - eB2n * numfluxB[2]);
784 numfluxBwd[2][k] =
785 (1.0 / denomBwd) *
786 (-1.0 * eB1t * numfluxB[1] + eB1n * numfluxB[2]);
787 }
788 }
789 break;
790
794 {
795 Array<OneD, NekDouble> numfluxF(nvariables * (m_shapedim + 1));
796 Array<OneD, NekDouble> numfluxB(nvariables * (m_shapedim + 1));
797
798 for (k = 0; k < nTraceNumPoints; ++k)
799 {
801 {
802 AverageFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
803 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
804 Bwd[2][k], numfluxF, numfluxB);
805 }
806
808 {
809 LaxFriedrichFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
810 Fwd[2][k], Bwd[0][k] + DepthFwd[k],
811 Bwd[1][k], Bwd[2][k], numfluxF, numfluxB);
812 }
813
815 {
816 RusanovFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
817 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
818 Bwd[2][k], numfluxF, numfluxB);
819 }
820
821 int indx;
822 NekDouble tmpF0, tmpF1, tmpB0, tmpB1;
823 for (i = 0; i < nvariables; ++i)
824 {
825 indx = i * (m_shapedim + 1);
826
827 tmpF0 = numfluxF[indx] * m_ncdotMFFwd[0][k];
828 tmpF1 = numfluxF[indx + 1] * m_ncdotMFFwd[1][k];
829
830 tmpB0 = numfluxB[indx] * m_ncdotMFBwd[0][k];
831 tmpB1 = numfluxB[indx + 1] * m_ncdotMFBwd[1][k];
832
833 numfluxFwd[i][k] = tmpF0 + tmpF1 + numfluxF[indx + 2];
834 numfluxBwd[i][k] = tmpB0 + tmpB1 + numfluxB[indx + 2];
835 }
836 }
837 }
838 break;
839
840 default:
841 {
843 }
844 break;
845 }
846}
void AverageFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:960
void RiemannSolverHLLC(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:848
void LaxFriedrichFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:1014
void RusanovFlux(const int index, NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, Array< OneD, NekDouble > &numfluxF, Array< OneD, NekDouble > &numfluxB)
Definition: MMFSWE.cpp:1089
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:740
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
Definition: MMFSystem.h:191
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFBwd
Definition: MMFSystem.h:192
SOLVER_UTILS_EXPORT void CopyBoundaryTrace(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const BoundaryCopyType BDCopyType, const int var=0, const std::string btype="NoUserDefined")
Definition: MMFSystem.cpp:821
References ASSERTL0, AverageFlux(), Nektar::SolverUtils::MMFSystem::CopyBoundaryTrace(), Nektar::SolverUtils::eAverage, Nektar::SolverUtils::eFwdEQBwd, Nektar::SolverUtils::eHLLC, Nektar::SolverUtils::eLaxFriedrich, Nektar::SolverUtils::eRusanov, Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), LaxFriedrichFlux(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd, Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd, Nektar::SolverUtils::MMFSystem::m_nperpcdotMFBwd, Nektar::SolverUtils::MMFSystem::m_nperpcdotMFFwd, Nektar::SolverUtils::MMFSystem::m_shapedim, Nektar::SolverUtils::MMFSystem::m_upwindType, RiemannSolverHLLC(), and RusanovFlux().
Referenced by WeakDGSWEDirDeriv().
◆ PrimitiveToConservative() [1/2]
Definition at line 2929 of file MMFSWE.cpp.
2930{
2932
2933 // h = \eta + d
2935 m_fields[0]->UpdatePhys(), 1);
2936
2937 // hu = h * u
2939 m_fields[1]->UpdatePhys(), 1);
2940
2941 // hv = h * v
2943 m_fields[2]->UpdatePhys(), 1);
2944}
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Vadd(), and Vmath::Vmul().
Referenced by DoOdeProjection(), v_DoInitialise(), and v_DoSolve().
◆ PrimitiveToConservative() [2/2]
|
protected |
Definition at line 2871 of file MMFSWE.cpp.
2874{
2875
2877
2878 if (physin[0].get() == physout[0].get())
2879 {
2880 // copy indata and work with tmp array
2881 Array<OneD, Array<OneD, NekDouble>> tmp(3);
2882 for (int i = 0; i < 3; ++i)
2883 {
2884 // deep copy
2885 tmp[i] = Array<OneD, NekDouble>(nq);
2886 Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2887 }
2888
2889 // h = \eta + d
2891
2892 // hu = h * u
2893 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
2894
2895 // hv = h * v
2896 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
2897 }
2898 else
2899 {
2900 // h = \eta + d
2902
2903 // hu = h * u
2904 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
2905
2906 // hv = h * v
2907 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
2908 }
2909}
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vmul().
◆ RiemannSolverHLLC()
|
protected |
Definition at line 848 of file MMFSWE.cpp.
852{
853 boost::ignore_unused(index);
854
856
859
860 NekDouble uRF, vRF, uLB, vLB;
861 NekDouble hstarF, hstarB;
862
863 // Temporary assignments
864 uRF = uR;
865 vRF = vR;
866 uLB = uL;
867 vLB = vL;
868
869 // the two-rarefaction wave assumption
870 hstarF = 0.5 * (cL + cR) + 0.25 * (uL - uRF);
871 hstarF *= hstarF;
872 hstarF *= (1.0 / g);
873
874 hstarB = 0.5 * (cL + cR) + 0.25 * (uLB - uR);
875 hstarB *= hstarB;
876 hstarB *= (1.0 / g);
877
878 NekDouble hfluxF, hufluxF, hvfluxF;
879 NekDouble hfluxB, hufluxB, hvfluxB;
880 Computehhuhvflux(hL, uL, vL, hR, uRF, vRF, hstarF, hfluxF, hufluxF,
881 hvfluxF);
882 Computehhuhvflux(hL, uLB, vLB, hR, uR, vR, hstarB, hfluxB, hufluxB,
883 hvfluxB);
884
885 numfluxF[0] = hfluxF;
886 numfluxF[1] = hufluxF;
887 numfluxF[2] = hvfluxF;
888
889 numfluxB[0] = hfluxB;
890 numfluxB[1] = hufluxB;
891 numfluxB[2] = hvfluxB;
892}
void Computehhuhvflux(NekDouble hL, NekDouble uL, NekDouble vL, NekDouble hR, NekDouble uR, NekDouble vR, NekDouble hstar, NekDouble &hflux, NekDouble &huflux, NekDouble &hvflux)
Definition: MMFSWE.cpp:894
References Computehhuhvflux(), m_g, and tinysimd::sqrt().
Referenced by NumericalSWEFlux().
◆ RossbyWave()
|
protected |
Definition at line 2592 of file MMFSWE.cpp.
2593{
2595 NekDouble uhat, vhat;
2596 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2597 NekDouble x0j, x1j, x2j;
2598 NekDouble Ath, Bth, Cth, tmp;
2599
2600 Array<OneD, NekDouble> eta(nq, 0.0);
2601 Array<OneD, NekDouble> u(nq, 0.0);
2602 Array<OneD, NekDouble> v(nq, 0.0);
2603
2604 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2606 {
2607 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2608 }
2609
2610 Array<OneD, NekDouble> x(nq);
2611 Array<OneD, NekDouble> y(nq);
2612 Array<OneD, NekDouble> z(nq);
2613
2615
2616 NekDouble R = 4.0;
2617 NekDouble cos2theta, cosRtheta, cos6theta, cos2Rtheta, cosRm1theta;
2618 NekDouble cos2phi, cos4phi, sin4phi, cos8phi;
2619
2620 // disturbancees of Rossby-Haurwitz Wave
2621 NekDouble x0d, y0d, z0d, phi0, theta0;
2622 // NekDouble rad_earth = 6.37122 * 1000000;
2623
2624 phi0 = 40.0 * m_pi / 180.0;
2625 theta0 = 50.0 * m_pi / 180.0;
2626
2627 x0d = cos(phi0) * cos(theta0);
2628 y0d = sin(phi0) * cos(theta0);
2629 z0d = sin(theta0);
2630
2631 for (int j = 0; j < nq; ++j)
2632 {
2633 x0j = x[j];
2634 x1j = y[j];
2635 x2j = z[j];
2636
2637 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2638 cos_theta);
2639
2640 // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2641 // \sin \alpha + \sin \theta \cos \alpha )^2
2642
2643 // tmp = cos^{2R} \theta, R = 4;
2644 cos2theta = cos_theta * cos_theta;
2645 cosRm1theta = cos_theta * cos2theta;
2646 cosRtheta = cos2theta * cos2theta;
2647 cos6theta = cos2theta * cosRtheta;
2648 cos2Rtheta = cosRtheta * cosRtheta;
2649
2651 Ath = tmp * cos2theta +
2653 ((R + 1.0) * cosRtheta + (2 * R * R - R - 2.0) * cos2theta -
2654 2 * R * R);
2655
2657 Bth = tmp * cosRtheta *
2658 ((R * R + 2 * R + 2) - (R + 1.0) * (R + 1.0) * cos2theta);
2659
2660 Cth =
2662
2663 // cos (2 \phi) = 2 * cos \phi * cos \phi - 1.0
2664 cos2phi = 2.0 * cos_varphi * cos_varphi - 1.0;
2665 cos4phi = 2.0 * cos2phi * cos2phi - 1.0;
2666 cos8phi = 2.0 * cos4phi * cos4phi - 1.0;
2667
2668 // sin (2 \phi) = 2 * cos \phi * sin \phi
2669 sin4phi = 4.0 * sin_varphi * cos_varphi * cos2phi;
2670
2672
2673 // disturbances is added
2675 {
2676 eta[j] = eta[j] *
2677 (1.0 + (1.0 / 40.0) * (x0j * x0d + x1j * y0d + x2j * z0d));
2678 }
2679
2680 // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2681 // || e^2 ||^2
2682 uhat = m_angfreq * cos_theta +
2683 m_K * cosRm1theta * (R * sin_theta * sin_theta - cos2theta) *
2684 cos4phi;
2685 vhat = -1.0 * m_K * R * cosRm1theta * sin_theta * sin4phi;
2686
2687 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2688 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2689 uvec[2][j] = vhat * cos_theta;
2690 }
2691
2692 // NekDouble etamin, etaaver;
2693 // etamin = Vmath::Vmin(nq, eta, 1)*rad_earth;
2694 // etaaver = Vmath::Vsum(nq, eta, 1)/nq*rad_earth;
2695
2696 // Projection of u onto the tangent plane with conserving the mag. of the
2697 // velocity.
2698 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2699
2701 {
2702 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2703 }
2704
2705 // u is projected on the tangent plane with preserving its length
2706 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2707
2708 // Change it to the coordinate of moving frames
2709 // CartesianToMovingframes(0,uvecproj,u);
2710 // CartesianToMovingframes(1,uvecproj,v);
2711
2712 u = CartesianToMovingframes(uvec, 0);
2713 v = CartesianToMovingframes(uvec, 1);
2714
2715 switch (field)
2716 {
2717 case (0):
2718 {
2719 outfield = eta;
2720 }
2721 break;
2722
2723 case (1):
2724 {
2725 outfield = u;
2726 }
2727 break;
2728
2729 case (2):
2730 {
2731 outfield = v;
2732 }
2733 break;
2734 }
2735}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_angfreq, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, m_K, m_Omega, Nektar::SolverUtils::MMFSystem::m_pi, m_RossbyDisturbance, Nektar::SolverUtils::EquationSystem::m_spacedim, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ RusanovFlux()
|
protected |
Definition at line 1089 of file MMFSWE.cpp.
1093{
1094 int nvariables = 3;
1095 NekDouble MageF1, MageF2, MageB1, MageB2;
1096 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1097 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1098
1100 NekDouble uRF, vRF, uLB, vLB;
1101 NekDouble velL, velR;
1102
1103 Array<OneD, NekDouble> EigF(nvariables);
1104 Array<OneD, NekDouble> EigB(nvariables);
1105
1106 // Compute Magnitude and Dot product of moving frames for the index
1107 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1108 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1109
1110 // Get the velocity in the normal to the edge
1113
1114 NekDouble SL, SR;
1115 SL = fabs(velL) + sqrt(g * hL);
1116 SR = fabs(velR) + sqrt(g * hR);
1117
1118 NekDouble S;
1119 if (SL > SR)
1120 S = SL;
1121 else
1122 S = SR;
1123
1124 // uRF = uR component in moving frames e^{Fwd}
1125 // vRF = vR component in moving frames e^{Fwd}
1126 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1127 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1128
1129 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1130 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1131 numfluxF[2] = 0.5 * S * (hL - hR);
1132
1133 numfluxF[3] =
1134 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
1135 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1136 numfluxF[5] = 0.5 * S * (uL * hL - uRF * hR);
1137
1138 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1139 numfluxF[7] =
1140 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1141 numfluxF[8] = 0.5 * S * (vL * hL - vRF * hR);
1142
1143 // uLB = uL component in moving frames e^{Bwd}
1144 // vLB = vL component in moving frames e^{Bwd}
1145 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1146 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1147
1148 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1149 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1150 numfluxB[2] = 0.5 * S * (hL - hR);
1151
1152 numfluxB[3] =
1153 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1154 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1155 numfluxB[5] = 0.5 * S * (uLB * hL - uR * hR);
1156
1157 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1158 numfluxB[7] =
1159 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1160 numfluxB[8] = 0.5 * S * (vLB * hL - vR * hR);
1161}
References ComputeMagAndDot(), m_g, Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd, Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd, and tinysimd::sqrt().
Referenced by NumericalSWEFlux().
◆ SetBoundaryConditions()
|
protected |
Definition at line 1377 of file MMFSWE.cpp.
1379{
1380
1382 int cnt = 0;
1383
1384 // loop over Boundary Regions
1386 {
1387
1388 // Zonal Boundary Condition
1390 {
1392 {
1394 }
1396 {
1397 // ZonalBoundary2D(n,cnt,inarray);
1398 }
1399 }
1400
1401 // Wall Boundary Condition
1403 {
1405 {
1407 }
1409 {
1410 WallBoundary2D(n, cnt, inarray);
1411 }
1412 }
1413
1414 // Time Dependent Boundary Condition (specified in meshfile)
1416 "eTimeDependent")
1417 {
1418 for (int i = 0; i < nvariables; ++i)
1419 {
1420 m_fields[i]->EvaluateBoundaryConditions(time);
1421 }
1422 }
1423 cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
1424 }
1425}
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Definition: MMFSWE.cpp:1427
References ASSERTL0, Nektar::SolverUtils::EquationSystem::m_expdim, Nektar::SolverUtils::EquationSystem::m_fields, and WallBoundary2D().
Referenced by DoOdeProjection().
◆ SteadyZonalFlow()
|
protected |
Definition at line 2050 of file MMFSWE.cpp.
2052{
2054 NekDouble uhat, vhat;
2055 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2056 NekDouble x0j, x1j, x2j, tmp;
2057
2058 Array<OneD, NekDouble> eta(nq, 0.0);
2059 Array<OneD, NekDouble> u(nq, 0.0);
2060 Array<OneD, NekDouble> v(nq, 0.0);
2061
2062 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2064 {
2065 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2066 }
2067
2068 Array<OneD, NekDouble> x(nq);
2069 Array<OneD, NekDouble> y(nq);
2070 Array<OneD, NekDouble> z(nq);
2071
2073
2074 for (int j = 0; j < nq; ++j)
2075 {
2076 x0j = x[j];
2077 x1j = y[j];
2078 x2j = z[j];
2079
2080 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2081 cos_theta);
2082
2083 // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2084 // \sin \alpha + \sin \theta \cos \alpha )^2
2085 tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
2086 sin_theta * cos(m_alpha);
2088
2089 // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2090 // || e^2 ||^2
2092 sin_theta * cos_varphi * sin(m_alpha));
2094
2095 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2096 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2097 uvec[2][j] = vhat * cos_theta;
2098 }
2099
2100 // Projection of u onto the tangent plane with conserving the mag. of the
2101 // velocity.
2102 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2103
2105 {
2106 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2107 }
2108
2109 // u is projected on the tangent plane with preserving its length
2110 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2111
2112 // Change it to the coordinate of moving frames
2113 // CartesianToMovingframes(0,uvecproj,u);
2114 // CartesianToMovingframes(1,uvecproj,v);
2115
2116 u = CartesianToMovingframes(uvec, 0);
2117 v = CartesianToMovingframes(uvec, 1);
2118
2119 switch (field)
2120 {
2121 case (0):
2122 {
2123 outfield = eta;
2124 }
2125 break;
2126
2127 case (1):
2128 {
2129 outfield = u;
2130 }
2131 break;
2132
2133 case (2):
2134 {
2135 outfield = v;
2136 }
2137 break;
2138 }
2139}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_H0, m_Hvar, Nektar::SolverUtils::EquationSystem::m_spacedim, m_u0, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ TestSWE2Dproblem()
|
private |
Definition at line 1985 of file MMFSWE.cpp.
1987{
1988 boost::ignore_unused(time);
1989
1991
1992 Array<OneD, NekDouble> x0(nq);
1993 Array<OneD, NekDouble> x1(nq);
1994 Array<OneD, NekDouble> x2(nq);
1995
1996 m_fields[0]->GetCoords(x0, x1, x2);
1997
1998 Array<OneD, NekDouble> eta0(nq);
1999 Array<OneD, NekDouble> u0(nq);
2000 Array<OneD, NekDouble> v0(nq);
2001
2002 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2003
2005 {
2006 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2007 }
2008
2009 for (int i = 0; i < nq; ++i)
2010 {
2011 eta0[i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2012 (1.0 / cosh(0.395 * x0[i]))) *
2013 (3.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2014 exp(-0.5 * x1[i] * x1[i]);
2015 uvec[0][i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2016 (1.0 / cosh(0.395 * x0[i]))) *
2017 (-9.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2018 exp(-0.5 * x1[i] * x1[i]);
2019 uvec[1][i] = (-2.0 * 0.395 * tanh(0.395 * x0[i])) *
2020 (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2021 (1.0 / cosh(0.395 * x0[i]))) *
2022 (2.0 * x1[i]) * exp(-0.5 * x1[i] * x1[i]);
2023 }
2024
2025 u0 = CartesianToMovingframes(uvec, 0);
2026 v0 = CartesianToMovingframes(uvec, 1);
2027
2028 switch (field)
2029 {
2030 case (0):
2031 {
2032 outfield = eta0;
2033 }
2034 break;
2035
2036 case (1):
2037 {
2038 outfield = u0;
2039 }
2040 break;
2041
2042 case (2):
2043 {
2044 outfield = v0;
2045 }
2046 break;
2047 }
2048}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::EquationSystem::m_fields, and Nektar::SolverUtils::EquationSystem::m_spacedim.
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ TestVorticityComputation()
Definition at line 2946 of file MMFSWE.cpp.
2947{
2948 // Construct beta
2949 int i, k;
2950 int n = 1, m = 1;
2952
2953 NekDouble alpha, beta_theta, beta_phi;
2954
2955 NekDouble xp, yp, zp, Re;
2956 NekDouble theta, phi, sin_theta, cos_theta, sin_varphi, cos_varphi;
2957 NekDouble cosntheta3;
2958
2959 NekDouble thetax, thetay, thetaz;
2960 NekDouble phix, phiy, phiz;
2961
2962 Array<OneD, NekDouble> x(nq);
2963 Array<OneD, NekDouble> y(nq);
2964 Array<OneD, NekDouble> z(nq);
2965
2966 Array<OneD, NekDouble> u(nq);
2967 Array<OneD, NekDouble> v(nq);
2968
2969 Array<OneD, NekDouble> vorticitycompt(nq);
2970 Array<OneD, NekDouble> vorticityexact(nq);
2971
2973
2974 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2976 {
2977 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2978 }
2979
2980 // Generate SphericalCoords
2981 for (k = 0; k < nq; ++k)
2982 {
2983 xp = x[k];
2984 yp = y[k];
2985 zp = z[k];
2986
2987 Re = sqrt(xp * xp + yp * yp + zp * zp);
2988
2989 CartesianToSpherical(xp, yp, zp, sin_varphi, cos_varphi, sin_theta,
2990 cos_theta);
2991
2992 alpha = sin_varphi;
2993
2994 theta = atan2(sin_theta, cos_theta);
2995 phi = atan2(sin_varphi, cos_varphi);
2996
2997 cosntheta3 = cos(n * theta) * cos(n * theta) * cos(n * theta);
2998
2999 beta_theta = -4.0 * n * cosntheta3 * cos(m * phi) * sin(n * theta) / Re;
3000 beta_phi = -m * cosntheta3 * sin(m * phi) / Re;
3001
3002 thetax = -1.0 * cos_varphi * sin_theta;
3003 thetay = -1.0 * sin_varphi * sin_theta;
3004 thetaz = cos_theta;
3005
3006 phix = -1.0 * sin_varphi;
3007 phiy = cos_varphi;
3008 phiz = 0.0;
3009
3010 uvec[0][k] = alpha * (beta_theta * thetax + beta_phi * phix);
3011 uvec[1][k] = alpha * (beta_theta * thetay + beta_phi * phiy);
3012 uvec[2][k] = alpha * (beta_theta * thetaz + beta_phi * phiz);
3013
3014 vorticityexact[k] = -4.0 * n / Re / Re * cos_theta * cos_theta *
3015 cos_varphi * cos(m * phi) * sin(n * theta);
3016 }
3017
3018 u = CartesianToMovingframes(uvec, 0);
3019 v = CartesianToMovingframes(uvec, 1);
3020
3021 std::cout << "chi migi1" << std::endl;
3022
3023 ComputeVorticity(u, v, vorticitycompt);
3024 /*for (int k=0; k < nq; k++)
3025 {
3026
3027 std::cout << "vorticitycompt[ " << k << "]"<< "\t"<<vorticitycompt[k]<<
3028 std::endl;
3029
3030 }*/
3031
3032 Vmath::Vsub(nq, vorticityexact, 1, vorticitycompt, 1, vorticitycompt, 1);
3033
3036 << std::endl;
3037}
References Nektar::SolverUtils::MMFSystem::AvgAbsInt(), Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), ComputeVorticity(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, tinysimd::sqrt(), Vmath::Vamax(), Vmath::Vsub(), and Nektar::UnitTests::z().
Referenced by v_InitObject().
◆ UnstableJetFlow()
|
protected |
Definition at line 2482 of file MMFSWE.cpp.
2484{
2485 boost::ignore_unused(time);
2486
2488
2489 NekDouble uhat, vhat;
2490 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2491 NekDouble x0j, x1j, x2j;
2492 NekDouble Ttheta, Tphi;
2493
2494 Array<OneD, NekDouble> eta(nq, 0.0);
2495 Array<OneD, NekDouble> u(nq, 0.0);
2496 Array<OneD, NekDouble> v(nq, 0.0);
2497
2498 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2500 {
2501 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2502 }
2503
2504 Array<OneD, NekDouble> x(nq);
2505 Array<OneD, NekDouble> y(nq);
2506 Array<OneD, NekDouble> z(nq);
2507
2509
2510 int Nint = 1000;
2511 NekDouble dth, intj;
2512 for (int j = 0; j < nq; ++j)
2513 {
2514 x0j = x[j];
2515 x1j = y[j];
2516 x2j = z[j];
2517
2518 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2519 cos_theta);
2520
2521 Ttheta = atan2(sin_theta, cos_theta);
2522 Tphi = atan2(sin_varphi, cos_varphi);
2523
2524 uhat = ComputeUnstableJetuphi(Ttheta);
2525 vhat = 0.0;
2526
2527 // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2528 dth = Ttheta / Nint;
2529 eta[j] = dth * 0.5 *
2531 for (int i = 1; i < Nint - 1; i++)
2532 {
2533 intj = i * dth;
2534 eta[j] = eta[j] + dth * ComputeUnstableJetEta(intj);
2535 }
2536 eta[j] = (-1.0 / m_g) * eta[j];
2537
2538 // Add perturbation
2540 {
2541 eta[j] = eta[j] + m_hbar * cos_theta * exp(-9.0 * Tphi * Tphi) *
2542 exp(-225.0 * (m_pi / 4.0 - Ttheta) *
2543 (m_pi / 4.0 - Ttheta));
2544 }
2545
2546 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2547 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2548 uvec[2][j] = vhat * cos_theta;
2549 }
2550
2551 // Projection of u onto the tangent plane with conserving the mag. of the
2552 // velocity.
2553 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2554
2556 {
2557 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2558 }
2559
2560 // u is projected on the tangent plane with preserving its length
2561 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2562
2563 // Change it to the coordinate of moving frames
2564 // CartesianToMovingframes(0,uvecproj,u);
2565 // CartesianToMovingframes(1,uvecproj,v);
2566
2567 u = CartesianToMovingframes(uvec, 0);
2568 v = CartesianToMovingframes(uvec, 1);
2569
2570 switch (field)
2571 {
2572 case (0):
2573 {
2574 outfield = eta;
2575 }
2576 break;
2577
2578 case (1):
2579 {
2580 outfield = u;
2581 }
2582 break;
2583
2584 case (2):
2585 {
2586 outfield = v;
2587 }
2588 break;
2589 }
2590}
NekDouble ComputeUnstableJetEta(const NekDouble theta)
Definition: MMFSWE.cpp:2737
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), ComputeUnstableJetEta(), ComputeUnstableJetuphi(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_hbar, Nektar::SolverUtils::MMFSystem::m_pi, m_PurturbedJet, Nektar::SolverUtils::EquationSystem::m_spacedim, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ UnsteadyZonalFlow()
|
protected |
Definition at line 2292 of file MMFSWE.cpp.
2294{
2296 NekDouble uhat, vhat;
2297 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2298 NekDouble x0j, x1j, x2j, tmp;
2299
2300 NekDouble TR, Ttheta;
2301
2302 Array<OneD, NekDouble> eta(nq, 0.0);
2303 Array<OneD, NekDouble> u(nq, 0.0);
2304 Array<OneD, NekDouble> v(nq, 0.0);
2305
2306 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2308 {
2309 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2310 }
2311
2312 Array<OneD, NekDouble> x(nq);
2313 Array<OneD, NekDouble> y(nq);
2314 Array<OneD, NekDouble> z(nq);
2315
2317
2318 for (int j = 0; j < nq; ++j)
2319 {
2320 x0j = x[j];
2321 x1j = y[j];
2322 x2j = z[j];
2323
2324 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2325 cos_theta);
2326
2327 // \eta = ( - ( u_0 ( - T_R sin \alpha \cos \theta + \cos \alpha \sin
2328 // \theta ) + \Omega \sin \theta )^2 + \Omega \sin \theta )^2 + (\Omega
2329 // \sin \theta )^2
2330 TR =
2332 Ttheta =
2335
2339 eta[j] = 0.5 * eta[j] / m_g;
2340
2341 // u = u_0*(TR*\sin \alpha * \sin \theta + \cos \alpha * \cos \theta
2342 // v = - u_0 Ttheta * \sin \alpha
2343 uhat =
2346
2347 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2348 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2349 uvec[2][j] = vhat * cos_theta;
2350 }
2351
2352 // Projection of u onto the tangent plane with conserving the mag. of the
2353 // velocity.
2354 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2355
2357 {
2358 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2359 }
2360
2361 // u is projected on the tangent plane with preserving its length
2362 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2363
2364 // Change it to the coordinate of moving frames
2365 // CartesianToMovingframes(0,uvecproj,u);
2366 // CartesianToMovingframes(1,uvecproj,v);
2367
2368 u = CartesianToMovingframes(uvec, 0);
2369 v = CartesianToMovingframes(uvec, 1);
2370
2371 switch (field)
2372 {
2373 case (0):
2374 {
2375 outfield = eta;
2376 }
2377 break;
2378
2379 case (1):
2380 {
2381 outfield = u;
2382 }
2383 break;
2384
2385 case (2):
2386 {
2387 outfield = v;
2388 }
2389 break;
2390 }
2391}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_Omega, Nektar::SolverUtils::EquationSystem::m_spacedim, m_u0, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ v_DoInitialise()
|
overrideprotectedvirtual |
Sets up initial conditions.
Sets the initial conditions.
Reimplemented from Nektar::SolverUtils::UnsteadySystem.
Definition at line 1540 of file MMFSWE.cpp.
1541{
1542 // Compute m_depth and m_Derivdepth
1543 EvaluateWaterDepth();
1544 EvaluateCoriolis();
1545 SetInitialConditions(0.0, dumpInitialConditions);
1546 PrimitiveToConservative();
1547
1548 // transfer the initial conditions to modal values
1550 {
1553 m_fields[i]->UpdateCoeffs());
1554 }
1555}
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
Definition: EquationSystem.h:675
References EvaluateCoriolis(), EvaluateWaterDepth(), Nektar::SolverUtils::EquationSystem::m_fields, PrimitiveToConservative(), and Nektar::SolverUtils::EquationSystem::SetInitialConditions().
◆ v_DoSolve()
Solves an unsteady problem.
Initialises the time integration scheme (as specified in the session file), and perform the time integration.
Reimplemented from Nektar::SolverUtils::UnsteadySystem.
Definition at line 249 of file MMFSWE.cpp.
250{
252
253 int i, nchk = 1;
254 int nvariables = 0;
257
259 {
260 for (i = 0; i < nfields; ++i)
261 {
262 m_intVariables.push_back(i);
263 }
264 nvariables = nfields;
265 }
266 else
267 {
268 nvariables = m_intVariables.size();
269 }
270
271 // Set up wrapper to fields data storage.
272 Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
273 Array<OneD, Array<OneD, NekDouble>> tmp(nvariables);
274
275 // Order storage to list time-integrated fields first.
276 for (i = 0; i < nvariables; ++i)
277 {
280 }
281
282 // Initialise time integration scheme
284
285 // Check uniqueness of checkpoint output
289 "Only one of IO_CheckTime and IO_CheckSteps "
290 "should be set!");
291
292 LibUtilities::Timer timer;
293 bool doCheckTime = false;
294 int step = 0;
295 NekDouble intTime = 0.0;
296 NekDouble cpuTime = 0.0;
297 NekDouble elapsed = 0.0;
298
299 NekDouble Mass = 0.0, Energy = 0.0, Enstrophy = 0.0, Vorticity = 0.0;
300 Array<OneD, NekDouble> zeta(nq);
301 Array<OneD, Array<OneD, NekDouble>> fieldsprimitive(nvariables);
302 for (int i = 0; i < nvariables; ++i)
303 {
304 fieldsprimitive[i] = Array<OneD, NekDouble>(nq);
305 }
306
308 {
309 timer.Start();
311 timer.Stop();
312
314 elapsed = timer.TimePerTest(1);
315 intTime += elapsed;
316 cpuTime += elapsed;
317
318 // Write out status information
320 m_session->GetComm()->GetRank() == 0)
321 {
322 std::cout << "Steps: " << std::setw(8) << std::left << step + 1
323 << " "
325
326 std::stringstream ss;
327 ss << cpuTime << "s";
328 std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
329 << std::endl;
330
331 // Printout Mass, Energy, Enstrophy
332 ConservativeToPrimitive(fields, fieldsprimitive);
333
334 // Vorticity zeta
335 ComputeVorticity(fieldsprimitive[1], fieldsprimitive[2], zeta);
337
338 // Masss = h^*
340
341 // Energy = 0.5*( h^*(u^2 + v^2) + g ( h^2 - h_s^s ) )
342 Energy = (ComputeEnergy(fieldsprimitive[0], fieldsprimitive[1],
343 fieldsprimitive[2]) -
344 m_Energy0) /
345 m_Energy0;
346
347 // Enstrophy = 0.5/h^* ( \mathbf{k} \cdot (\nabla \times \mathbf{v}
348 // ) + f )^2
349 Enstrophy =
350 (ComputeEnstrophy(fieldsprimitive[0], fieldsprimitive[1],
351 fieldsprimitive[2]) -
352 m_Enstrophy0) /
353 m_Enstrophy0;
354
355 std::cout << "dMass = " << std::setw(8) << std::left << Mass << " "
356 << ", dEnergy = " << std::setw(8) << std::left << Energy
357 << " "
358 << ", dEnstrophy = " << std::setw(8) << std::left
359 << Enstrophy << " "
360 << ", dVorticity = " << std::setw(8) << std::left
361 << Vorticity << std::endl
362 << std::endl;
363
364 cpuTime = 0.0;
365 }
366
367 // (h, hu, hv) -> (\eta, u, v)
368 ConservativeToPrimitive();
369
370 // Transform data into coefficient space
371 for (i = 0; i < nvariables; ++i)
372 {
377 }
378
379 // Write out checkpoint files
381 doCheckTime)
382 {
383 Checkpoint_Output(nchk++);
384 doCheckTime = false;
385
386 // (\eta, u, v) -> (h, hu, hv)
387 PrimitiveToConservative();
388 }
389
390 // Step advance
391 ++step;
392 }
393
394 // Print out summary statistics
396 {
398 {
400 << std::endl
402 }
403
405 {
406 std::cout << "Time-integration : " << intTime << "s" << std::endl;
407 }
408 }
409
410 for (i = 0; i < nvariables; ++i)
411 {
414 }
415
416 // (h, hu, hv) -> (\eta, u, v)
417 ConservativeToPrimitive();
418
419 for (i = 0; i < nvariables; ++i)
420 {
422 m_fields[i]->UpdateCoeffs());
423 }
424}
NekDouble ComputeEnstrophy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2180
NekDouble ComputeMass(const Array< OneD, const NekDouble > &eta)
Definition: MMFSWE.cpp:2141
NekDouble ComputeEnergy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2151
int m_infosteps
Number of time steps between outputting status information.
Definition: EquationSystem.h:394
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition: EquationSystem.cpp:1176
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:79
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
Definition: UnsteadySystem.h:87
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
Definition: UnsteadySystem.h:77
std::vector< int > m_intVariables
Definition: UnsteadySystem.h:84
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
References tinysimd::abs(), ASSERTL0, Nektar::SolverUtils::EquationSystem::Checkpoint_Output(), ComputeEnergy(), ComputeEnstrophy(), ComputeMass(), ComputeVorticity(), ConservativeToPrimitive(), Nektar::NekConstants::kNekZeroTol, Nektar::SolverUtils::UnsteadySystem::m_cflSafetyFactor, Nektar::SolverUtils::EquationSystem::m_checksteps, Nektar::SolverUtils::EquationSystem::m_checktime, m_Energy0, m_Enstrophy0, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_fintime, Nektar::SolverUtils::EquationSystem::m_infosteps, Nektar::SolverUtils::UnsteadySystem::m_intScheme, Nektar::SolverUtils::UnsteadySystem::m_intVariables, m_Mass0, Nektar::SolverUtils::UnsteadySystem::m_ode, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_steps, Nektar::SolverUtils::EquationSystem::m_time, Nektar::SolverUtils::EquationSystem::m_timestep, m_Vorticity0, PrimitiveToConservative(), Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), and Nektar::LibUtilities::Timer::TimePerTest().
◆ v_EvaluateExactSolution()
|
overrideprotectedvirtual |
Reimplemented from Nektar::SolverUtils::EquationSystem.
Definition at line 3228 of file MMFSWE.cpp.
3231{
3233 {
3235 {
3236 SteadyZonalFlow(field, outfield);
3237 }
3238 break;
3239
3241 {
3242 UnsteadyZonalFlow(field, time, outfield);
3243 }
3244 break;
3245
3247 {
3248 IsolatedMountainFlow(field, time, outfield);
3249 }
3250 break;
3251
3253 {
3254 UnstableJetFlow(field, time, outfield);
3255 }
3256 break;
3257
3259 {
3260 RossbyWave(field, outfield);
3261 }
3262 break;
3263
3265 {
3266 TestSWE2Dproblem(time, field, outfield);
3267 }
3268 break;
3269
3270 default:
3271 break;
3272 }
3273}
void UnstableJetFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2482
void RossbyWave(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2592
void IsolatedMountainFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2393
void UnsteadyZonalFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2292
void TestSWE2Dproblem(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:1985
void SteadyZonalFlow(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2050
References Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, IsolatedMountainFlow(), m_TestType, RossbyWave(), SteadyZonalFlow(), TestSWE2Dproblem(), UnstableJetFlow(), and UnsteadyZonalFlow().
Referenced by v_L2Error().
◆ v_GenerateSummary()
|
overrideprotectedvirtual |
Print Summary.
Reimplemented from Nektar::SolverUtils::MMFSystem.
Definition at line 3275 of file MMFSWE.cpp.
3276{
3277 MMFSystem::v_GenerateSummary(s);
3279
3281 {
3283 {
3285 }
3286 break;
3287
3289 {
3291 m_RossbyDisturbance);
3292 }
3293 break;
3294
3296 {
3298 }
3299 break;
3300
3301 default:
3302 break;
3303 }
3304}
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
Definition: MMFSystem.cpp:2452
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:49
References Nektar::SolverUtils::AddSummaryItem(), Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, m_alpha, m_PurturbedJet, m_RossbyDisturbance, m_TestType, Nektar::TestTypeMap, and Nektar::SolverUtils::MMFSystem::v_GenerateSummary().
◆ v_InitObject()
|
overrideprotectedvirtual |
Initialise the object.
Initialisation object for the unsteady linear advection equation.
Reimplemented from Nektar::SolverUtils::UnsteadySystem.
Definition at line 62 of file MMFSWE.cpp.
63{
64 // Call to the initialisation object
65 UnsteadySystem::v_InitObject(DeclareFields);
66
69 Array<OneD, Array<OneD, NekDouble>> Anisotropy(shapedim);
70 for (int j = 0; j < shapedim; ++j)
71 {
72 Anisotropy[j] = Array<OneD, NekDouble>(nq, 1.0);
73 }
74
75 MMFSystem::MMFInitObject(Anisotropy);
76
77 // Load acceleration of gravity
79
80 // Add Coriolois effects
82
83 // Add Rotation of the sphere along the pole
85
88
89 // Define TestType
91 "No TESTTYPE defined in session.");
94 {
96 {
98 break;
99 }
100 }
101
102 // Variable Setting for each test problem
103 NekDouble gms = 9.80616;
104
106 {
108 {
109 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
110 NekDouble rad_earth = 6.37122 * 1000000;
111 NekDouble Omegams;
112
113 // Nondimensionalized coeffs.
114 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
115
119 m_Omega = Omegams * SecondToDay;
120
123
125 }
126 break;
127
129 {
130 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
131 NekDouble rad_earth = 6.37122 * 1000000;
132 NekDouble Omegams;
133
134 // Nondimensionalized coeffs.
135 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
136
140 m_Omega = Omegams * SecondToDay;
141
144 }
145 break;
146
148 {
149 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
150 NekDouble rad_earth = 6.37122 * 1000000;
151 NekDouble Omegams;
152
153 // Nondimensionalized coeffs.
154 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
155
157
160
162 m_Omega = Omegams * SecondToDay;
163
164 m_H0 = 5960.0 / rad_earth;
165 m_hs0 = 2000.0 / rad_earth;
166 }
167 break;
168
170 {
171 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
172 NekDouble rad_earth = 6.37122 * 1000000;
173 NekDouble Omegams;
174
175 // Nondimensionalized coeffs.
176 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
177
180
182 m_Omega = Omegams * SecondToDay;
183
186
187 m_uthetamax = 80 * SecondToDay / rad_earth;
191 m_hbar = 120.0 / rad_earth;
192
196 }
197 break;
198
200 {
201 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
202 NekDouble rad_earth = 6.37122 * 1000000;
203 NekDouble Omegams;
204
205 // Nondimensionalized coeffs.
206 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
207
209 m_Omega = Omegams * SecondToDay;
210
213
214 m_angfreq = 7.848 * 0.000001 * SecondToDay;
215 m_K = 7.848 * 0.000001 * SecondToDay;
216 }
217 break;
218
219 default:
220 break;
221 }
222
223 // TestVorticityComputation
225 {
226 TestVorticityComputation();
227 }
228
229 // If explicit it computes RHS and PROJECTION for the time integration
231 {
234 }
235 // Otherwise it gives an error (no implicit integration)
236 else
237 {
239 }
240}
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:96
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:88
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the RHS.
Definition: MMFSWE.cpp:426
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Compute the projection.
Definition: MMFSWE.cpp:1358
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:53
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
Definition: UnsteadySystem.h:99
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:87
References ASSERTL0, Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineProjection(), DoOdeProjection(), DoOdeRhs(), Nektar::SolverUtils::eSphere, Nektar::eTestIsolatedMountain, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, m_AddCoriolis, m_AddRotation, m_alpha, m_angfreq, m_en, Nektar::SolverUtils::UnsteadySystem::m_explicitAdvection, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, m_hbar, m_hs0, m_Hvar, m_K, m_k2, Nektar::SolverUtils::UnsteadySystem::m_ode, m_Omega, Nektar::SolverUtils::MMFSystem::m_pi, m_PurturbedJet, m_RossbyDisturbance, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::MMFSystem::m_surfaceType, m_TestType, m_theta0, m_theta1, m_u0, m_uthetamax, Nektar::SolverUtils::MMFSystem::MMFInitObject(), Nektar::SIZE_TestType, Nektar::TestTypeMap, TestVorticityComputation(), and Nektar::SolverUtils::UnsteadySystem::v_InitObject().
◆ v_L2Error()
|
overrideprotectedvirtual |
Virtual function for the L_2 error computation between fields and a given exact solution.
Compute the error in the L2-norm.
- Parameters
-
field The field to compare. exactsoln The exact solution to compare with. Normalised Normalise L2-error.
- Returns
- Error in the L2-norm.
Reimplemented from Nektar::SolverUtils::EquationSystem.
Definition at line 3039 of file MMFSWE.cpp.
3042{
3043 boost::ignore_unused(exactsoln);
3044
3046 NekDouble L2error = -1.0;
3047
3049 {
3051 {
3053 m_fields[field]->UpdatePhys());
3054 }
3055
3056 switch (field)
3057 {
3058 case (0):
3059 {
3060 // I(h) = I (h - h_exact ) / I (h_exact)
3061 Array<OneD, NekDouble> exactsolution(nq);
3063
3064 // exactsoln = u - u_T so that L2 compute u_T
3066
3068 &exactsolution[0], 1, &exactsolution[0], 1);
3069 Vmath::Vabs(nq, exactsolution, 1, exactsolution, 1);
3070
3071 L2error = (m_fields[0]->Integral(exactsolution)) / L2exact;
3072 }
3073 break;
3074
3075 case (1):
3076 {
3077 // I2 (u) = I( (u - u_ext)^2 + (v - v_ext)^2 )^{1/2} / I(
3078 // u_ext^2 + v_ext^2 )^{1/2}
3079 Array<OneD, NekDouble> exactu(nq);
3080 Array<OneD, NekDouble> exactv(nq);
3081 Array<OneD, NekDouble> tmp(nq);
3082
3083 // L2exact = \int (\sqrt{exactu*exactu+exactv*exactv})
3084 NekDouble L2exact;
3087 Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3088 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3089 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3090
3091 L2exact = m_fields[1]->Integral(tmp);
3092
3093 // L2exact = \int
3094 // (\sqrt{(u-exactu)*(u-exactu)+(v-exactv)*(v-exactv)})
3096 &exactu[0], 1);
3098 &exactv[0], 1);
3099 Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3100 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3101 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3102
3103 L2error = (m_fields[1]->Integral(tmp)) / L2exact;
3104 }
3105 break;
3106
3107 case (2):
3108 {
3109 L2error = 0.0;
3110 }
3111 break;
3112
3113 default:
3114 break;
3115 }
3116
3117 if (Normalised == true)
3118 {
3120
3123
3124 L2error = sqrt(L2error * L2error / Vol);
3125 }
3126 }
3127 else
3128 {
3129 Array<OneD, NekDouble> L2INF(2);
3130 L2INF = ErrorExtraPoints(field);
3131 L2error = L2INF[0];
3132 }
3133
3134 return L2error;
3135}
virtual void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
Definition: MMFSWE.cpp:3228
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > ErrorExtraPoints(unsigned int field)
Compute error (L2 and L_inf) over an larger set of quadrature points return [L2 Linf].
Definition: EquationSystem.cpp:895
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:775
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Definition: EquationSystem.h:443
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:529
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:548
References Nektar::SolverUtils::EquationSystem::ErrorExtraPoints(), Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::SolverUtils::EquationSystem::m_comm, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError, Nektar::SolverUtils::EquationSystem::m_time, Nektar::LibUtilities::ReduceSum, tinysimd::sqrt(), v_EvaluateExactSolution(), Vmath::Vabs(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().
◆ v_LinfError()
|
overrideprotectedvirtual |
Virtual function for the L_inf error computation between fields and a given exact solution.
Compute the error in the L_inf-norm
- Parameters
-
field The field to compare. exactsoln The exact solution to compare with.
- Returns
- Error in the L_inft-norm.
Reimplemented from Nektar::SolverUtils::EquationSystem.
Definition at line 3137 of file MMFSWE.cpp.
3139{
3140 boost::ignore_unused(exactsoln);
3141
3143
3145 {
3147 m_fields[field]->UpdatePhys());
3148 }
3149
3151
3152 // Obtain \vec{u} in cartesian coordinate
3153 Array<OneD, NekDouble> x(nq);
3154 Array<OneD, NekDouble> y(nq);
3155 Array<OneD, NekDouble> z(nq);
3156
3158
3159 switch (field)
3160 {
3161 case (0):
3162 {
3163 NekDouble Letaint;
3164
3165 Array<OneD, NekDouble> exactsolution(nq);
3166
3169 exactsolution);
3170
3171 Letaint = Vmath::Vamax(nq, exactsolution, 1);
3172
3174 1, &exactsolution[0], 1);
3175
3177 }
3178 break;
3179
3180 case (1):
3181 {
3182 Array<OneD, NekDouble> exactu(nq);
3183 Array<OneD, NekDouble> exactv(nq);
3184 Array<OneD, NekDouble> tmpu(nq);
3185 Array<OneD, NekDouble> tmpv(nq);
3186 Array<OneD, NekDouble> Lerr(nq);
3187 Array<OneD, NekDouble> uT(nq);
3188
3191
3192 // Compute max[sqrt{(u-uex)^2 + (v-vex)^2}]
3195
3196 Vmath::Vsub(nq, &exactu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3197 Vmath::Vsub(nq, &exactv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3198
3199 Vmath::Vmul(nq, &tmpu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3200 Vmath::Vmul(nq, &tmpv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3201
3202 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &Lerr[0], 1);
3203 Vmath::Vsqrt(nq, &Lerr[0], 1, &Lerr[0], 1);
3204
3205 // uT = max[sqrt( u_T^2 + v_T^2 ) ]
3206 Vmath::Vmul(nq, &exactu[0], 1, &exactu[0], 1, &tmpu[0], 1);
3207 Vmath::Vmul(nq, &exactv[0], 1, &exactv[0], 1, &tmpv[0], 1);
3208 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &uT[0], 1);
3209 Vmath::Vsqrt(nq, &uT[0], 1, &uT[0], 1);
3210
3212 }
3213 break;
3214
3215 case (2):
3216 {
3217 LinfError = 0.0;
3218 }
3219 break;
3220
3221 default:
3222 break;
3223 }
3224
3226}
SOLVER_UTILS_EXPORT NekDouble LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Linf error computation.
Definition: EquationSystem.h:600
SOLVER_UTILS_EXPORT void EvaluateExactSolution(int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Evaluates an exact solution.
Definition: EquationSystem.h:683
References Nektar::SolverUtils::EquationSystem::EvaluateExactSolution(), Nektar::SolverUtils::EquationSystem::LinfError(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_time, Vmath::Vadd(), Vmath::Vamax(), Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Nektar::UnitTests::z().
◆ v_SetInitialConditions()
|
overrideprotectedvirtual |
Set the physical fields based on a restart file, or a function describing the initial condition given in the session.
- Parameters
-
initialtime Time at which to evaluate the function. dumpInitialConditions Write the initial condition to file?
Reimplemented from Nektar::SolverUtils::EquationSystem.
Definition at line 1786 of file MMFSWE.cpp.
1789{
1790 boost::ignore_unused(domain);
1791
1793
1795 {
1797 {
1798 Array<OneD, NekDouble> eta0(nq);
1799 Array<OneD, NekDouble> u0(nq);
1800 Array<OneD, NekDouble> v0(nq);
1801
1802 TestSWE2Dproblem(initialtime, 0, eta0);
1803 m_fields[0]->SetPhys(eta0);
1804
1805 TestSWE2Dproblem(initialtime, 1, u0);
1806 m_fields[1]->SetPhys(u0);
1807
1808 TestSWE2Dproblem(initialtime, 2, v0);
1809 m_fields[2]->SetPhys(v0);
1810
1811 // forward transform to fill the modal coeffs
1813 {
1816 m_fields[i]->UpdateCoeffs());
1817 }
1818 }
1819 break;
1820
1822 {
1823 Array<OneD, NekDouble> eta0(nq);
1824 Array<OneD, NekDouble> u0(nq);
1825 Array<OneD, NekDouble> v0(nq);
1826 Array<OneD, NekDouble> zeta0(nq);
1827
1828 SteadyZonalFlow(0, eta0);
1829 m_fields[0]->SetPhys(eta0);
1830
1831 SteadyZonalFlow(1, u0);
1832 m_fields[1]->SetPhys(u0);
1833
1834 SteadyZonalFlow(2, v0);
1835 m_fields[2]->SetPhys(v0);
1836
1837 // ComputeVorticity(u0, v0, zeta0);
1839
1843
1844 // forward transform to fill the modal coeffs
1846 {
1849 m_fields[i]->UpdateCoeffs());
1850 }
1851 }
1852 break;
1853
1855 {
1856 Array<OneD, NekDouble> eta0(nq);
1857 Array<OneD, NekDouble> u0(nq);
1858 Array<OneD, NekDouble> v0(nq);
1859
1860 UnsteadyZonalFlow(0, initialtime, eta0);
1861 m_fields[0]->SetPhys(eta0);
1862
1863 UnsteadyZonalFlow(1, initialtime, u0);
1864 m_fields[1]->SetPhys(u0);
1865
1866 UnsteadyZonalFlow(2, initialtime, v0);
1867 m_fields[2]->SetPhys(v0);
1868
1872
1873 // forward transform to fill the modal coeffs
1875 {
1878 m_fields[i]->UpdateCoeffs());
1879 }
1880 }
1881 break;
1882
1884 {
1885 Array<OneD, NekDouble> eta0(nq);
1886 Array<OneD, NekDouble> u0(nq);
1887 Array<OneD, NekDouble> v0(nq);
1888
1889 IsolatedMountainFlow(0, initialtime, eta0);
1890 m_fields[0]->SetPhys(eta0);
1891
1892 IsolatedMountainFlow(1, initialtime, u0);
1893 m_fields[1]->SetPhys(u0);
1894
1895 IsolatedMountainFlow(2, initialtime, v0);
1896 m_fields[2]->SetPhys(v0);
1897
1901
1902 // forward transform to fill the modal coeffs
1904 {
1907 m_fields[i]->UpdateCoeffs());
1908 }
1909 }
1910 break;
1911
1913 {
1914 Array<OneD, NekDouble> eta0(nq);
1915 Array<OneD, NekDouble> u0(nq);
1916 Array<OneD, NekDouble> v0(nq);
1917
1918 UnstableJetFlow(0, initialtime, eta0);
1919 m_fields[0]->SetPhys(eta0);
1920
1921 UnstableJetFlow(1, initialtime, u0);
1922 m_fields[1]->SetPhys(u0);
1923
1924 UnstableJetFlow(2, initialtime, v0);
1925 m_fields[2]->SetPhys(v0);
1926
1930
1931 // forward transform to fill the modal coeffs
1933 {
1936 m_fields[i]->UpdateCoeffs());
1937 }
1938 }
1939 break;
1940
1942 {
1943 Array<OneD, NekDouble> eta0(nq);
1944 Array<OneD, NekDouble> u0(nq);
1945 Array<OneD, NekDouble> v0(nq);
1946
1947 RossbyWave(0, eta0);
1948 m_fields[0]->SetPhys(eta0);
1949
1950 RossbyWave(1, u0);
1951 m_fields[1]->SetPhys(u0);
1952
1953 RossbyWave(2, v0);
1954 m_fields[2]->SetPhys(v0);
1955
1959
1960 // forward transform to fill the modal coeffs
1962 {
1965 m_fields[i]->UpdateCoeffs());
1966 }
1967 }
1968 break;
1969
1970 default:
1971 break;
1972 }
1973
1974 if (dumpInitialConditions)
1975 {
1976 // dump initial conditions to file
1978 WriteFld(outname);
1979
1981 Checkpoint_Output_Cartesian(outname);
1982 }
1983}
void Checkpoint_Output_Cartesian(std::string outname)
Definition: MMFSWE.cpp:2767
References Checkpoint_Output_Cartesian(), ComputeEnergy(), ComputeEnstrophy(), ComputeMass(), Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, Nektar::SolverUtils::EquationSystem::GetTotPoints(), IsolatedMountainFlow(), m_Energy0, m_Enstrophy0, Nektar::SolverUtils::EquationSystem::m_fields, m_Mass0, Nektar::SolverUtils::EquationSystem::m_sessionName, m_TestType, m_Vorticity0, RossbyWave(), SteadyZonalFlow(), TestSWE2Dproblem(), UnstableJetFlow(), UnsteadyZonalFlow(), and Nektar::SolverUtils::EquationSystem::WriteFld().
◆ WallBoundary2D()
|
protected |
Definition at line 1427 of file MMFSWE.cpp.
1429{
1430
1431 int i;
1433 int nvariables = physarray.size();
1434
1435 // get physical values of the forward trace
1436 Array<OneD, Array<OneD, NekDouble>> Fwd0(nvariables);
1437 for (i = 0; i < nvariables; ++i)
1438 {
1439 Fwd0[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1440 m_fields[i]->ExtractTracePhys(physarray[i], Fwd0[i]);
1441 }
1442
1443 Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
1444 for (i = 0; i < nvariables; ++i)
1445 {
1446 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1447 m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1448 }
1449
1450 // Adjust the physical values of the trace to take
1451 // user defined boundaries into account
1452 int e, id1, id2, npts;
1453
1455 ++e)
1456 {
1457 npts = m_fields[0]
1458 ->GetBndCondExpansions()[bcRegion]
1459 ->GetExp(e)
1460 ->GetTotPoints();
1461 id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
1462 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
1463 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1464
1466 {
1467 case 1:
1468 {
1469 // negate the forward flux
1470 Vmath::Neg(npts, &Fwd[1][id2], 1);
1471 }
1472 break;
1473 case 2:
1474 {
1475 Array<OneD, NekDouble> tmp_n(npts);
1476 Array<OneD, NekDouble> tmp_t(npts);
1477
1479 &tmp_n[0], 1);
1481 &tmp_n[0], 1, &tmp_n[0], 1);
1482
1484 &tmp_t[0], 1);
1486 1, &tmp_t[0], 1, &tmp_t[0], 1);
1487
1488 // negate the normal flux
1489 Vmath::Neg(npts, tmp_n, 1);
1490
1491 Array<OneD, NekDouble> denom(npts);
1492 Array<OneD, NekDouble> tmp_u(npts);
1493 Array<OneD, NekDouble> tmp_v(npts);
1494
1495 // denom = (e^1 \cdot n ) (e^2 \cdot t) - (e^2 \cdot n ) (e^1
1496 // \cdot t)
1498 &m_nperpcdotMFFwd[0][id2], 1, &denom[0], 1);
1500 &m_nperpcdotMFFwd[1][id2], 1, &denom[0], 1,
1501 &denom[0], 1);
1502
1504 &tmp_u[0], 1);
1506 &tmp_u[0], 1, &tmp_u[0], 1);
1507 Vmath::Vdiv(npts, &tmp_u[0], 1, &denom[0], 1, &tmp_u[0], 1);
1508
1509 Vmath::Vcopy(npts, &tmp_u[0], 1, &Fwd[1][id2], 1);
1510
1512 &tmp_v[0], 1);
1514 &tmp_v[0], 1, &tmp_v[0], 1);
1515 Vmath::Vdiv(npts, &tmp_v[0], 1, &denom[0], 1, &tmp_v[0], 1);
1516
1517 Vmath::Vcopy(npts, &tmp_v[0], 1, &Fwd[2][id2], 1);
1518 }
1519 break;
1520
1521 default:
1523 }
1524
1525 // copy boundary adjusted values into the boundary expansion
1526 for (i = 0; i < nvariables; ++i)
1527 {
1528 Vmath::Vcopy(npts, &Fwd[i][id2], 1,
1529 &(m_fields[i]
1530 ->GetBndCondExpansions()[bcRegion]
1531 ->UpdatePhys())[id1],
1532 1);
1533 }
1534 }
1535}
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:750
void Vvtvm(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)
vvtvm (vector times vector minus vector): z = w*x - y
Definition: Vmath.cpp:593
References ASSERTL0, Nektar::SolverUtils::EquationSystem::GetExpSize(), Nektar::SolverUtils::EquationSystem::GetTraceTotPoints(), Nektar::SolverUtils::EquationSystem::m_expdim, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd, Nektar::SolverUtils::MMFSystem::m_nperpcdotMFFwd, Vmath::Neg(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vvtvm(), and Vmath::Vvtvp().
Referenced by SetBoundaryConditions().
◆ WeakDGSWEDirDeriv()
|
protected |
Definition at line 482 of file MMFSWE.cpp.
485{
486 int i;
491
492 Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
493 Array<OneD, Array<OneD, NekDouble>> physfield(nvariables);
494
496 {
497 fluxvector[i] = Array<OneD, NekDouble>(nq);
498 }
499
500 // InField is Primitive
501 for (i = 0; i < nvariables; ++i)
502 {
503 physfield[i] = InField[i];
504 }
505
506 // Get the ith component of the flux vector in (physical space)
507 // fluxvector[0] = component for e^1 cdot \nabla \varphi
508 // fluxvector[1] = component for e^2 cdot \nabla \varphi
509 Array<OneD, NekDouble> fluxtmp(nq);
510 Array<OneD, NekDouble> tmp(ncoeffs);
511
512 // Compute Divergence Components
513 for (i = 0; i < nvariables; ++i)
514 {
515 GetSWEFluxVector(i, physfield, fluxvector);
516
517 OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
519 {
520 // Directional derivation with respect to the j'th moving frame
521 // tmp_j = ( \nabla \phi, fluxvector[j] \mathbf{e}^j )
523 fluxvector[j], tmp);
524 Vmath::Vadd(ncoeffs, &tmp[0], 1, &OutField[i][0], 1,
525 &OutField[i][0], 1);
526 }
527 }
528
529 // Numerical Flux
530 Array<OneD, Array<OneD, NekDouble>> numfluxFwd(nvariables);
531 Array<OneD, Array<OneD, NekDouble>> numfluxBwd(nvariables);
532
533 for (i = 0; i < nvariables; ++i)
534 {
535 numfluxFwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
536 numfluxBwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
537 }
538
539 NumericalSWEFlux(physfield, numfluxFwd, numfluxBwd);
540
541 // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
542 for (i = 0; i < nvariables; ++i)
543 {
544 Vmath::Neg(ncoeffs, OutField[i], 1);
545 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
546 OutField[i]);
548 }
549}
void NumericalSWEFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd)
Definition: MMFSWE.cpp:677
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:745
References Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetNpoints(), GetSWEFluxVector(), Nektar::SolverUtils::EquationSystem::GetTraceNpoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Neg(), NumericalSWEFlux(), and Vmath::Vadd().
Referenced by DoOdeRhs().
Friends And Related Function Documentation
◆ MemoryManager< MMFSWE >
|
friend |
Member Data Documentation
◆ className
|
static |
Initial value:
=
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: MMFSWE.h:67
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:103
Name of class.
◆ m_AddCoriolis
|
protected |
Definition at line 92 of file MMFSWE.h.
Referenced by DoOdeRhs(), and v_InitObject().
◆ m_AddRotation
|
protected |
Definition at line 92 of file MMFSWE.h.
Referenced by DoOdeRhs(), and v_InitObject().
◆ m_alpha
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by EvaluateCoriolisForZonalFlow(), SteadyZonalFlow(), UnsteadyZonalFlow(), v_GenerateSummary(), and v_InitObject().
◆ m_angfreq
|
protected |
Definition at line 99 of file MMFSWE.h.
Referenced by RossbyWave(), and v_InitObject().
◆ m_coriolis
Coriolis force.
Definition at line 90 of file MMFSWE.h.
Referenced by AddCoriolis(), ComputeEnstrophy(), and EvaluateCoriolis().
◆ m_depth
Still water depth.
Definition at line 86 of file MMFSWE.h.
Referenced by AddCoriolis(), AddElevationEffect(), AddRotation(), Checkpoint_Output_Cartesian(), ComputeEnergy(), ComputeEnstrophy(), ComputeMass(), ConservativeToPrimitive(), EvaluateWaterDepth(), GetSWEFluxVector(), NumericalSWEFlux(), and PrimitiveToConservative().
◆ m_Derivdepth
Definition at line 87 of file MMFSWE.h.
Referenced by AddElevationEffect(), and EvaluateWaterDepth().
◆ m_en
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_Energy0
|
protected |
Definition at line 101 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_Enstrophy0
|
protected |
Definition at line 101 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_g
|
protected |
Definition at line 94 of file MMFSWE.h.
Referenced by AddElevationEffect(), AverageFlux(), ComputeEnergy(), Computehhuhvflux(), EvaluateWaterDepth(), GetSWEFluxVector(), IsolatedMountainFlow(), LaxFriedrichFlux(), RiemannSolverHLLC(), RossbyWave(), RusanovFlux(), UnstableJetFlow(), UnsteadyZonalFlow(), and v_InitObject().
◆ m_H0
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by Checkpoint_Output_Cartesian(), ComputeEnergy(), EvaluateWaterDepth(), RossbyWave(), SteadyZonalFlow(), and v_InitObject().
◆ m_hbar
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by UnstableJetFlow(), and v_InitObject().
◆ m_hs0
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by EvaluateWaterDepth(), and v_InitObject().
◆ m_Hvar
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by SteadyZonalFlow(), and v_InitObject().
◆ m_K
|
protected |
Definition at line 99 of file MMFSWE.h.
Referenced by RossbyWave(), and v_InitObject().
◆ m_k2
|
protected |
Definition at line 96 of file MMFSWE.h.
Referenced by EvaluateWaterDepth(), and v_InitObject().
◆ m_Mass0
|
protected |
Definition at line 101 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_Omega
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by ComputeUnstableJetEta(), EvaluateCoriolisForZonalFlow(), EvaluateStandardCoriolis(), EvaluateWaterDepth(), IsolatedMountainFlow(), RossbyWave(), UnsteadyZonalFlow(), and v_InitObject().
◆ m_planeNumber
|
protected |
◆ m_primitive
|
protected |
◆ m_PurturbedJet
|
private |
Definition at line 276 of file MMFSWE.h.
Referenced by UnstableJetFlow(), v_GenerateSummary(), and v_InitObject().
◆ m_RossbyDisturbance
|
private |
Definition at line 275 of file MMFSWE.h.
Referenced by RossbyWave(), v_GenerateSummary(), and v_InitObject().
◆ m_TestType
| TestType Nektar::MMFSWE::m_TestType |
Definition at line 79 of file MMFSWE.h.
Referenced by EvaluateCoriolis(), EvaluateWaterDepth(), v_EvaluateExactSolution(), v_GenerateSummary(), v_InitObject(), and v_SetInitialConditions().
◆ m_theta0
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_theta1
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_traceVn
◆ m_u0
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by IsolatedMountainFlow(), SteadyZonalFlow(), UnsteadyZonalFlow(), and v_InitObject().
◆ m_uthetamax
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_veldotMF
◆ m_vellc
◆ m_velocity
Advection velocity.
Definition at line 107 of file MMFSWE.h.
Referenced by ComputeNablaCdotVelocity().
◆ m_Vorticity0
|
protected |
Definition at line 101 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
