#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... | |
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 | |
| TestType | m_TestType |
| 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, 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_movingFrameData |
| Moving reference frame status in the inertial frame X, Y, Z, Theta_x, Theta_y, Theta_z, U, V, W, Omega_x, Omega_y, Omega_z, A_x, A_y, A_z, DOmega_x, DOmega_y, DOmega_z, pivot_x, pivot_y, pivot_z. More... | |
| std::vector< std::string > | m_strFrameData |
| variable name in m_movingFrameData 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... | |
| Protected Attributes inherited from Nektar::SolverUtils::ALEHelper | |
| Array< OneD, MultiRegions::ExpListSharedPtr > | m_fieldsALE |
| Array< OneD, Array< OneD, NekDouble > > | m_gridVelocity |
| Array< OneD, Array< OneD, NekDouble > > | m_gridVelocityTrace |
| std::vector< ALEBaseShPtr > | m_ALEs |
| bool | m_ALESolver = false |
| bool | m_ImplicitALESolver = false |
| NekDouble | m_prevStageTime = 0.0 |
| int | m_spaceDim |
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 > |
Detailed Description
Constructor & Destructor Documentation
◆ MMFSWE()
|
protected |
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:41
SOLVER_UTILS_EXPORT UnsteadySystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Initialises UnsteadySystem class members.
Definition: UnsteadySystem.cpp:74
References m_planeNumber.
◆ ~MMFSWE()
|
overrideprotecteddefault |
Member Function Documentation
◆ AddCoriolis()
|
protected |
Definition at line 1207 of file MMFSWE.cpp.
1209{
1210 int ncoeffs = outarray[0].size();
1211 int nq = physarray[0].size();
1212
1213 Array<OneD, NekDouble> h(nq);
1214 Array<OneD, NekDouble> tmp(nq);
1215 Array<OneD, NekDouble> tmpc(ncoeffs);
1216
1217 // physarray is primitive
1218 // conservative formulation compute h
1219 // h = \eta + d
1221
1222 int indx = 0;
1224 {
1225 if (j == 0)
1226 {
1227 indx = 2;
1228 }
1229
1230 else if (j == 1)
1231 {
1232 indx = 1;
1233 }
1234
1235 // add to hu equation
1237 Vmath::Vmul(nq, h, 1, tmp, 1, tmp, 1);
1238
1239 if (j == 1)
1240 {
1241 Vmath::Neg(nq, tmp, 1);
1242 }
1243
1244 // N \cdot (e^1 \times e^2 )
1245 // Vmath::Vmul(nq, &m_MF1crossMF2dotSN[0], 1, &tmp[0], 1, &tmp[0], 1);
1246 m_fields[0]->IProductWRTBase(tmp, tmpc);
1247 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1248 }
1249}
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:467
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.hpp:72
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.hpp:180
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 543 of file MMFSWE.cpp.
545{
546 // routine works for both primitive and conservative formulations
547 int ncoeffs = outarray[0].size();
548 int nq = physarray[0].size();
549
550 Array<OneD, NekDouble> h(nq);
551 Array<OneD, NekDouble> tmp(nq);
552 Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
554 {
555 fluxvector[i] = Array<OneD, NekDouble>(nq);
556 }
557
558 // Get 0.5 g H*H / || e^m ||^2
559 GetSWEFluxVector(3, physarray, fluxvector);
560
561 Array<OneD, NekDouble> tmpc(ncoeffs);
563 {
565
566 Vmath::Neg(nq, &tmp[0], 1);
567 m_fields[0]->IProductWRTBase(tmp, tmpc);
568
569 Vmath::Vadd(ncoeffs, outarray[j + 1], 1, tmpc, 1, outarray[j + 1], 1);
570 }
571}
void GetSWEFluxVector(const int i, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: MMFSWE.cpp:573
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Definition: MMFSystem.h:192
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 1252 of file MMFSWE.cpp.
1254{
1255 int ncoeffs = outarray[0].size();
1256 int nq = physarray[0].size();
1257
1258 Array<OneD, NekDouble> h(nq);
1259 Array<OneD, NekDouble> tmp(nq);
1260 Array<OneD, NekDouble> tmpc(ncoeffs);
1261
1262 // physarray is primitive
1263 // conservative formulation compute h
1264 // h = \eta + d
1266
1268 {
1271
1272 m_fields[0]->IProductWRTBase(tmp, tmpc);
1273
1274 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1275 }
1276}
Array< OneD, Array< OneD, NekDouble > > m_Derivdepth
Definition: MMFSWE.h:83
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.hpp:100
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 1281 of file MMFSWE.cpp.
1283{
1284 // routine works for both primitive and conservative formulations
1285 int ncoeffs = outarray[0].size();
1286 int nq = physarray[0].size();
1287
1288 // Compute h
1289 Array<OneD, NekDouble> h(nq);
1291
1292 Array<OneD, NekDouble> de0dt_cdot_e0;
1293 Array<OneD, NekDouble> de0dt_cdot_e1;
1294 Array<OneD, NekDouble> de1dt_cdot_e0;
1295 Array<OneD, NekDouble> de1dt_cdot_e1;
1296 Compute_demdt_cdot_ek(0, 0, physarray, de0dt_cdot_e0);
1297 Compute_demdt_cdot_ek(1, 0, physarray, de1dt_cdot_e0);
1298 Compute_demdt_cdot_ek(0, 1, physarray, de0dt_cdot_e1);
1299 Compute_demdt_cdot_ek(1, 1, physarray, de1dt_cdot_e1);
1300
1301 Array<OneD, NekDouble> Rott1(nq);
1302 Array<OneD, NekDouble> Rott2(nq);
1303 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e0, 1, Rott1, 1);
1304 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e1, 1, Rott2, 1);
1305 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e0, 1, Rott1, 1, Rott1, 1);
1306 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e1, 1, Rott2, 1, Rott2, 1);
1307
1308 // Multiply H and \partial \phi / \partial t which is assumed to be u_{\phi}
1309 Vmath::Vmul(nq, &h[0], 1, &Rott1[0], 1, &Rott1[0], 1);
1310 Vmath::Vmul(nq, &h[0], 1, &Rott2[0], 1, &Rott2[0], 1);
1311
1312 Vmath::Neg(nq, Rott1, 1);
1313 Vmath::Neg(nq, Rott2, 1);
1314
1315 Array<OneD, NekDouble> tmpc1(ncoeffs);
1316 Array<OneD, NekDouble> tmpc2(ncoeffs);
1317 m_fields[0]->IProductWRTBase(Rott1, tmpc1);
1318 m_fields[0]->IProductWRTBase(Rott2, tmpc2);
1319
1320 Vmath::Vadd(ncoeffs, tmpc1, 1, outarray[1], 1, outarray[1], 1);
1321 Vmath::Vadd(ncoeffs, tmpc2, 1, outarray[2], 1, outarray[2], 1);
1322}
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:1328
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.hpp:366
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 962 of file MMFSWE.cpp.
966{
967 NekDouble MageF1, MageF2, MageB1, MageB2;
968 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
969 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
970
972 NekDouble uRF, vRF, uLB, vLB;
973
974 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
975 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
976
977 // uRF = uR component in moving frames e^{Fwd}
978 // vRF = vR component in moving frames e^{Fwd}
979 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
980 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
981
982 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
983 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
984 numfluxF[2] = 0.0;
985
986 numfluxF[3] =
987 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
988 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
989 numfluxF[5] = 0.0;
990
991 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
992 numfluxF[7] =
993 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
994 numfluxF[8] = 0.0;
995
996 // uLB = uL component in moving frames e^{Bwd}
997 // vLB = vL component in moving frames e^{Bwd}
998 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
999 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1000
1001 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1002 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1003 numfluxB[2] = 0.0;
1004
1005 numfluxB[3] =
1006 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1007 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1008 numfluxB[5] = 0.0;
1009
1010 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1011 numfluxB[7] =
1012 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1013 numfluxB[8] = 0.0;
1014}
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:1169
References ComputeMagAndDot(), and m_g.
Referenced by NumericalSWEFlux().
◆ Checkpoint_Output_Cartesian()
|
private |
Definition at line 2768 of file MMFSWE.cpp.
2769{
2772
2773 NekDouble rad_earth = 6.37122 * 1000000;
2774
2775 int nvariables = 7;
2776
2777 // vector u in Cartesian coordinates
2778 std::vector<std::string> variables(nvariables);
2779
2780 variables[0] = "eta";
2781 variables[1] = "hstar";
2782 variables[2] = "vorticity";
2783 variables[3] = "ux";
2784 variables[4] = "uy";
2785 variables[5] = "uz";
2786 variables[6] = "null";
2787
2788 // Obtain \vec{u} in cartesian coordinate
2789 Array<OneD, Array<OneD, NekDouble>> fieldphys(nvariables);
2790 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvariables);
2791 for (int i = 0; i < nvariables; ++i)
2792 {
2793 fieldphys[i] = Array<OneD, NekDouble>(nq, 0.0);
2794 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2795 }
2796
2798 &fieldphys[0][0], 1);
2799 // Vmath::Vadd(nq, &(m_fields[0]->GetPhys())[0], 1, &m_depth[0], 1,
2800 // &fieldphys[1][0], 1);
2801
2803 Vmath::Neg(nq, &fieldphys[1][0], 1);
2805 Vmath::Smul(nq, rad_earth, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2806
2807 Array<OneD, NekDouble> utmp(nq);
2808 Array<OneD, NekDouble> vtmp(nq);
2809
2812
2813 ComputeVorticity(utmp, vtmp, fieldphys[2]);
2814 // u_x = u e^1_x + v e^2_x
2815 int indx = 3;
2817 {
2819 &fieldphys[k + indx][0], 1);
2821 &fieldphys[k + indx][0], 1, &fieldphys[k + indx][0], 1);
2822 }
2823
2824 for (int i = 0; i < nvariables; ++i)
2825 {
2826 m_fields[0]->FwdTrans(fieldphys[i], fieldcoeffs[i]);
2827 }
2828
2830}
void ComputeVorticity(const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
Definition: MMFSWE.cpp:2214
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1315
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:183
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.hpp:194
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
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 1328 of file MMFSWE.cpp.
1332{
1333 int j, k;
1335
1336 Array<OneD, NekDouble> tmp(nq, 0.0);
1337 Array<OneD, NekDouble> tmpderiv(nq);
1338
1339 outarray = Array<OneD, NekDouble>(nq, 0.0);
1341 {
1343 {
1344 // Compute d e^m / d \xi_1 and d e^m / d \xi_2
1347
1348 Vmath::Vmul(nq, &physarray[j + 1][0], 1, &tmpderiv[0], 1,
1349 &tmpderiv[0], 1);
1350
1352 &outarray[0], 1, &outarray[0], 1);
1353 }
1354 }
1355}
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 2154 of file MMFSWE.cpp.
2157{
2159
2160 Array<OneD, NekDouble> tmp(nq);
2161 Array<OneD, NekDouble> htmp(nq);
2162 Array<OneD, NekDouble> hstmp(nq);
2163
2164 Vmath::Vmul(nq, u, 1, u, 1, tmp, 1);
2165 Vmath::Vvtvp(nq, v, 1, v, 1, tmp, 1, tmp, 1);
2166 Vmath::Vmul(nq, eta, 1, tmp, 1, tmp, 1);
2167
2169 Vmath::Vmul(nq, htmp, 1, htmp, 1, htmp, 1);
2170
2172 Vmath::Vmul(nq, hstmp, 1, hstmp, 1, hstmp, 1);
2173
2174 Vmath::Vsub(nq, htmp, 1, hstmp, 1, htmp, 1);
2176
2177 Vmath::Vadd(nq, htmp, 1, tmp, 1, tmp, 1);
2178 Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2179
2181}
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.hpp:220
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 2183 of file MMFSWE.cpp.
2186{
2188
2189 Array<OneD, NekDouble> hstartmp(nq);
2190 Array<OneD, NekDouble> tmp(nq);
2191
2193
2194 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2195 Array<OneD, Array<OneD, NekDouble>> Curlu(m_spacedim);
2197 {
2198 uvec[i] = Array<OneD, NekDouble>(nq);
2199 Curlu[i] = Array<OneD, NekDouble>(nq);
2200 }
2201
2202 ComputeVorticity(u, v, tmp);
2203
2205 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
2206 Vmath::Vdiv(nq, tmp, 1, hstartmp, 1, tmp, 1);
2207 Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2208
2210}
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.hpp:126
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 884 of file MMFSWE.cpp.
888{
890
891 NekDouble hC, huC, hvC, SL, SR, Sstar;
894
895 // Compute SL
896 if (hstar > hL)
897 {
898 SL = uL - cL * sqrt(0.5 * ((hstar * hstar + hstar * hL) / (hL * hL)));
899 }
900 else
901 {
902 SL = uL - cL;
903 }
904
905 // Compute SR
906 if (hstar > hR)
907 {
908 SR = uR + cR * sqrt(0.5 * ((hstar * hstar + hstar * hR) / (hR * hR)));
909 }
910 else
911 {
912 SR = uR + cR;
913 }
914
915 if (fabs(hR * (uR - SR) - hL * (uL - SL)) <= 1.0e-15)
916 {
917 Sstar = 0.0;
918 }
919 else
920 {
921 Sstar = (SL * hR * (uR - SR) - SR * hL * (uL - SL)) /
922 (hR * (uR - SR) - hL * (uL - SL));
923 }
924
925 if (SL >= 0)
926 {
927 hflux = hL * uL;
928 huflux = uL * uL * hL + 0.5 * g * hL * hL;
929 hvflux = hL * uL * vL;
930 }
931 else if (SR <= 0)
932 {
933 hflux = hR * uR;
934 huflux = uR * uR * hR + 0.5 * g * hR * hR;
935 hvflux = hR * uR * vR;
936 }
937 else
938 {
939 if ((SL < 0) && (Sstar >= 0))
940 {
941 hC = hL * ((SL - uL) / (SL - Sstar));
942 huC = hC * Sstar;
943 hvC = hC * vL;
944
945 hflux = hL * uL + SL * (hC - hL);
946 huflux = (uL * uL * hL + 0.5 * g * hL * hL) + SL * (huC - hL * uL);
947 hvflux = (uL * vL * hL) + SL * (hvC - hL * vL);
948 }
949 else
950 {
951 hC = hR * ((SR - uR) / (SR - Sstar));
952 huC = hC * Sstar;
953 hvC = hC * vR;
954
955 hflux = hR * uR + SR * (hC - hR);
956 huflux = (uR * uR * hR + 0.5 * g * hR * hR) + SR * (huC - hR * uR);
957 hvflux = (uR * vR * hR) + SR * (hvC - hR * vR);
958 }
959 }
960}
References m_g, and tinysimd::sqrt().
Referenced by RiemannSolverHLLC().
◆ ComputeMagAndDot()
|
protected |
Definition at line 1169 of file MMFSWE.cpp.
1174{
1175 NekDouble MF1x, MF1y, MF1z, MF2x, MF2y, MF2z;
1176 NekDouble MB1x, MB1y, MB1z, MB2x, MB2y, MB2z;
1177
1178 MF1x = m_MFtraceFwd[0][0][index];
1179 MF1y = m_MFtraceFwd[0][1][index];
1180 MF1z = m_MFtraceFwd[0][2][index];
1181
1182 MF2x = m_MFtraceFwd[1][0][index];
1183 MF2y = m_MFtraceFwd[1][1][index];
1184 MF2z = m_MFtraceFwd[1][2][index];
1185
1186 MB1x = m_MFtraceBwd[0][0][index];
1187 MB1y = m_MFtraceBwd[0][1][index];
1188 MB1z = m_MFtraceBwd[0][2][index];
1189
1190 MB2x = m_MFtraceBwd[1][0][index];
1191 MB2y = m_MFtraceBwd[1][1][index];
1192 MB2z = m_MFtraceBwd[1][2][index];
1193
1194 // MFtrace = MFtrace [ j*spacedim + k ], j = shape, k = sapce
1195 MageF1 = MF1x * MF1x + MF1y * MF1y + MF1z * MF1z;
1196 MageF2 = MF2x * MF2x + MF2y * MF2y + MF2z * MF2z;
1197 MageB1 = MB1x * MB1x + MB1y * MB1y + MB1z * MB1z;
1198 MageB2 = MB2x * MB2x + MB2y * MB2y + MB2z * MB2z;
1199
1200 eF1_cdot_eB1 = MF1x * MB1x + MF1y * MB1y + MF1z * MB1z;
1201 eF1_cdot_eB2 = MF1x * MB2x + MF1y * MB2y + MF1z * MB2z;
1202 eF2_cdot_eB1 = MF2x * MB1x + MF2y * MB1y + MF2z * MB1z;
1203 eF2_cdot_eB2 = MF2x * MB2x + MF2y * MB2y + MF2z * MB2z;
1204}
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
Definition: MMFSystem.h:196
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
Definition: MMFSystem.h:197
References Nektar::SolverUtils::MMFSystem::m_MFtraceBwd, and Nektar::SolverUtils::MMFSystem::m_MFtraceFwd.
Referenced by AverageFlux(), LaxFriedrichFlux(), and RusanovFlux().
◆ ComputeMass()
Definition at line 2144 of file MMFSWE.cpp.
2145{
2147
2148 Array<OneD, NekDouble> tmp(nq);
2150
2152}
References m_depth, Nektar::SolverUtils::EquationSystem::m_fields, and Vmath::Vadd().
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ ComputeNablaCdotVelocity()
Definition at line 2235 of file MMFSWE.cpp.
2236{
2238
2239 Array<OneD, NekDouble> velcoeff(nq, 0.0);
2240
2241 Array<OneD, NekDouble> Dtmp0(nq);
2242 Array<OneD, NekDouble> Dtmp1(nq);
2243 Array<OneD, NekDouble> Dtmp2(nq);
2244 Array<OneD, NekDouble> Drv(nq);
2245
2246 vellc = Array<OneD, NekDouble>(nq, 0.0);
2247
2248 // m_vellc = \nabla m_vel \cdot tan_i
2249 Array<OneD, NekDouble> tmp(nq);
2250 Array<OneD, NekDouble> vessel(nq);
2251
2253 {
2254 Vmath::Zero(nq, velcoeff, 1);
2256 {
2257 // a_j = tan_j cdot m_vel
2259 1, &velcoeff[0], 1, &velcoeff[0], 1);
2260 }
2261
2262 // d a_j / d x^k
2263 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
2264
2266 {
2267 // tan_j^k ( d a_j / d x^k )
2268 switch (k)
2269 {
2270 case (0):
2271 {
2273 1, &vellc[0], 1, &vellc[0], 1);
2274 }
2275 break;
2276
2277 case (1):
2278 {
2280 1, &vellc[0], 1, &vellc[0], 1);
2281 }
2282 break;
2283
2284 case (2):
2285 {
2287 1, &vellc[0], 1, &vellc[0], 1);
2288 }
2289 break;
2290 }
2291 }
2292 }
2293}
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 2738 of file MMFSWE.cpp.
2739{
2740 NekDouble uphi, f, dh;
2741
2742 uphi = ComputeUnstableJetuphi(theta);
2743 f = 2.0 * m_Omega * sin(theta);
2744
2745 dh = f * uphi + tan(theta) * uphi * uphi;
2746
2747 return dh;
2748}
NekDouble ComputeUnstableJetuphi(const NekDouble theta)
Definition: MMFSWE.cpp:2750
References ComputeUnstableJetuphi(), and m_Omega.
Referenced by UnstableJetFlow().
◆ ComputeUnstableJetuphi()
Definition at line 2750 of file MMFSWE.cpp.
2751{
2752 NekDouble uphi;
2753
2755 {
2758 }
2759
2760 else
2761 {
2762 uphi = 0.0;
2763 }
2764
2765 return uphi;
2766}
References m_en, m_theta0, m_theta1, and m_uthetamax.
Referenced by ComputeUnstableJetEta(), and UnstableJetFlow().
◆ ComputeVorticity()
|
protected |
Definition at line 2214 of file MMFSWE.cpp.
2217{
2219
2220 Array<OneD, NekDouble> tmp(nq);
2221
2222 Vorticity = Array<OneD, NekDouble>(nq, 0.0);
2223
2226 &Vorticity[0], 1);
2227
2229 Vmath::Neg(nq, tmp, 1);
2231
2232 Vmath::Vadd(nq, tmp, 1, Vorticity, 1, Vorticity, 1);
2233}
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Definition: MMFSystem.h:193
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 2913 of file MMFSWE.cpp.
2914{
2916
2917 // u = hu/h
2919 m_fields[1]->UpdatePhys(), 1);
2920
2921 // v = hv/ v
2923 m_fields[2]->UpdatePhys(), 1);
2924
2925 // \eta = h - d
2927 m_fields[0]->UpdatePhys(), 1);
2928}
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:341
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 2833 of file MMFSWE.cpp.
2836{
2838
2839 if (physin[0].data() == physout[0].data())
2840 {
2841 // copy indata and work with tmp array
2842 Array<OneD, Array<OneD, NekDouble>> tmp(3);
2843 for (int i = 0; i < 3; ++i)
2844 {
2845 // deep copy
2846 tmp[i] = Array<OneD, NekDouble>(nq);
2847 Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2848 }
2849
2850 // \eta = h - d
2852
2853 // u = hu/h
2854 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
2855
2856 // v = hv/h
2857 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
2858 }
2859 else
2860 {
2861 // \eta = h - d
2863
2864 // u = hu/h
2865 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
2866
2867 // v = hv/h
2868 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
2869 }
2870}
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 66 of file MMFSWE.h.
69 {
71 MemoryManager<MMFSWE>::AllocateSharedPtr(pSession, pGraph);
72 p->InitObject();
74 }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:166
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition: EquationSystem.h:65
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 1364 of file MMFSWE.cpp.
1367{
1369 {
1371 {
1372 ConservativeToPrimitive(inarray, outarray);
1373 SetBoundaryConditions(outarray, time);
1374 PrimitiveToConservative(outarray, outarray);
1375 }
1376 break;
1377 default:
1379 break;
1380 }
1381}
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &inarray, NekDouble time)
Definition: MMFSWE.cpp:1383
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:526
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 419 of file MMFSWE.cpp.
422{
423 int i;
424 int nvariables = inarray.size();
427
428 // inarray in physical space
429 Array<OneD, Array<OneD, NekDouble>> physarray(nvariables);
430 Array<OneD, Array<OneD, NekDouble>> modarray(nvariables);
431
432 for (i = 0; i < nvariables; ++i)
433 {
434 physarray[i] = Array<OneD, NekDouble>(nq);
435 modarray[i] = Array<OneD, NekDouble>(ncoeffs);
436 }
437
438 // (h, hu, hv) -> (\eta, u, v)
439 ConservativeToPrimitive(inarray, physarray);
440
441 // Weak Directional Derivative
442 WeakDGSWEDirDeriv(physarray, modarray);
443 AddDivForGradient(physarray, modarray);
444
445 for (i = 0; i < nvariables; ++i)
446 {
447 Vmath::Neg(ncoeffs, modarray[i], 1);
448 }
449
450 // coriolis forcing
452 {
453 AddCoriolis(physarray, modarray);
454 }
455
456 // Bottom Elevation Effect
457 // Add g H \nabla m_depth
458 AddElevationEffect(physarray, modarray);
459
460 // Add terms concerning the rotation of the moving frame
462 {
463 AddRotation(physarray, modarray);
464 }
465
466 for (i = 0; i < nvariables; ++i)
467 {
468 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
469 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
470 }
471}
void AddDivForGradient(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:543
void AddCoriolis(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1207
void WeakDGSWEDirDeriv(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField)
Definition: MMFSWE.cpp:473
void AddRotation(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1281
void AddElevationEffect(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1252
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:283
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 1703 of file MMFSWE.cpp.
1704{
1706 {
1708 {
1710 }
1711 break;
1712
1714 {
1716 }
1717 break;
1718
1723 {
1725 }
1726 break;
1727
1728 default:
1729 break;
1730 }
1731}
void EvaluateStandardCoriolis(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1765
void EvaluateCoriolisForZonalFlow(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1733
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:741
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 1733 of file MMFSWE.cpp.
1734{
1736 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1737
1738 Array<OneD, NekDouble> x(nq);
1739 Array<OneD, NekDouble> y(nq);
1740 Array<OneD, NekDouble> z(nq);
1741
1743
1744 NekDouble x0j, x1j, x2j;
1745 NekDouble tmp;
1746
1747 outarray = Array<OneD, NekDouble>(nq, 0.0);
1748 for (int j = 0; j < nq; ++j)
1749 {
1750 x0j = x[j];
1751 x1j = y[j];
1752 x2j = z[j];
1753
1754 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1755 cos_theta);
1756
1757 // H = 2 \Omega *(- \cos \phi \cos \theta \sin \alpha + \sin \theta \cos
1758 // \alpha )
1759 tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
1760 sin_theta * cos(m_alpha);
1761 outarray[j] = 2.0 * m_Omega * tmp;
1762 }
1763}
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:792
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 1765 of file MMFSWE.cpp.
1766{
1768
1769 NekDouble x0j, x1j, x2j;
1770 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1771
1772 Array<OneD, NekDouble> x(nq);
1773 Array<OneD, NekDouble> y(nq);
1774 Array<OneD, NekDouble> z(nq);
1775
1777
1778 outarray = Array<OneD, NekDouble>(nq, 0.0);
1779 for (int j = 0; j < nq; ++j)
1780 {
1781 x0j = x[j];
1782 x1j = y[j];
1783 x2j = z[j];
1784
1785 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1786 cos_theta);
1787
1788 outarray[j] = 2.0 * m_Omega * sin_theta;
1789 }
1790}
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 1563 of file MMFSWE.cpp.
1564{
1566
1567 m_depth = Array<OneD, NekDouble>(nq, 0.0);
1568
1570 {
1572 {
1574 }
1575 break;
1576
1578 {
1579 // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1580 Array<OneD, NekDouble> x(nq);
1581 Array<OneD, NekDouble> y(nq);
1582 Array<OneD, NekDouble> z(nq);
1583
1585
1586 NekDouble x0j, x1j, x2j;
1587 NekDouble sin_varphi, cos_varphi, sin_theta, cos_theta;
1588
1589 for (int j = 0; j < nq; ++j)
1590 {
1591 x0j = x[j];
1592 x1j = y[j];
1593 x2j = z[j];
1594
1595 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1596 sin_theta, cos_theta);
1597 m_depth[j] =
1600 }
1601 }
1602 break;
1603
1605 {
1606 // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1607 Array<OneD, NekDouble> x(nq);
1608 Array<OneD, NekDouble> y(nq);
1609 Array<OneD, NekDouble> z(nq);
1610
1612
1613 int indx = 0;
1614 NekDouble x0j, x1j, x2j, dist2;
1615 NekDouble phi, theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
1616 NekDouble hRad, phic, thetac;
1617 NekDouble Tol = 0.000001;
1618
1619 hRad = m_pi / 9.0;
1620 phic = -m_pi / 2.0;
1621 thetac = m_pi / 6.0;
1622
1623 for (int j = 0; j < nq; ++j)
1624 {
1625 x0j = x[j];
1626 x1j = y[j];
1627 x2j = z[j];
1628
1629 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1630 sin_theta, cos_theta);
1631
1633 std::abs(sin(thetac) - sin_theta)) < Tol)
1634 {
1635 std::cout << "A point " << j
1636 << " is coincient with the singularity "
1637 << std::endl;
1638 indx = 1;
1639 }
1640
1641 phi = atan2(sin_varphi, cos_varphi);
1642 theta = atan2(sin_theta, cos_theta);
1643
1644 // Compute r
1645 dist2 = (phi - phic) * (phi - phic) +
1646 (theta - thetac) * (theta - thetac);
1647
1648 if (dist2 > hRad * hRad)
1649 {
1650 dist2 = hRad * hRad;
1651 }
1652
1654 }
1655
1656 if (!indx)
1657 {
1658 std::cout << "No point is coincident with the singularity point"
1659 << std::endl;
1660 }
1661 }
1662 break;
1663
1665 {
1666 for (int j = 0; j < nq; ++j)
1667 {
1669 }
1670 }
1671 break;
1672
1675 {
1677 }
1678 break;
1679
1680 default:
1681 {
1683 }
1684 break;
1685 }
1686
1687 // Comptue \nabla m_depth \cdot e^i
1689
1691 {
1692 m_Derivdepth[j] = Array<OneD, NekDouble>(nq);
1694 m_Derivdepth[j]);
1695 }
1696
1697 std::cout << "Water Depth (m_depth) was generated with mag = "
1701}
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2307
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
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.hpp:685
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().
◆ GetSWEFluxVector()
|
protected |
Definition at line 573 of file MMFSWE.cpp.
576{
578
579 switch (i)
580 {
581 // flux function for the h equation = [(\eta + d ) u, (\eta + d) v ]
582 case 0:
583 {
584 // h in flux 1
586
587 // hu in flux 0
588 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
589
590 // hv in flux 1
591 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
592 }
593 break;
594
595 // flux function for the hu equation = [Hu^2 + 0.5 g H^2, Huv]
596 case 1:
597 {
598 Array<OneD, NekDouble> tmp(nq);
599
600 // h in tmp
602
603 // hu in flux 1
604 Vmath::Vmul(nq, tmp, 1, physfield[1], 1, flux[1], 1);
605
606 // huu in flux 0
607 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
608
609 // hh in tmp
610 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
611
612 // huu + 0.5 g hh in flux 0
613 // Daxpy overwrites flux[0] on exit
615
616 // huv in flux 1
617 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
618 }
619 break;
620
621 // flux function for the hv equation = [Huv, Hv^2]
622 case 2:
623 {
624 Array<OneD, NekDouble> tmp(nq);
625
626 // h in tmp
628
629 // hv in flux 0
630 Vmath::Vmul(nq, tmp, 1, physfield[2], 1, flux[0], 1);
631
632 // hvv in flux 1
633 Vmath::Vmul(nq, flux[0], 1, physfield[2], 1, flux[1], 1);
634
635 // huv in flux 0
636 Vmath::Vmul(nq, flux[0], 1, physfield[1], 1, flux[0], 1);
637
638 // hh in tmp
639 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
640
641 // hvv + 0.5 g hh in flux 1
643 }
644 break;
645
646 // flux function 0.5 g h * h
647 case 3:
648 {
649 Array<OneD, NekDouble> h(nq);
650 flux[0] = Array<OneD, NekDouble>(nq, 0.0);
651
652 // h in tmp
654
655 // hh in tmp
656 Vmath::Vmul(nq, h, 1, h, 1, h, 1);
657
658 // 0.5 g hh in flux 0
660 }
661 break;
662
663 default:
664 break;
665 }
666}
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:135
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 2396 of file MMFSWE.cpp.
2399{
2401
2402 NekDouble uhat, vhat;
2403 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2404 NekDouble x0j, x1j, x2j;
2405
2406 Array<OneD, NekDouble> eta(nq, 0.0);
2407 Array<OneD, NekDouble> u(nq, 0.0);
2408 Array<OneD, NekDouble> v(nq, 0.0);
2409
2410 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2412 {
2413 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2414 }
2415
2416 Array<OneD, NekDouble> x(nq);
2417 Array<OneD, NekDouble> y(nq);
2418 Array<OneD, NekDouble> z(nq);
2419
2421
2422 for (int j = 0; j < nq; ++j)
2423 {
2424 x0j = x[j];
2425 x1j = y[j];
2426 x2j = z[j];
2427
2428 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2429 cos_theta);
2430
2431 // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2433 sin_theta * sin_theta;
2434
2435 uhat = m_u0 * cos_theta;
2436 vhat = 0.0;
2437
2438 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2439 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2440 uvec[2][j] = vhat * cos_theta;
2441 }
2442
2443 // Projection of u onto the tangent plane with conserving the mag. of the
2444 // velocity.
2445 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2446
2448 {
2449 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2450 }
2451
2452 // u is projected on the tangent plane with preserving its length
2453 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2454
2455 // Change it to the coordinate of moving frames
2456 // CartesianToMovingframes(0,uvecproj,u);
2457 // CartesianToMovingframes(1,uvecproj,v);
2458
2459 u = CartesianToMovingframes(uvec, 0);
2460 v = CartesianToMovingframes(uvec, 1);
2461
2463 {
2464 case (0):
2465 {
2466 outfield = eta;
2467 }
2468 break;
2469
2470 case (1):
2471 {
2472 outfield = u;
2473 }
2474 break;
2475
2476 case (2):
2477 {
2478 outfield = v;
2479 }
2480 break;
2481 }
2482}
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble > > &uvec, unsigned int field)
Definition: MMFSystem.cpp:771
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), FilterPython_Function::field, 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 1016 of file MMFSWE.cpp.
1020{
1021 int nvariables = 3;
1022 NekDouble MageF1, MageF2, MageB1, MageB2;
1023 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1024 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1025
1027 NekDouble uRF, vRF, uLB, vLB;
1028 NekDouble velL, velR, lambdaF, lambdaB;
1029
1030 Array<OneD, NekDouble> EigF(nvariables);
1031 Array<OneD, NekDouble> EigB(nvariables);
1032
1033 // Compute Magnitude and Dot product of moving frames for the index
1034 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1035 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1036
1037 // Get the velocity in the normal to the edge
1040
1041 EigF[0] = velL - sqrt(g * hL);
1042 EigF[1] = velL;
1043 EigF[2] = velL + sqrt(g * hL);
1044
1045 EigB[0] = velR - sqrt(g * hR);
1046 EigB[1] = velR;
1047 EigB[2] = velR + sqrt(g * hR);
1048
1049 lambdaF = Vmath::Vamax(nvariables, EigF, 1);
1050 lambdaB = Vmath::Vamax(nvariables, EigB, 1);
1051
1052 // uRF = uR component in moving frames e^{Fwd}
1053 // vRF = vR component in moving frames e^{Fwd}
1054 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1055 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1056
1057 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1058 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1059 numfluxF[2] = 0.5 * lambdaF * (hL - hR);
1060
1061 numfluxF[3] = 0.5 * (hL * uL * uL * MageF1 + hR * uRF * uRF * MageB1 +
1062 0.5 * g * (hL * hL + hR * hR));
1063 numfluxF[4] = 0.5 * (hL * uL * vL * MageF1 + hR * uRF * vRF * MageB1);
1064 numfluxF[5] = 0.5 * lambdaF * (uL * hL - uRF * hR);
1065
1066 numfluxF[6] = 0.5 * (hL * uL * vL * MageF2 + hR * uRF * vRF * MageB2);
1067 numfluxF[7] = 0.5 * (hL * vL * vL * MageF2 + hR * vRF * vRF * MageB2 +
1068 0.5 * g * (hL * hL + hR * hR));
1069 numfluxF[8] = 0.5 * lambdaF * (vL * hL - vRF * hR);
1070
1071 // uLB = uL component in moving frames e^{Bwd}
1072 // vLB = vL component in moving frames e^{Bwd}
1073 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1074 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1075
1076 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1077 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1078 numfluxB[2] = 0.5 * lambdaB * (hL - hR);
1079
1080 numfluxB[3] = 0.5 * (hR * uR * uR * MageB1 + hR * uLB * uLB * MageF1 +
1081 0.5 * g * (hR * hR + hL * hL));
1082 numfluxB[4] = 0.5 * (hR * uR * vR * MageB1 + hR * uLB * vLB * MageF1);
1083 numfluxB[5] = 0.5 * lambdaB * (uLB * hL - uR * hR);
1084
1085 numfluxB[6] = 0.5 * (hR * uR * vR * MageB2 + hR * uLB * vLB * MageF2);
1086 numfluxB[7] = 0.5 * (hR * vR * vR * MageB2 + hR * vLB * vLB * MageF2 +
1087 0.5 * g * (hR * hR + hL * hL));
1088 numfluxB[8] = 0.5 * lambdaB * (vLB * hL - vR * hR);
1089}
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:187
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:186
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 668 of file MMFSWE.cpp.
671{
672
673 int i, k;
675 int nvariables = 3; // only the dependent variables
676
677 // get temporary arrays
678 Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
679 Array<OneD, Array<OneD, NekDouble>> Bwd(nvariables);
680 Array<OneD, NekDouble> DepthFwd(nTraceNumPoints);
681 Array<OneD, NekDouble> DepthBwd(nTraceNumPoints);
682
683 for (i = 0; i < nvariables; ++i)
684 {
685 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
686 Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
687 }
688
689 // get the physical values at the trace
690 for (i = 0; i < nvariables; ++i)
691 {
692 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
693
694 // Copy Fwd to Bwd at boundaries
696 }
699
700 // note that we are using the same depth - i.e. the depth is assumed
701 // continuous...
703 {
705 {
706 // rotate the values to the normal direction
707 NekDouble tmpX, tmpY;
708
709 // Fwd[1] = hu^+ = ( h^1^+ (e^1 \cdot n) + h^2^+ ( e^2 \cdot n) )
710 // Fwd[2] = hv^+ = ( h^1^+ (e^1 \cdot n^{\perp}) + h^2^+ ( e^2 \cdot
711 // n^{\perp}) )
712 // Bwd[1] = hu^- = ( h^1^- (e^1 \cdot n) + h^2^- ( e^2 \cdot n) )
713 // Bwd[2] = hv^- = ( h^1^- (e^1 \cdot n^{\perp}) + h^2^- ( e^2 \cdot
714 // n^{\perp}) )
715
716 for (k = 0; k < nTraceNumPoints; ++k)
717 {
718 tmpX = Fwd[1][k] * m_ncdotMFFwd[0][k] +
719 Fwd[2][k] * m_ncdotMFFwd[1][k];
720 tmpY = Fwd[1][k] * m_nperpcdotMFFwd[0][k] +
721 Fwd[2][k] * m_nperpcdotMFFwd[1][k];
722 Fwd[1][k] = tmpX;
723 Fwd[2][k] = tmpY;
724
725 tmpX = Bwd[1][k] * m_ncdotMFBwd[0][k] +
726 Bwd[2][k] * m_ncdotMFBwd[1][k];
727 tmpY = Bwd[1][k] * m_nperpcdotMFBwd[0][k] +
728 Bwd[2][k] * m_nperpcdotMFBwd[1][k];
729 Bwd[1][k] = tmpX;
730 Bwd[2][k] = tmpY;
731 }
732
733 // Solve the Riemann problem
734 NekDouble denomFwd, denomBwd;
735 Array<OneD, NekDouble> numfluxF(nvariables);
736 Array<OneD, NekDouble> numfluxB(nvariables);
737
738 NekDouble eF1n, eF2n, eB1n, eB2n;
739 NekDouble eF1t, eF2t, eB1t, eB2t;
740 for (k = 0; k < nTraceNumPoints; ++k)
741 {
742 RiemannSolverHLLC(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
743 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
744 Bwd[2][k], numfluxF, numfluxB);
745
746 // uflux = 1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
747 // \cdot n)(e^1 \cdot n^{\perp} ) ] ( e^2 \cdot n^{\perp} huflux
748 // - e^2 \cdot n hvflux
749 // vflux = -1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
750 // \cdot n)(e^1 \cdot n^{\perp} ) ] ( -e^1 \cdot n^{\perp}
751 // huflux + e^1 \cdot n hvflux
752 eF1n = m_ncdotMFFwd[0][k];
753 eF2n = m_ncdotMFFwd[1][k];
754 eB1n = m_ncdotMFBwd[0][k];
755 eB2n = m_ncdotMFBwd[1][k];
756
757 eF1t = m_nperpcdotMFFwd[0][k];
758 eF2t = m_nperpcdotMFFwd[1][k];
759 eB1t = m_nperpcdotMFBwd[0][k];
760 eB2t = m_nperpcdotMFBwd[1][k];
761
762 denomFwd = eF1n * eF2t - eF2n * eF1t;
763 denomBwd = eB1n * eB2t - eB2n * eB1t;
764
765 numfluxFwd[0][k] = numfluxF[0];
766 numfluxFwd[1][k] = (1.0 / denomFwd) *
767 (eF2t * numfluxF[1] - eF2n * numfluxF[2]);
768 numfluxFwd[2][k] =
769 (1.0 / denomFwd) *
770 (-1.0 * eF1t * numfluxF[1] + eF1n * numfluxF[2]);
771
772 numfluxBwd[0][k] = 1.0 * numfluxB[0];
773 numfluxBwd[1][k] = (1.0 / denomBwd) *
774 (eB2t * numfluxB[1] - eB2n * numfluxB[2]);
775 numfluxBwd[2][k] =
776 (1.0 / denomBwd) *
777 (-1.0 * eB1t * numfluxB[1] + eB1n * numfluxB[2]);
778 }
779 }
780 break;
781
785 {
786 Array<OneD, NekDouble> numfluxF(nvariables * (m_shapedim + 1));
787 Array<OneD, NekDouble> numfluxB(nvariables * (m_shapedim + 1));
788
789 for (k = 0; k < nTraceNumPoints; ++k)
790 {
792 {
793 AverageFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
794 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
795 Bwd[2][k], numfluxF, numfluxB);
796 }
797
799 {
800 LaxFriedrichFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
801 Fwd[2][k], Bwd[0][k] + DepthFwd[k],
802 Bwd[1][k], Bwd[2][k], numfluxF, numfluxB);
803 }
804
806 {
807 RusanovFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
808 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
809 Bwd[2][k], numfluxF, numfluxB);
810 }
811
812 int indx;
813 NekDouble tmpF0, tmpF1, tmpB0, tmpB1;
814 for (i = 0; i < nvariables; ++i)
815 {
816 indx = i * (m_shapedim + 1);
817
818 tmpF0 = numfluxF[indx] * m_ncdotMFFwd[0][k];
819 tmpF1 = numfluxF[indx + 1] * m_ncdotMFFwd[1][k];
820
821 tmpB0 = numfluxB[indx] * m_ncdotMFBwd[0][k];
822 tmpB1 = numfluxB[indx + 1] * m_ncdotMFBwd[1][k];
823
824 numfluxFwd[i][k] = tmpF0 + tmpF1 + numfluxF[indx + 2];
825 numfluxBwd[i][k] = tmpB0 + tmpB1 + numfluxB[indx + 2];
826 }
827 }
828 }
829 break;
830
831 default:
832 {
834 }
835 break;
836 }
837}
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:962
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:839
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:1016
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:1091
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:316
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
Definition: MMFSystem.h:189
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFBwd
Definition: MMFSystem.h:190
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:835
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 2930 of file MMFSWE.cpp.
2931{
2933
2934 // h = \eta + d
2936 m_fields[0]->UpdatePhys(), 1);
2937
2938 // hu = h * u
2940 m_fields[1]->UpdatePhys(), 1);
2941
2942 // hv = h * v
2944 m_fields[2]->UpdatePhys(), 1);
2945}
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 2872 of file MMFSWE.cpp.
2875{
2876
2878
2879 if (physin[0].data() == physout[0].data())
2880 {
2881 // copy indata and work with tmp array
2882 Array<OneD, Array<OneD, NekDouble>> tmp(3);
2883 for (int i = 0; i < 3; ++i)
2884 {
2885 // deep copy
2886 tmp[i] = Array<OneD, NekDouble>(nq);
2887 Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2888 }
2889
2890 // h = \eta + d
2892
2893 // hu = h * u
2894 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
2895
2896 // hv = h * v
2897 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
2898 }
2899 else
2900 {
2901 // h = \eta + d
2903
2904 // hu = h * u
2905 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
2906
2907 // hv = h * v
2908 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
2909 }
2910}
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vmul().
◆ RiemannSolverHLLC()
|
protected |
Definition at line 839 of file MMFSWE.cpp.
844{
846
849
850 NekDouble uRF, vRF, uLB, vLB;
851 NekDouble hstarF, hstarB;
852
853 // Temporary assignments
854 uRF = uR;
855 vRF = vR;
856 uLB = uL;
857 vLB = vL;
858
859 // the two-rarefaction wave assumption
860 hstarF = 0.5 * (cL + cR) + 0.25 * (uL - uRF);
861 hstarF *= hstarF;
862 hstarF *= (1.0 / g);
863
864 hstarB = 0.5 * (cL + cR) + 0.25 * (uLB - uR);
865 hstarB *= hstarB;
866 hstarB *= (1.0 / g);
867
868 NekDouble hfluxF, hufluxF, hvfluxF;
869 NekDouble hfluxB, hufluxB, hvfluxB;
870 Computehhuhvflux(hL, uL, vL, hR, uRF, vRF, hstarF, hfluxF, hufluxF,
871 hvfluxF);
872 Computehhuhvflux(hL, uLB, vLB, hR, uR, vR, hstarB, hfluxB, hufluxB,
873 hvfluxB);
874
875 numfluxF[0] = hfluxF;
876 numfluxF[1] = hufluxF;
877 numfluxF[2] = hvfluxF;
878
879 numfluxB[0] = hfluxB;
880 numfluxB[1] = hufluxB;
881 numfluxB[2] = hvfluxB;
882}
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:884
References Computehhuhvflux(), m_g, and tinysimd::sqrt().
Referenced by NumericalSWEFlux().
◆ RossbyWave()
|
protected |
Definition at line 2593 of file MMFSWE.cpp.
2594{
2596 NekDouble uhat, vhat;
2597 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2598 NekDouble x0j, x1j, x2j;
2599 NekDouble Ath, Bth, Cth, tmp;
2600
2601 Array<OneD, NekDouble> eta(nq, 0.0);
2602 Array<OneD, NekDouble> u(nq, 0.0);
2603 Array<OneD, NekDouble> v(nq, 0.0);
2604
2605 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2607 {
2608 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2609 }
2610
2611 Array<OneD, NekDouble> x(nq);
2612 Array<OneD, NekDouble> y(nq);
2613 Array<OneD, NekDouble> z(nq);
2614
2616
2617 NekDouble R = 4.0;
2618 NekDouble cos2theta, cosRtheta, cos6theta, cos2Rtheta, cosRm1theta;
2619 NekDouble cos2phi, cos4phi, sin4phi, cos8phi;
2620
2621 // disturbancees of Rossby-Haurwitz Wave
2622 NekDouble x0d, y0d, z0d, phi0, theta0;
2623 // NekDouble rad_earth = 6.37122 * 1000000;
2624
2625 phi0 = 40.0 * m_pi / 180.0;
2626 theta0 = 50.0 * m_pi / 180.0;
2627
2628 x0d = cos(phi0) * cos(theta0);
2629 y0d = sin(phi0) * cos(theta0);
2630 z0d = sin(theta0);
2631
2632 for (int j = 0; j < nq; ++j)
2633 {
2634 x0j = x[j];
2635 x1j = y[j];
2636 x2j = z[j];
2637
2638 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2639 cos_theta);
2640
2641 // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2642 // \sin \alpha + \sin \theta \cos \alpha )^2
2643
2644 // tmp = cos^{2R} \theta, R = 4;
2645 cos2theta = cos_theta * cos_theta;
2646 cosRm1theta = cos_theta * cos2theta;
2647 cosRtheta = cos2theta * cos2theta;
2648 cos6theta = cos2theta * cosRtheta;
2649 cos2Rtheta = cosRtheta * cosRtheta;
2650
2652 Ath = tmp * cos2theta +
2654 ((R + 1.0) * cosRtheta + (2 * R * R - R - 2.0) * cos2theta -
2655 2 * R * R);
2656
2658 Bth = tmp * cosRtheta *
2659 ((R * R + 2 * R + 2) - (R + 1.0) * (R + 1.0) * cos2theta);
2660
2661 Cth =
2663
2664 // cos (2 \phi) = 2 * cos \phi * cos \phi - 1.0
2665 cos2phi = 2.0 * cos_varphi * cos_varphi - 1.0;
2666 cos4phi = 2.0 * cos2phi * cos2phi - 1.0;
2667 cos8phi = 2.0 * cos4phi * cos4phi - 1.0;
2668
2669 // sin (2 \phi) = 2 * cos \phi * sin \phi
2670 sin4phi = 4.0 * sin_varphi * cos_varphi * cos2phi;
2671
2673
2674 // disturbances is added
2676 {
2677 eta[j] = eta[j] *
2678 (1.0 + (1.0 / 40.0) * (x0j * x0d + x1j * y0d + x2j * z0d));
2679 }
2680
2681 // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2682 // || e^2 ||^2
2683 uhat = m_angfreq * cos_theta +
2684 m_K * cosRm1theta * (R * sin_theta * sin_theta - cos2theta) *
2685 cos4phi;
2686 vhat = -1.0 * m_K * R * cosRm1theta * sin_theta * sin4phi;
2687
2688 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2689 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2690 uvec[2][j] = vhat * cos_theta;
2691 }
2692
2693 // NekDouble etamin, etaaver;
2694 // etamin = Vmath::Vmin(nq, eta, 1)*rad_earth;
2695 // etaaver = Vmath::Vsum(nq, eta, 1)/nq*rad_earth;
2696
2697 // Projection of u onto the tangent plane with conserving the mag. of the
2698 // velocity.
2699 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2700
2702 {
2703 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2704 }
2705
2706 // u is projected on the tangent plane with preserving its length
2707 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2708
2709 // Change it to the coordinate of moving frames
2710 // CartesianToMovingframes(0,uvecproj,u);
2711 // CartesianToMovingframes(1,uvecproj,v);
2712
2713 u = CartesianToMovingframes(uvec, 0);
2714 v = CartesianToMovingframes(uvec, 1);
2715
2717 {
2718 case (0):
2719 {
2720 outfield = eta;
2721 }
2722 break;
2723
2724 case (1):
2725 {
2726 outfield = u;
2727 }
2728 break;
2729
2730 case (2):
2731 {
2732 outfield = v;
2733 }
2734 break;
2735 }
2736}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), FilterPython_Function::field, 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 1091 of file MMFSWE.cpp.
1095{
1096 int nvariables = 3;
1097 NekDouble MageF1, MageF2, MageB1, MageB2;
1098 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1099 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1100
1102 NekDouble uRF, vRF, uLB, vLB;
1103 NekDouble velL, velR;
1104
1105 Array<OneD, NekDouble> EigF(nvariables);
1106 Array<OneD, NekDouble> EigB(nvariables);
1107
1108 // Compute Magnitude and Dot product of moving frames for the index
1109 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1110 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1111
1112 // Get the velocity in the normal to the edge
1115
1116 NekDouble SL, SR;
1117 SL = fabs(velL) + sqrt(g * hL);
1118 SR = fabs(velR) + sqrt(g * hR);
1119
1120 NekDouble S;
1121 if (SL > SR)
1122 {
1123 S = SL;
1124 }
1125 else
1126 {
1127 S = SR;
1128 }
1129
1130 // uRF = uR component in moving frames e^{Fwd}
1131 // vRF = vR component in moving frames e^{Fwd}
1132 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1133 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1134
1135 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1136 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1137 numfluxF[2] = 0.5 * S * (hL - hR);
1138
1139 numfluxF[3] =
1140 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
1141 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1142 numfluxF[5] = 0.5 * S * (uL * hL - uRF * hR);
1143
1144 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1145 numfluxF[7] =
1146 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1147 numfluxF[8] = 0.5 * S * (vL * hL - vRF * hR);
1148
1149 // uLB = uL component in moving frames e^{Bwd}
1150 // vLB = vL component in moving frames e^{Bwd}
1151 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1152 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1153
1154 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1155 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1156 numfluxB[2] = 0.5 * S * (hL - hR);
1157
1158 numfluxB[3] =
1159 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1160 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1161 numfluxB[5] = 0.5 * S * (uLB * hL - uR * hR);
1162
1163 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1164 numfluxB[7] =
1165 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1166 numfluxB[8] = 0.5 * S * (vLB * hL - vR * hR);
1167}
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 1383 of file MMFSWE.cpp.
1385{
1386
1388 int cnt = 0;
1389
1390 // loop over Boundary Regions
1392 {
1393
1394 // Zonal Boundary Condition
1396 {
1398 {
1400 }
1402 {
1403 // ZonalBoundary2D(n,cnt,inarray);
1404 }
1405 }
1406
1407 // Wall Boundary Condition
1409 {
1411 {
1413 }
1415 {
1416 WallBoundary2D(n, cnt, inarray);
1417 }
1418 }
1419
1420 // Time Dependent Boundary Condition (specified in meshfile)
1422 "eTimeDependent")
1423 {
1424 for (int i = 0; i < nvariables; ++i)
1425 {
1426 m_fields[i]->EvaluateBoundaryConditions(time);
1427 }
1428 }
1429 cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
1430 }
1431}
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Definition: MMFSWE.cpp:1433
References ASSERTL0, Nektar::SolverUtils::EquationSystem::m_expdim, Nektar::SolverUtils::EquationSystem::m_fields, and WallBoundary2D().
Referenced by DoOdeProjection().
◆ SteadyZonalFlow()
|
protected |
Definition at line 2053 of file MMFSWE.cpp.
2055{
2057 NekDouble uhat, vhat;
2058 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2059 NekDouble x0j, x1j, x2j, tmp;
2060
2061 Array<OneD, NekDouble> eta(nq, 0.0);
2062 Array<OneD, NekDouble> u(nq, 0.0);
2063 Array<OneD, NekDouble> v(nq, 0.0);
2064
2065 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2067 {
2068 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2069 }
2070
2071 Array<OneD, NekDouble> x(nq);
2072 Array<OneD, NekDouble> y(nq);
2073 Array<OneD, NekDouble> z(nq);
2074
2076
2077 for (int j = 0; j < nq; ++j)
2078 {
2079 x0j = x[j];
2080 x1j = y[j];
2081 x2j = z[j];
2082
2083 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2084 cos_theta);
2085
2086 // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2087 // \sin \alpha + \sin \theta \cos \alpha )^2
2088 tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
2089 sin_theta * cos(m_alpha);
2091
2092 // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2093 // || e^2 ||^2
2095 sin_theta * cos_varphi * sin(m_alpha));
2097
2098 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2099 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2100 uvec[2][j] = vhat * cos_theta;
2101 }
2102
2103 // Projection of u onto the tangent plane with conserving the mag. of the
2104 // velocity.
2105 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2106
2108 {
2109 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2110 }
2111
2112 // u is projected on the tangent plane with preserving its length
2113 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2114
2115 // Change it to the coordinate of moving frames
2116 // CartesianToMovingframes(0,uvecproj,u);
2117 // CartesianToMovingframes(1,uvecproj,v);
2118
2119 u = CartesianToMovingframes(uvec, 0);
2120 v = CartesianToMovingframes(uvec, 1);
2121
2123 {
2124 case (0):
2125 {
2126 outfield = eta;
2127 }
2128 break;
2129
2130 case (1):
2131 {
2132 outfield = u;
2133 }
2134 break;
2135
2136 case (2):
2137 {
2138 outfield = v;
2139 }
2140 break;
2141 }
2142}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), FilterPython_Function::field, 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 1989 of file MMFSWE.cpp.
1992{
1994
1995 Array<OneD, NekDouble> x0(nq);
1996 Array<OneD, NekDouble> x1(nq);
1997 Array<OneD, NekDouble> x2(nq);
1998
1999 m_fields[0]->GetCoords(x0, x1, x2);
2000
2001 Array<OneD, NekDouble> eta0(nq);
2002 Array<OneD, NekDouble> u0(nq);
2003 Array<OneD, NekDouble> v0(nq);
2004
2005 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2006
2008 {
2009 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2010 }
2011
2012 for (int i = 0; i < nq; ++i)
2013 {
2014 eta0[i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2015 (1.0 / cosh(0.395 * x0[i]))) *
2016 (3.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2017 exp(-0.5 * x1[i] * x1[i]);
2018 uvec[0][i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2019 (1.0 / cosh(0.395 * x0[i]))) *
2020 (-9.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2021 exp(-0.5 * x1[i] * x1[i]);
2022 uvec[1][i] = (-2.0 * 0.395 * tanh(0.395 * x0[i])) *
2023 (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2024 (1.0 / cosh(0.395 * x0[i]))) *
2025 (2.0 * x1[i]) * exp(-0.5 * x1[i] * x1[i]);
2026 }
2027
2028 u0 = CartesianToMovingframes(uvec, 0);
2029 v0 = CartesianToMovingframes(uvec, 1);
2030
2032 {
2033 case (0):
2034 {
2035 outfield = eta0;
2036 }
2037 break;
2038
2039 case (1):
2040 {
2041 outfield = u0;
2042 }
2043 break;
2044
2045 case (2):
2046 {
2047 outfield = v0;
2048 }
2049 break;
2050 }
2051}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::m_fields, and Nektar::SolverUtils::EquationSystem::m_spacedim.
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ TestVorticityComputation()
Definition at line 2947 of file MMFSWE.cpp.
2948{
2949 // Construct beta
2950 int i, k;
2951 int n = 1, m = 1;
2953
2954 NekDouble alpha, beta_theta, beta_phi;
2955
2956 NekDouble xp, yp, zp, Re;
2957 NekDouble theta, phi, sin_theta, cos_theta, sin_varphi, cos_varphi;
2958 NekDouble cosntheta3;
2959
2960 NekDouble thetax, thetay, thetaz;
2961 NekDouble phix, phiy, phiz;
2962
2963 Array<OneD, NekDouble> x(nq);
2964 Array<OneD, NekDouble> y(nq);
2965 Array<OneD, NekDouble> z(nq);
2966
2967 Array<OneD, NekDouble> u(nq);
2968 Array<OneD, NekDouble> v(nq);
2969
2970 Array<OneD, NekDouble> vorticitycompt(nq);
2971 Array<OneD, NekDouble> vorticityexact(nq);
2972
2974
2975 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2977 {
2978 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2979 }
2980
2981 // Generate SphericalCoords
2982 for (k = 0; k < nq; ++k)
2983 {
2984 xp = x[k];
2985 yp = y[k];
2986 zp = z[k];
2987
2988 Re = sqrt(xp * xp + yp * yp + zp * zp);
2989
2990 CartesianToSpherical(xp, yp, zp, sin_varphi, cos_varphi, sin_theta,
2991 cos_theta);
2992
2993 alpha = sin_varphi;
2994
2995 theta = atan2(sin_theta, cos_theta);
2996 phi = atan2(sin_varphi, cos_varphi);
2997
2998 cosntheta3 = cos(n * theta) * cos(n * theta) * cos(n * theta);
2999
3000 beta_theta = -4.0 * n * cosntheta3 * cos(m * phi) * sin(n * theta) / Re;
3001 beta_phi = -m * cosntheta3 * sin(m * phi) / Re;
3002
3003 thetax = -1.0 * cos_varphi * sin_theta;
3004 thetay = -1.0 * sin_varphi * sin_theta;
3005 thetaz = cos_theta;
3006
3007 phix = -1.0 * sin_varphi;
3008 phiy = cos_varphi;
3009 phiz = 0.0;
3010
3011 uvec[0][k] = alpha * (beta_theta * thetax + beta_phi * phix);
3012 uvec[1][k] = alpha * (beta_theta * thetay + beta_phi * phiy);
3013 uvec[2][k] = alpha * (beta_theta * thetaz + beta_phi * phiz);
3014
3015 vorticityexact[k] = -4.0 * n / Re / Re * cos_theta * cos_theta *
3016 cos_varphi * cos(m * phi) * sin(n * theta);
3017 }
3018
3019 u = CartesianToMovingframes(uvec, 0);
3020 v = CartesianToMovingframes(uvec, 1);
3021
3022 std::cout << "chi migi1" << std::endl;
3023
3024 ComputeVorticity(u, v, vorticitycompt);
3025
3026 Vmath::Vsub(nq, vorticityexact, 1, vorticitycompt, 1, vorticitycompt, 1);
3027
3030 << std::endl;
3031}
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 2484 of file MMFSWE.cpp.
2487{
2489
2490 NekDouble uhat, vhat;
2491 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2492 NekDouble x0j, x1j, x2j;
2493 NekDouble Ttheta, Tphi;
2494
2495 Array<OneD, NekDouble> eta(nq, 0.0);
2496 Array<OneD, NekDouble> u(nq, 0.0);
2497 Array<OneD, NekDouble> v(nq, 0.0);
2498
2499 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2501 {
2502 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2503 }
2504
2505 Array<OneD, NekDouble> x(nq);
2506 Array<OneD, NekDouble> y(nq);
2507 Array<OneD, NekDouble> z(nq);
2508
2510
2511 int Nint = 1000;
2512 NekDouble dth, intj;
2513 for (int j = 0; j < nq; ++j)
2514 {
2515 x0j = x[j];
2516 x1j = y[j];
2517 x2j = z[j];
2518
2519 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2520 cos_theta);
2521
2522 Ttheta = atan2(sin_theta, cos_theta);
2523 Tphi = atan2(sin_varphi, cos_varphi);
2524
2525 uhat = ComputeUnstableJetuphi(Ttheta);
2526 vhat = 0.0;
2527
2528 // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2529 dth = Ttheta / Nint;
2530 eta[j] = dth * 0.5 *
2532 for (int i = 1; i < Nint - 1; i++)
2533 {
2534 intj = i * dth;
2535 eta[j] = eta[j] + dth * ComputeUnstableJetEta(intj);
2536 }
2537 eta[j] = (-1.0 / m_g) * eta[j];
2538
2539 // Add perturbation
2541 {
2542 eta[j] = eta[j] + m_hbar * cos_theta * exp(-9.0 * Tphi * Tphi) *
2543 exp(-225.0 * (m_pi / 4.0 - Ttheta) *
2544 (m_pi / 4.0 - Ttheta));
2545 }
2546
2547 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2548 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2549 uvec[2][j] = vhat * cos_theta;
2550 }
2551
2552 // Projection of u onto the tangent plane with conserving the mag. of the
2553 // velocity.
2554 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2555
2557 {
2558 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2559 }
2560
2561 // u is projected on the tangent plane with preserving its length
2562 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2563
2564 // Change it to the coordinate of moving frames
2565 // CartesianToMovingframes(0,uvecproj,u);
2566 // CartesianToMovingframes(1,uvecproj,v);
2567
2568 u = CartesianToMovingframes(uvec, 0);
2569 v = CartesianToMovingframes(uvec, 1);
2570
2572 {
2573 case (0):
2574 {
2575 outfield = eta;
2576 }
2577 break;
2578
2579 case (1):
2580 {
2581 outfield = u;
2582 }
2583 break;
2584
2585 case (2):
2586 {
2587 outfield = v;
2588 }
2589 break;
2590 }
2591}
NekDouble ComputeUnstableJetEta(const NekDouble theta)
Definition: MMFSWE.cpp:2738
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), ComputeUnstableJetEta(), ComputeUnstableJetuphi(), FilterPython_Function::field, 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 2295 of file MMFSWE.cpp.
2297{
2299 NekDouble uhat, vhat;
2300 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2301 NekDouble x0j, x1j, x2j, tmp;
2302
2303 NekDouble TR, Ttheta;
2304
2305 Array<OneD, NekDouble> eta(nq, 0.0);
2306 Array<OneD, NekDouble> u(nq, 0.0);
2307 Array<OneD, NekDouble> v(nq, 0.0);
2308
2309 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2311 {
2312 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2313 }
2314
2315 Array<OneD, NekDouble> x(nq);
2316 Array<OneD, NekDouble> y(nq);
2317 Array<OneD, NekDouble> z(nq);
2318
2320
2321 for (int j = 0; j < nq; ++j)
2322 {
2323 x0j = x[j];
2324 x1j = y[j];
2325 x2j = z[j];
2326
2327 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2328 cos_theta);
2329
2330 // \eta = ( - ( u_0 ( - T_R sin \alpha \cos \theta + \cos \alpha \sin
2331 // \theta ) + \Omega \sin \theta )^2 + \Omega \sin \theta )^2 + (\Omega
2332 // \sin \theta )^2
2333 TR =
2335 Ttheta =
2338
2342 eta[j] = 0.5 * eta[j] / m_g;
2343
2344 // u = u_0*(TR*\sin \alpha * \sin \theta + \cos \alpha * \cos \theta
2345 // v = - u_0 Ttheta * \sin \alpha
2346 uhat =
2349
2350 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2351 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2352 uvec[2][j] = vhat * cos_theta;
2353 }
2354
2355 // Projection of u onto the tangent plane with conserving the mag. of the
2356 // velocity.
2357 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2358
2360 {
2361 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2362 }
2363
2364 // u is projected on the tangent plane with preserving its length
2365 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2366
2367 // Change it to the coordinate of moving frames
2368 // CartesianToMovingframes(0,uvecproj,u);
2369 // CartesianToMovingframes(1,uvecproj,v);
2370
2371 u = CartesianToMovingframes(uvec, 0);
2372 v = CartesianToMovingframes(uvec, 1);
2373
2375 {
2376 case (0):
2377 {
2378 outfield = eta;
2379 }
2380 break;
2381
2382 case (1):
2383 {
2384 outfield = u;
2385 }
2386 break;
2387
2388 case (2):
2389 {
2390 outfield = v;
2391 }
2392 break;
2393 }
2394}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), FilterPython_Function::field, 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 1546 of file MMFSWE.cpp.
1547{
1548 // Compute m_depth and m_Derivdepth
1549 EvaluateWaterDepth();
1550 EvaluateCoriolis();
1551 SetInitialConditions(0.0, dumpInitialConditions);
1552 PrimitiveToConservative();
1553
1554 // transfer the initial conditions to modal values
1556 {
1559 m_fields[i]->UpdateCoeffs());
1560 }
1561}
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:171
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 242 of file MMFSWE.cpp.
243{
245
246 int i, nchk = 1;
247 int nvariables = 0;
250
252 {
253 for (i = 0; i < nfields; ++i)
254 {
255 m_intVariables.push_back(i);
256 }
257 nvariables = nfields;
258 }
259 else
260 {
261 nvariables = m_intVariables.size();
262 }
263
264 // Set up wrapper to fields data storage.
265 Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
266 Array<OneD, Array<OneD, NekDouble>> tmp(nvariables);
267
268 // Order storage to list time-integrated fields first.
269 for (i = 0; i < nvariables; ++i)
270 {
273 }
274
275 // Initialise time integration scheme
277
278 // Check uniqueness of checkpoint output
282 "Only one of IO_CheckTime and IO_CheckSteps "
283 "should be set!");
284
285 LibUtilities::Timer timer;
286 bool doCheckTime = false;
287 int step = 0;
288 NekDouble intTime = 0.0;
289 NekDouble cpuTime = 0.0;
290 NekDouble elapsed = 0.0;
291
292 NekDouble Mass = 0.0, Energy = 0.0, Enstrophy = 0.0, Vorticity = 0.0;
293 Array<OneD, NekDouble> zeta(nq);
294 Array<OneD, Array<OneD, NekDouble>> fieldsprimitive(nvariables);
295 for (int i = 0; i < nvariables; ++i)
296 {
297 fieldsprimitive[i] = Array<OneD, NekDouble>(nq);
298 }
299
301 {
302 timer.Start();
304 timer.Stop();
305
307 elapsed = timer.TimePerTest(1);
308 intTime += elapsed;
309 cpuTime += elapsed;
310
311 // Write out status information
313 m_session->GetComm()->GetRank() == 0)
314 {
315 std::cout << "Steps: " << std::setw(8) << std::left << step + 1
316 << " "
318
319 std::stringstream ss;
320 ss << cpuTime << "s";
321 std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
322 << std::endl;
323
324 // Printout Mass, Energy, Enstrophy
325 ConservativeToPrimitive(fields, fieldsprimitive);
326
327 // Vorticity zeta
328 ComputeVorticity(fieldsprimitive[1], fieldsprimitive[2], zeta);
330
331 // Masss = h^*
333
334 // Energy = 0.5*( h^*(u^2 + v^2) + g ( h^2 - h_s^s ) )
335 Energy = (ComputeEnergy(fieldsprimitive[0], fieldsprimitive[1],
336 fieldsprimitive[2]) -
337 m_Energy0) /
338 m_Energy0;
339
340 // Enstrophy = 0.5/h^* ( \mathbf{k} \cdot (\nabla \times \mathbf{v}
341 // ) + f )^2
342 Enstrophy =
343 (ComputeEnstrophy(fieldsprimitive[0], fieldsprimitive[1],
344 fieldsprimitive[2]) -
345 m_Enstrophy0) /
346 m_Enstrophy0;
347
348 std::cout << "dMass = " << std::setw(8) << std::left << Mass << " "
349 << ", dEnergy = " << std::setw(8) << std::left << Energy
350 << " "
351 << ", dEnstrophy = " << std::setw(8) << std::left
352 << Enstrophy << " "
353 << ", dVorticity = " << std::setw(8) << std::left
354 << Vorticity << std::endl
355 << std::endl;
356
357 cpuTime = 0.0;
358 }
359
360 // (h, hu, hv) -> (\eta, u, v)
361 ConservativeToPrimitive();
362
363 // Transform data into coefficient space
364 for (i = 0; i < nvariables; ++i)
365 {
370 }
371
372 // Write out checkpoint files
374 doCheckTime)
375 {
376 Checkpoint_Output(nchk++);
377 doCheckTime = false;
378
379 // (\eta, u, v) -> (h, hu, hv)
380 PrimitiveToConservative();
381 }
382
383 // Step advance
384 ++step;
385 }
386
387 // Print out summary statistics
389 {
391 {
393 << std::endl
395 }
396
398 {
399 std::cout << "Time-integration : " << intTime << "s" << std::endl;
400 }
401 }
402
403 for (i = 0; i < nvariables; ++i)
404 {
407 }
408
409 // (h, hu, hv) -> (\eta, u, v)
410 ConservativeToPrimitive();
411
412 for (i = 0; i < nvariables; ++i)
413 {
415 m_fields[i]->UpdateCoeffs());
416 }
417}
NekDouble ComputeEnstrophy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2183
NekDouble ComputeMass(const Array< OneD, const NekDouble > &eta)
Definition: MMFSWE.cpp:2144
NekDouble ComputeEnergy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2154
int m_infosteps
Number of time steps between outputting status information.
Definition: EquationSystem.h:498
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition: EquationSystem.cpp:1253
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:89
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
Definition: UnsteadySystem.h:97
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
Definition: UnsteadySystem.h:86
std::vector< int > m_intVariables
Definition: UnsteadySystem.h:94
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:49
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 3219 of file MMFSWE.cpp.
3222{
3224 {
3226 {
3228 }
3229 break;
3230
3232 {
3234 }
3235 break;
3236
3238 {
3240 }
3241 break;
3242
3244 {
3246 }
3247 break;
3248
3250 {
3252 }
3253 break;
3254
3256 {
3258 }
3259 break;
3260
3261 default:
3262 break;
3263 }
3264}
void UnstableJetFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2484
void RossbyWave(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2593
void IsolatedMountainFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2396
void UnsteadyZonalFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2295
void TestSWE2Dproblem(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:1989
void SteadyZonalFlow(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2053
References Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, FilterPython_Function::field, 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 3266 of file MMFSWE.cpp.
3267{
3268 MMFSystem::v_GenerateSummary(s);
3270
3272 {
3274 {
3276 }
3277 break;
3278
3280 {
3282 m_RossbyDisturbance);
3283 }
3284 break;
3285
3287 {
3289 }
3290 break;
3291
3292 default:
3293 break;
3294 }
3295}
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Virtual function for generating summary information.
Definition: MMFSystem.cpp:2463
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:47
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:92
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:84
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:419
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:1364
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:51
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
Definition: UnsteadySystem.h:109
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:85
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 3033 of file MMFSWE.cpp.
3036{
3038 NekDouble L2error = -1.0;
3039
3041 {
3043 {
3046 }
3047
3049 {
3050 case (0):
3051 {
3052 // I(h) = I (h - h_exact ) / I (h_exact)
3053 Array<OneD, NekDouble> exactsolution(nq);
3055
3056 // exactsoln = u - u_T so that L2 compute u_T
3058
3060 &exactsolution[0], 1, &exactsolution[0], 1);
3061 Vmath::Vabs(nq, exactsolution, 1, exactsolution, 1);
3062
3063 L2error = (m_fields[0]->Integral(exactsolution)) / L2exact;
3064 }
3065 break;
3066
3067 case (1):
3068 {
3069 // I2 (u) = I( (u - u_ext)^2 + (v - v_ext)^2 )^{1/2} / I(
3070 // u_ext^2 + v_ext^2 )^{1/2}
3071 Array<OneD, NekDouble> exactu(nq);
3072 Array<OneD, NekDouble> exactv(nq);
3073 Array<OneD, NekDouble> tmp(nq);
3074
3075 // L2exact = \int (\sqrt{exactu*exactu+exactv*exactv})
3076 NekDouble L2exact;
3079 Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3080 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3081 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3082
3083 L2exact = m_fields[1]->Integral(tmp);
3084
3085 // L2exact = \int
3086 // (\sqrt{(u-exactu)*(u-exactu)+(v-exactv)*(v-exactv)})
3088 &exactu[0], 1);
3090 &exactv[0], 1);
3091 Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3092 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3093 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3094
3095 L2error = (m_fields[1]->Integral(tmp)) / L2exact;
3096 }
3097 break;
3098
3099 case (2):
3100 {
3101 L2error = 0.0;
3102 }
3103 break;
3104
3105 default:
3106 break;
3107 }
3108
3109 if (Normalised == true)
3110 {
3112
3115
3116 L2error = sqrt(L2error * L2error / Vol);
3117 }
3118 }
3119 else
3120 {
3121 Array<OneD, NekDouble> L2INF(2);
3123 L2error = L2INF[0];
3124 }
3125
3126 return L2error;
3127}
void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
Definition: MMFSWE.cpp:3219
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:351
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:974
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Definition: EquationSystem.h:545
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.hpp:352
References Nektar::SolverUtils::EquationSystem::ErrorExtraPoints(), FilterPython_Function::field, 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 3129 of file MMFSWE.cpp.
3132{
3134
3136 {
3139 }
3140
3142
3143 // Obtain \vec{u} in cartesian coordinate
3144 Array<OneD, NekDouble> x(nq);
3145 Array<OneD, NekDouble> y(nq);
3146 Array<OneD, NekDouble> z(nq);
3147
3149
3151 {
3152 case (0):
3153 {
3154 NekDouble Letaint;
3155
3156 Array<OneD, NekDouble> exactsolution(nq);
3157
3160 exactsolution);
3161
3162 Letaint = Vmath::Vamax(nq, exactsolution, 1);
3163
3165 1, &exactsolution[0], 1);
3166
3168 }
3169 break;
3170
3171 case (1):
3172 {
3173 Array<OneD, NekDouble> exactu(nq);
3174 Array<OneD, NekDouble> exactv(nq);
3175 Array<OneD, NekDouble> tmpu(nq);
3176 Array<OneD, NekDouble> tmpv(nq);
3177 Array<OneD, NekDouble> Lerr(nq);
3178 Array<OneD, NekDouble> uT(nq);
3179
3182
3183 // Compute max[sqrt{(u-uex)^2 + (v-vex)^2}]
3186
3187 Vmath::Vsub(nq, &exactu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3188 Vmath::Vsub(nq, &exactv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3189
3190 Vmath::Vmul(nq, &tmpu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3191 Vmath::Vmul(nq, &tmpv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3192
3193 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &Lerr[0], 1);
3194 Vmath::Vsqrt(nq, &Lerr[0], 1, &Lerr[0], 1);
3195
3196 // uT = max[sqrt( u_T^2 + v_T^2 ) ]
3197 Vmath::Vmul(nq, &exactu[0], 1, &exactu[0], 1, &tmpu[0], 1);
3198 Vmath::Vmul(nq, &exactv[0], 1, &exactv[0], 1, &tmpv[0], 1);
3199 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &uT[0], 1);
3200 Vmath::Vsqrt(nq, &uT[0], 1, &uT[0], 1);
3201
3203 }
3204 break;
3205
3206 case (2):
3207 {
3208 LinfError = 0.0;
3209 }
3210 break;
3211
3212 default:
3213 break;
3214 }
3215
3217}
SOLVER_UTILS_EXPORT NekDouble LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Linf error computation.
Definition: EquationSystem.h:203
SOLVER_UTILS_EXPORT void EvaluateExactSolution(int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Evaluates an exact solution.
Definition: EquationSystem.h:179
References Nektar::SolverUtils::EquationSystem::EvaluateExactSolution(), FilterPython_Function::field, 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 1792 of file MMFSWE.cpp.
1795{
1797
1799 {
1801 {
1802 Array<OneD, NekDouble> eta0(nq);
1803 Array<OneD, NekDouble> u0(nq);
1804 Array<OneD, NekDouble> v0(nq);
1805
1806 TestSWE2Dproblem(initialtime, 0, eta0);
1807 m_fields[0]->SetPhys(eta0);
1808
1809 TestSWE2Dproblem(initialtime, 1, u0);
1810 m_fields[1]->SetPhys(u0);
1811
1812 TestSWE2Dproblem(initialtime, 2, v0);
1813 m_fields[2]->SetPhys(v0);
1814
1815 // forward transform to fill the modal coeffs
1817 {
1820 m_fields[i]->UpdateCoeffs());
1821 }
1822 }
1823 break;
1824
1826 {
1827 Array<OneD, NekDouble> eta0(nq);
1828 Array<OneD, NekDouble> u0(nq);
1829 Array<OneD, NekDouble> v0(nq);
1830 Array<OneD, NekDouble> zeta0(nq);
1831
1832 SteadyZonalFlow(0, eta0);
1833 m_fields[0]->SetPhys(eta0);
1834
1835 SteadyZonalFlow(1, u0);
1836 m_fields[1]->SetPhys(u0);
1837
1838 SteadyZonalFlow(2, v0);
1839 m_fields[2]->SetPhys(v0);
1840
1841 // ComputeVorticity(u0, v0, zeta0);
1843
1847
1848 // forward transform to fill the modal coeffs
1850 {
1853 m_fields[i]->UpdateCoeffs());
1854 }
1855 }
1856 break;
1857
1859 {
1860 Array<OneD, NekDouble> eta0(nq);
1861 Array<OneD, NekDouble> u0(nq);
1862 Array<OneD, NekDouble> v0(nq);
1863
1864 UnsteadyZonalFlow(0, initialtime, eta0);
1865 m_fields[0]->SetPhys(eta0);
1866
1867 UnsteadyZonalFlow(1, initialtime, u0);
1868 m_fields[1]->SetPhys(u0);
1869
1870 UnsteadyZonalFlow(2, initialtime, v0);
1871 m_fields[2]->SetPhys(v0);
1872
1876
1877 // forward transform to fill the modal coeffs
1879 {
1882 m_fields[i]->UpdateCoeffs());
1883 }
1884 }
1885 break;
1886
1888 {
1889 Array<OneD, NekDouble> eta0(nq);
1890 Array<OneD, NekDouble> u0(nq);
1891 Array<OneD, NekDouble> v0(nq);
1892
1893 IsolatedMountainFlow(0, initialtime, eta0);
1894 m_fields[0]->SetPhys(eta0);
1895
1896 IsolatedMountainFlow(1, initialtime, u0);
1897 m_fields[1]->SetPhys(u0);
1898
1899 IsolatedMountainFlow(2, initialtime, v0);
1900 m_fields[2]->SetPhys(v0);
1901
1905
1906 // forward transform to fill the modal coeffs
1908 {
1911 m_fields[i]->UpdateCoeffs());
1912 }
1913 }
1914 break;
1915
1917 {
1918 Array<OneD, NekDouble> eta0(nq);
1919 Array<OneD, NekDouble> u0(nq);
1920 Array<OneD, NekDouble> v0(nq);
1921
1922 UnstableJetFlow(0, initialtime, eta0);
1923 m_fields[0]->SetPhys(eta0);
1924
1925 UnstableJetFlow(1, initialtime, u0);
1926 m_fields[1]->SetPhys(u0);
1927
1928 UnstableJetFlow(2, initialtime, v0);
1929 m_fields[2]->SetPhys(v0);
1930
1934
1935 // forward transform to fill the modal coeffs
1937 {
1940 m_fields[i]->UpdateCoeffs());
1941 }
1942 }
1943 break;
1944
1946 {
1947 Array<OneD, NekDouble> eta0(nq);
1948 Array<OneD, NekDouble> u0(nq);
1949 Array<OneD, NekDouble> v0(nq);
1950
1951 RossbyWave(0, eta0);
1952 m_fields[0]->SetPhys(eta0);
1953
1954 RossbyWave(1, u0);
1955 m_fields[1]->SetPhys(u0);
1956
1957 RossbyWave(2, v0);
1958 m_fields[2]->SetPhys(v0);
1959
1963
1964 // forward transform to fill the modal coeffs
1966 {
1969 m_fields[i]->UpdateCoeffs());
1970 }
1971 }
1972 break;
1973
1974 default:
1975 break;
1976 }
1977
1978 if (dumpInitialConditions)
1979 {
1980 // dump initial conditions to file
1982 WriteFld(outname);
1983
1985 Checkpoint_Output_Cartesian(outname);
1986 }
1987}
void Checkpoint_Output_Cartesian(std::string outname)
Definition: MMFSWE.cpp:2768
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 1433 of file MMFSWE.cpp.
1435{
1436
1437 int i;
1439 int nvariables = physarray.size();
1440
1441 // get physical values of the forward trace
1442 Array<OneD, Array<OneD, NekDouble>> Fwd0(nvariables);
1443 for (i = 0; i < nvariables; ++i)
1444 {
1445 Fwd0[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1446 m_fields[i]->ExtractTracePhys(physarray[i], Fwd0[i]);
1447 }
1448
1449 Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
1450 for (i = 0; i < nvariables; ++i)
1451 {
1452 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1453 m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1454 }
1455
1456 // Adjust the physical values of the trace to take
1457 // user defined boundaries into account
1458 int e, id1, id2, npts;
1459
1461 ++e)
1462 {
1463 npts = m_fields[0]
1464 ->GetBndCondExpansions()[bcRegion]
1465 ->GetExp(e)
1466 ->GetTotPoints();
1467 id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
1468 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
1469 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1470
1472 {
1473 case 1:
1474 {
1475 // negate the forward flux
1476 Vmath::Neg(npts, &Fwd[1][id2], 1);
1477 }
1478 break;
1479 case 2:
1480 {
1481 Array<OneD, NekDouble> tmp_n(npts);
1482 Array<OneD, NekDouble> tmp_t(npts);
1483
1485 &tmp_n[0], 1);
1487 &tmp_n[0], 1, &tmp_n[0], 1);
1488
1490 &tmp_t[0], 1);
1492 1, &tmp_t[0], 1, &tmp_t[0], 1);
1493
1494 // negate the normal flux
1495 Vmath::Neg(npts, tmp_n, 1);
1496
1497 Array<OneD, NekDouble> denom(npts);
1498 Array<OneD, NekDouble> tmp_u(npts);
1499 Array<OneD, NekDouble> tmp_v(npts);
1500
1501 // denom = (e^1 \cdot n ) (e^2 \cdot t) - (e^2 \cdot n ) (e^1
1502 // \cdot t)
1504 &m_nperpcdotMFFwd[0][id2], 1, &denom[0], 1);
1506 &m_nperpcdotMFFwd[1][id2], 1, &denom[0], 1,
1507 &denom[0], 1);
1508
1510 &tmp_u[0], 1);
1512 &tmp_u[0], 1, &tmp_u[0], 1);
1513 Vmath::Vdiv(npts, &tmp_u[0], 1, &denom[0], 1, &tmp_u[0], 1);
1514
1515 Vmath::Vcopy(npts, &tmp_u[0], 1, &Fwd[1][id2], 1);
1516
1518 &tmp_v[0], 1);
1520 &tmp_v[0], 1, &tmp_v[0], 1);
1521 Vmath::Vdiv(npts, &tmp_v[0], 1, &denom[0], 1, &tmp_v[0], 1);
1522
1523 Vmath::Vcopy(npts, &tmp_v[0], 1, &Fwd[2][id2], 1);
1524 }
1525 break;
1526
1527 default:
1529 }
1530
1531 // copy boundary adjusted values into the boundary expansion
1532 for (i = 0; i < nvariables; ++i)
1533 {
1534 Vmath::Vcopy(npts, &Fwd[i][id2], 1,
1535 &(m_fields[i]
1536 ->GetBndCondExpansions()[bcRegion]
1537 ->UpdatePhys())[id1],
1538 1);
1539 }
1540 }
1541}
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:326
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.hpp:381
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 473 of file MMFSWE.cpp.
476{
477 int i;
482
483 Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
484 Array<OneD, Array<OneD, NekDouble>> physfield(nvariables);
485
487 {
488 fluxvector[i] = Array<OneD, NekDouble>(nq);
489 }
490
491 // InField is Primitive
492 for (i = 0; i < nvariables; ++i)
493 {
494 physfield[i] = InField[i];
495 }
496
497 // Get the ith component of the flux vector in (physical space)
498 // fluxvector[0] = component for e^1 cdot \nabla \varphi
499 // fluxvector[1] = component for e^2 cdot \nabla \varphi
500 Array<OneD, NekDouble> fluxtmp(nq);
501 Array<OneD, NekDouble> tmp(ncoeffs);
502
503 // Compute Divergence Components
504 for (i = 0; i < nvariables; ++i)
505 {
506 GetSWEFluxVector(i, physfield, fluxvector);
507
508 OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
510 {
511 // Directional derivation with respect to the j'th moving frame
512 // tmp_j = ( \nabla \phi, fluxvector[j] \mathbf{e}^j )
514 fluxvector[j], tmp);
515 Vmath::Vadd(ncoeffs, &tmp[0], 1, &OutField[i][0], 1,
516 &OutField[i][0], 1);
517 }
518 }
519
520 // Numerical Flux
521 Array<OneD, Array<OneD, NekDouble>> numfluxFwd(nvariables);
522 Array<OneD, Array<OneD, NekDouble>> numfluxBwd(nvariables);
523
524 for (i = 0; i < nvariables; ++i)
525 {
526 numfluxFwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
527 numfluxBwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
528 }
529
530 NumericalSWEFlux(physfield, numfluxFwd, numfluxBwd);
531
532 // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
533 for (i = 0; i < nvariables; ++i)
534 {
535 Vmath::Neg(ncoeffs, OutField[i], 1);
536 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
537 OutField[i]);
539 }
540}
void NumericalSWEFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd)
Definition: MMFSWE.cpp:668
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:321
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: BasicUtils/NekFactory.hpp:197
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: MMFSWE.h:66
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:99
Name of class.
◆ m_AddCoriolis
|
protected |
Definition at line 88 of file MMFSWE.h.
Referenced by DoOdeRhs(), and v_InitObject().
◆ m_AddRotation
|
protected |
Definition at line 88 of file MMFSWE.h.
Referenced by DoOdeRhs(), and v_InitObject().
◆ m_alpha
|
protected |
Definition at line 91 of file MMFSWE.h.
Referenced by EvaluateCoriolisForZonalFlow(), SteadyZonalFlow(), UnsteadyZonalFlow(), v_GenerateSummary(), and v_InitObject().
◆ m_angfreq
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by RossbyWave(), and v_InitObject().
◆ m_coriolis
Coriolis force.
Definition at line 86 of file MMFSWE.h.
Referenced by AddCoriolis(), ComputeEnstrophy(), and EvaluateCoriolis().
◆ m_depth
Still water depth.
Definition at line 82 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 83 of file MMFSWE.h.
Referenced by AddElevationEffect(), and EvaluateWaterDepth().
◆ m_en
|
protected |
Definition at line 94 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_Energy0
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_Enstrophy0
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_g
|
protected |
Definition at line 90 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 91 of file MMFSWE.h.
Referenced by Checkpoint_Output_Cartesian(), ComputeEnergy(), EvaluateWaterDepth(), RossbyWave(), SteadyZonalFlow(), and v_InitObject().
◆ m_hbar
|
protected |
Definition at line 94 of file MMFSWE.h.
Referenced by UnstableJetFlow(), and v_InitObject().
◆ m_hs0
|
protected |
Definition at line 93 of file MMFSWE.h.
Referenced by EvaluateWaterDepth(), and v_InitObject().
◆ m_Hvar
|
protected |
Definition at line 91 of file MMFSWE.h.
Referenced by SteadyZonalFlow(), and v_InitObject().
◆ m_K
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by RossbyWave(), and v_InitObject().
◆ m_k2
|
protected |
Definition at line 92 of file MMFSWE.h.
Referenced by EvaluateWaterDepth(), and v_InitObject().
◆ m_Mass0
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_Omega
|
protected |
Definition at line 91 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 268 of file MMFSWE.h.
Referenced by UnstableJetFlow(), v_GenerateSummary(), and v_InitObject().
◆ m_RossbyDisturbance
|
private |
Definition at line 267 of file MMFSWE.h.
Referenced by RossbyWave(), v_GenerateSummary(), and v_InitObject().
◆ m_TestType
|
protected |
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 94 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_theta1
|
protected |
Definition at line 94 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_u0
|
protected |
Definition at line 91 of file MMFSWE.h.
Referenced by IsolatedMountainFlow(), SteadyZonalFlow(), UnsteadyZonalFlow(), and v_InitObject().
◆ m_uthetamax
|
protected |
Definition at line 94 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_veldotMF
◆ m_vellc
◆ m_velocity
Advection velocity.
Definition at line 103 of file MMFSWE.h.
Referenced by ComputeNablaCdotVelocity().
◆ m_Vorticity0
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
