#include <MMFSWE.h>
Inheritance diagram for Nektar::MMFSWE:
Static Public Member Functions | |
| static SolverUtils::EquationSystemSharedPtr | create (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph) |
| Creates an instance of this class. More... | |
Public Attributes | |
| TestType | m_TestType |
 Public Attributes inherited from Nektar::SolverUtils::MMFSystem | |
| NekDouble | m_pi |
| int | m_shapedim |
| SurfaceType | m_surfaceType |
| UpwindType | m_upwindType |
| TestMaxwellType | m_TestMaxwellType |
| PolType | m_PolType |
| IncType | m_IncType |
| Array< OneD, NekDouble > | m_MMFfactors |
Static Public Attributes | |
| static std::string | className |
| Name of class. More... | |
| Static Public Attributes inherited from Nektar::SolverUtils::UnsteadySystem | |
| static std::string | cmdSetStartTime |
| static std::string | cmdSetStartChkNum |
Protected Attributes | |
| Array< OneD, NekDouble > | m_depth |
| Still water depth. More... | |
| Array< OneD, Array< OneD, NekDouble > > | m_Derivdepth |
| Array< OneD, NekDouble > | m_coriolis |
| Coriolis force. More... | |
| int | m_AddCoriolis |
| int | m_AddRotation |
| NekDouble | m_g |
| NekDouble | m_alpha |
| NekDouble | m_u0 |
| NekDouble | m_Omega |
| NekDouble | m_H0 |
| NekDouble | m_Hvar |
| NekDouble | m_k2 |
| NekDouble | m_hs0 |
| NekDouble | m_uthetamax |
| NekDouble | m_theta0 |
| NekDouble | m_theta1 |
| NekDouble | m_en |
| NekDouble | m_hbar |
| NekDouble | m_angfreq |
| NekDouble | m_K |
| NekDouble | m_Vorticity0 |
| NekDouble | m_Mass0 |
| NekDouble | m_Energy0 |
| NekDouble | m_Enstrophy0 |
| bool | m_primitive |
| Indicates if variables are primitive or conservative. More... | |
| Array< OneD, Array< OneD, NekDouble > > | m_velocity |
| Advection velocity. More... | |
| Array< OneD, Array< OneD, NekDouble > > | m_veldotMF |
| Array< OneD, NekDouble > | m_vellc |
| int | m_planeNumber |
| Protected Attributes inherited from Nektar::SolverUtils::MMFSystem | |
| NekDouble | m_alpha |
| NekDouble | m_Incfreq |
| int | m_SmoothFactor |
| NekDouble | m_SFinit |
| Array< OneD, Array< OneD, NekDouble > > | m_movingframes |
| Array< OneD, Array< OneD, NekDouble > > | m_surfaceNormal |
| Array< OneD, Array< OneD, NekDouble > > | m_ncdotMFFwd |
| Array< OneD, Array< OneD, NekDouble > > | m_ncdotMFBwd |
| Array< OneD, Array< OneD, NekDouble > > | m_nperpcdotMFFwd |
| Array< OneD, Array< OneD, NekDouble > > | m_nperpcdotMFBwd |
| Array< OneD, Array< OneD, NekDouble > > | m_DivMF |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_CurlMF |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_MFtraceFwd |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_MFtraceBwd |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_ntimesMFFwd |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_ntimesMFBwd |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_ntimes_ntimesMFFwd |
| Array< OneD, Array< OneD, Array< OneD, NekDouble > > > | m_ntimes_ntimesMFBwd |
| Array< OneD, Array< OneD, NekDouble > > | m_ZimFwd |
| Array< OneD, Array< OneD, NekDouble > > | m_ZimBwd |
| Array< OneD, Array< OneD, NekDouble > > | m_YimFwd |
| Array< OneD, Array< OneD, NekDouble > > | m_YimBwd |
| Array< OneD, Array< OneD, NekDouble > > | m_epsvec |
| Array< OneD, Array< OneD, NekDouble > > | m_muvec |
| Array< OneD, Array< OneD, NekDouble > > | m_negepsvecminus1 |
| Array< OneD, Array< OneD, NekDouble > > | m_negmuvecminus1 |
| Array< OneD, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > > | m_dedxi_cdot_e |
| SpatialDomains::GeomMMF | m_MMFdir |
| Protected Attributes inherited from Nektar::SolverUtils::UnsteadySystem | |
| LibUtilities::TimeIntegrationSchemeSharedPtr | m_intScheme |
| Wrapper to the time integration scheme. More... | |
| LibUtilities::TimeIntegrationSchemeOperators | m_ode |
| The time integration scheme operators to use. More... | |
| Array< OneD, Array< OneD, NekDouble > > | m_previousSolution |
| Storage for previous solution for steady-state check. More... | |
| std::vector< int > | m_intVariables |
| NekDouble | m_cflSafetyFactor |
| CFL safety factor (comprise between 0 to 1). More... | |
| NekDouble | m_CFLGrowth |
| CFL growth rate. More... | |
| NekDouble | m_CFLEnd |
| Maximun cfl in cfl growth. More... | |
| int | m_abortSteps |
| Number of steps between checks for abort conditions. More... | |
| bool | m_explicitDiffusion |
| Indicates if explicit or implicit treatment of diffusion is used. More... | |
| bool | m_explicitAdvection |
| Indicates if explicit or implicit treatment of advection is used. More... | |
| bool | m_explicitReaction |
| Indicates if explicit or implicit treatment of reaction is used. More... | |
| int | m_steadyStateSteps |
| Check for steady state at step interval. More... | |
| NekDouble | m_steadyStateTol |
| Tolerance to which steady state should be evaluated at. More... | |
| int | m_filtersInfosteps |
| Number of time steps between outputting filters information. More... | |
| std::vector< std::pair< std::string, FilterSharedPtr > > | m_filters |
| bool | m_homoInitialFwd |
| Flag to determine if simulation should start in homogeneous forward transformed state. More... | |
| std::ofstream | m_errFile |
| NekDouble | m_epsilon |
| Diffusion coefficient. More... | |
| Protected Attributes inherited from Nektar::SolverUtils::EquationSystem | |
| LibUtilities::CommSharedPtr | m_comm |
| Communicator. More... | |
| bool | m_verbose |
| LibUtilities::SessionReaderSharedPtr | m_session |
| The session reader. More... | |
| std::map< std::string, SolverUtils::SessionFunctionSharedPtr > | m_sessionFunctions |
| Map of known SessionFunctions. More... | |
| LibUtilities::FieldIOSharedPtr | m_fld |
| Field input/output. More... | |
| Array< OneD, MultiRegions::ExpListSharedPtr > | m_fields |
| Array holding all dependent variables. More... | |
| SpatialDomains::BoundaryConditionsSharedPtr | m_boundaryConditions |
| Pointer to boundary conditions object. More... | |
| SpatialDomains::MeshGraphSharedPtr | m_graph |
| Pointer to graph defining mesh. More... | |
| std::string | m_sessionName |
| Name of the session. More... | |
| NekDouble | m_time |
| Current time of simulation. More... | |
| int | m_initialStep |
| Number of the step where the simulation should begin. More... | |
| NekDouble | m_fintime |
| Finish time of the simulation. More... | |
| NekDouble | m_timestep |
| Time step size. More... | |
| NekDouble | m_lambda |
| Lambda constant in real system if one required. More... | |
| NekDouble | m_checktime |
| Time between checkpoints. More... | |
| NekDouble | m_lastCheckTime |
| NekDouble | m_TimeIncrementFactor |
| int | m_nchk |
| Number of checkpoints written so far. More... | |
| int | m_steps |
| Number of steps to take. More... | |
| int | m_checksteps |
| Number of steps between checkpoints. More... | |
| int | m_infosteps |
| Number of time steps between outputting status information. More... | |
| int | m_iterPIT = 0 |
| Number of parallel-in-time time iteration. More... | |
| int | m_windowPIT = 0 |
| Index of windows for parallel-in-time time iteration. More... | |
| int | m_spacedim |
| Spatial dimension (>= expansion dim). More... | |
| int | m_expdim |
| Expansion dimension. More... | |
| bool | m_singleMode |
| Flag to determine if single homogeneous mode is used. More... | |
| bool | m_halfMode |
| Flag to determine if half homogeneous mode is used. More... | |
| bool | m_multipleModes |
| Flag to determine if use multiple homogenenous modes are used. More... | |
| bool | m_useFFT |
| Flag to determine if FFT is used for homogeneous transform. More... | |
| bool | m_homogen_dealiasing |
| Flag to determine if dealiasing is used for homogeneous simulations. More... | |
| bool | m_specHP_dealiasing |
| Flag to determine if dealisising is usde for the Spectral/hp element discretisation. More... | |
| enum MultiRegions::ProjectionType | m_projectionType |
| Type of projection; e.g continuous or discontinuous. More... | |
| Array< OneD, Array< OneD, NekDouble > > | m_traceNormals |
| Array holding trace normals for DG simulations in the forwards direction. More... | |
| Array< OneD, bool > | m_checkIfSystemSingular |
| Flag to indicate if the fields should be checked for singularity. More... | |
| LibUtilities::FieldMetaDataMap | m_fieldMetaDataMap |
| Map to identify relevant solver info to dump in output fields. More... | |
| Array< OneD, NekDouble > | m_movingFrameVelsxyz |
| Moving frame of reference velocities (u, v, w, omega_x, omega_y, omega_z, a_x, a_y, a_z, domega_x, domega_y, domega_z) More... | |
| Array< OneD, NekDouble > | m_movingFrameData |
| Moving frame of reference angles with respect to the. More... | |
| boost::numeric::ublas::matrix< NekDouble > | m_movingFrameProjMat |
| Projection matrix for transformation between inertial and moving. More... | |
| int | m_NumQuadPointsError |
| Number of Quadrature points used to work out the error. More... | |
| enum HomogeneousType | m_HomogeneousType |
| NekDouble | m_LhomX |
| physical length in X direction (if homogeneous) More... | |
| NekDouble | m_LhomY |
| physical length in Y direction (if homogeneous) More... | |
| NekDouble | m_LhomZ |
| physical length in Z direction (if homogeneous) More... | |
| int | m_npointsX |
| number of points in X direction (if homogeneous) More... | |
| int | m_npointsY |
| number of points in Y direction (if homogeneous) More... | |
| int | m_npointsZ |
| number of points in Z direction (if homogeneous) More... | |
| int | m_HomoDirec |
| number of homogenous directions More... | |
Private Member Functions | |
| void | TestSWE2Dproblem (const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield) |
| void | Checkpoint_Output_Cartesian (std::string outname) |
Private Attributes | |
| int | m_RossbyDisturbance |
| int | m_PurturbedJet |
Friends | |
| class | MemoryManager< MMFSWE > |
Additional Inherited Members | |
| Protected Types inherited from Nektar::SolverUtils::EquationSystem | |
| enum | HomogeneousType { eHomogeneous1D , eHomogeneous2D , eHomogeneous3D , eNotHomogeneous } |
| Parameter for homogeneous expansions. More... | |
| Static Protected Attributes inherited from Nektar::SolverUtils::EquationSystem | |
| static std::string | equationSystemTypeLookupIds [] |
| static std::string | projectionTypeLookupIds [] |
Detailed Description
Constructor & Destructor Documentation
◆ ~MMFSWE()
|
override |
Destructor.
Unsteady linear advection equation destructor.
Definition at line 244 of file MMFSWE.cpp.
245{
246}
◆ MMFSWE()
|
protected |
Session reader.
Definition at line 51 of file MMFSWE.cpp.
54{
55 m_planeNumber = 0;
56}
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.
Member Function Documentation
◆ AddCoriolis()
|
protected |
Definition at line 1213 of file MMFSWE.cpp.
1215{
1216 int ncoeffs = outarray[0].size();
1217 int nq = physarray[0].size();
1218
1219 Array<OneD, NekDouble> h(nq);
1220 Array<OneD, NekDouble> tmp(nq);
1221 Array<OneD, NekDouble> tmpc(ncoeffs);
1222
1223 // physarray is primitive
1224 // conservative formulation compute h
1225 // h = \eta + d
1227
1228 int indx = 0;
1230 {
1231 if (j == 0)
1232 {
1233 indx = 2;
1234 }
1235
1236 else if (j == 1)
1237 {
1238 indx = 1;
1239 }
1240
1241 // add to hu equation
1243 Vmath::Vmul(nq, h, 1, tmp, 1, tmp, 1);
1244
1245 if (j == 1)
1246 {
1247 Vmath::Neg(nq, tmp, 1);
1248 }
1249
1250 // N \cdot (e^1 \times e^2 )
1251 // Vmath::Vmul(nq, &m_MF1crossMF2dotSN[0], 1, &tmp[0], 1, &tmp[0], 1);
1252 m_fields[0]->IProductWRTBase(tmp, tmpc);
1253 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1254 }
1255}
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:359
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 549 of file MMFSWE.cpp.
551{
552 // routine works for both primitive and conservative formulations
553 int ncoeffs = outarray[0].size();
554 int nq = physarray[0].size();
555
556 Array<OneD, NekDouble> h(nq);
557 Array<OneD, NekDouble> tmp(nq);
558 Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
560 {
561 fluxvector[i] = Array<OneD, NekDouble>(nq);
562 }
563
564 // Get 0.5 g H*H / || e^m ||^2
565 GetSWEFluxVector(3, physarray, fluxvector);
566
567 Array<OneD, NekDouble> tmpc(ncoeffs);
569 {
571
572 Vmath::Neg(nq, &tmp[0], 1);
573 m_fields[0]->IProductWRTBase(tmp, tmpc);
574
575 Vmath::Vadd(ncoeffs, outarray[j + 1], 1, tmpc, 1, outarray[j + 1], 1);
576 }
577}
void GetSWEFluxVector(const int i, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: MMFSWE.cpp:579
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 1258 of file MMFSWE.cpp.
1260{
1261 int ncoeffs = outarray[0].size();
1262 int nq = physarray[0].size();
1263
1264 Array<OneD, NekDouble> h(nq);
1265 Array<OneD, NekDouble> tmp(nq);
1266 Array<OneD, NekDouble> tmpc(ncoeffs);
1267
1268 // physarray is primitive
1269 // conservative formulation compute h
1270 // h = \eta + d
1272
1274 {
1277
1278 m_fields[0]->IProductWRTBase(tmp, tmpc);
1279
1280 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1281 }
1282}
Array< OneD, Array< OneD, NekDouble > > m_Derivdepth
Definition: MMFSWE.h:87
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.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 1287 of file MMFSWE.cpp.
1289{
1290 // routine works for both primitive and conservative formulations
1291 int ncoeffs = outarray[0].size();
1292 int nq = physarray[0].size();
1293
1294 // Compute h
1295 Array<OneD, NekDouble> h(nq);
1297
1298 Array<OneD, NekDouble> de0dt_cdot_e0;
1299 Array<OneD, NekDouble> de0dt_cdot_e1;
1300 Array<OneD, NekDouble> de1dt_cdot_e0;
1301 Array<OneD, NekDouble> de1dt_cdot_e1;
1302 Compute_demdt_cdot_ek(0, 0, physarray, de0dt_cdot_e0);
1303 Compute_demdt_cdot_ek(1, 0, physarray, de1dt_cdot_e0);
1304 Compute_demdt_cdot_ek(0, 1, physarray, de0dt_cdot_e1);
1305 Compute_demdt_cdot_ek(1, 1, physarray, de1dt_cdot_e1);
1306
1307 Array<OneD, NekDouble> Rott1(nq);
1308 Array<OneD, NekDouble> Rott2(nq);
1309 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e0, 1, Rott1, 1);
1310 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e1, 1, Rott2, 1);
1311 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e0, 1, Rott1, 1, Rott1, 1);
1312 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e1, 1, Rott2, 1, Rott2, 1);
1313
1314 // Multiply H and \partial \phi / \partial t which is assumed to be u_{\phi}
1315 Vmath::Vmul(nq, &h[0], 1, &Rott1[0], 1, &Rott1[0], 1);
1316 Vmath::Vmul(nq, &h[0], 1, &Rott2[0], 1, &Rott2[0], 1);
1317
1318 Vmath::Neg(nq, Rott1, 1);
1319 Vmath::Neg(nq, Rott2, 1);
1320
1321 Array<OneD, NekDouble> tmpc1(ncoeffs);
1322 Array<OneD, NekDouble> tmpc2(ncoeffs);
1323 m_fields[0]->IProductWRTBase(Rott1, tmpc1);
1324 m_fields[0]->IProductWRTBase(Rott2, tmpc2);
1325
1326 Vmath::Vadd(ncoeffs, tmpc1, 1, outarray[1], 1, outarray[1], 1);
1327 Vmath::Vadd(ncoeffs, tmpc2, 1, outarray[2], 1, outarray[2], 1);
1328}
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:1334
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 968 of file MMFSWE.cpp.
972{
973 NekDouble MageF1, MageF2, MageB1, MageB2;
974 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
975 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
976
978 NekDouble uRF, vRF, uLB, vLB;
979
980 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
981 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
982
983 // uRF = uR component in moving frames e^{Fwd}
984 // vRF = vR component in moving frames e^{Fwd}
985 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
986 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
987
988 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
989 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
990 numfluxF[2] = 0.0;
991
992 numfluxF[3] =
993 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
994 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
995 numfluxF[5] = 0.0;
996
997 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
998 numfluxF[7] =
999 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1000 numfluxF[8] = 0.0;
1001
1002 // uLB = uL component in moving frames e^{Bwd}
1003 // vLB = vL component in moving frames e^{Bwd}
1004 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1005 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1006
1007 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1008 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1009 numfluxB[2] = 0.0;
1010
1011 numfluxB[3] =
1012 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1013 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1014 numfluxB[5] = 0.0;
1015
1016 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1017 numfluxB[7] =
1018 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1019 numfluxB[8] = 0.0;
1020}
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:1175
References ComputeMagAndDot(), and m_g.
Referenced by NumericalSWEFlux().
◆ Checkpoint_Output_Cartesian()
|
private |
Definition at line 2774 of file MMFSWE.cpp.
2775{
2778
2779 NekDouble rad_earth = 6.37122 * 1000000;
2780
2781 int nvariables = 7;
2782
2783 // vector u in Cartesian coordinates
2784 std::vector<std::string> variables(nvariables);
2785
2786 variables[0] = "eta";
2787 variables[1] = "hstar";
2788 variables[2] = "vorticity";
2789 variables[3] = "ux";
2790 variables[4] = "uy";
2791 variables[5] = "uz";
2792 variables[6] = "null";
2793
2794 // Obtain \vec{u} in cartesian coordinate
2795 Array<OneD, Array<OneD, NekDouble>> fieldphys(nvariables);
2796 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvariables);
2797 for (int i = 0; i < nvariables; ++i)
2798 {
2799 fieldphys[i] = Array<OneD, NekDouble>(nq, 0.0);
2800 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2801 }
2802
2804 &fieldphys[0][0], 1);
2805 // Vmath::Vadd(nq, &(m_fields[0]->GetPhys())[0], 1, &m_depth[0], 1,
2806 // &fieldphys[1][0], 1);
2807
2809 Vmath::Neg(nq, &fieldphys[1][0], 1);
2811 Vmath::Smul(nq, rad_earth, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2812
2813 Array<OneD, NekDouble> utmp(nq);
2814 Array<OneD, NekDouble> vtmp(nq);
2815
2818
2819 ComputeVorticity(utmp, vtmp, fieldphys[2]);
2820 // u_x = u e^1_x + v e^2_x
2821 int indx = 3;
2823 {
2825 &fieldphys[k + indx][0], 1);
2827 &fieldphys[k + indx][0], 1, &fieldphys[k + indx][0], 1);
2828 }
2829
2830 for (int i = 0; i < nvariables; ++i)
2831 {
2832 m_fields[0]->FwdTrans(fieldphys[i], fieldcoeffs[i]);
2833 }
2834
2836}
void ComputeVorticity(const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
Definition: MMFSWE.cpp:2220
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1238
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 1334 of file MMFSWE.cpp.
1338{
1339 int j, k;
1341
1342 Array<OneD, NekDouble> tmp(nq, 0.0);
1343 Array<OneD, NekDouble> tmpderiv(nq);
1344
1345 outarray = Array<OneD, NekDouble>(nq, 0.0);
1347 {
1349 {
1350 // Compute d e^m / d \xi_1 and d e^m / d \xi_2
1353
1354 Vmath::Vmul(nq, &physarray[j + 1][0], 1, &tmpderiv[0], 1,
1355 &tmpderiv[0], 1);
1356
1358 &outarray[0], 1, &outarray[0], 1);
1359 }
1360 }
1361}
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 2160 of file MMFSWE.cpp.
2163{
2165
2166 Array<OneD, NekDouble> tmp(nq);
2167 Array<OneD, NekDouble> htmp(nq);
2168 Array<OneD, NekDouble> hstmp(nq);
2169
2170 Vmath::Vmul(nq, u, 1, u, 1, tmp, 1);
2171 Vmath::Vvtvp(nq, v, 1, v, 1, tmp, 1, tmp, 1);
2172 Vmath::Vmul(nq, eta, 1, tmp, 1, tmp, 1);
2173
2175 Vmath::Vmul(nq, htmp, 1, htmp, 1, htmp, 1);
2176
2178 Vmath::Vmul(nq, hstmp, 1, hstmp, 1, hstmp, 1);
2179
2180 Vmath::Vsub(nq, htmp, 1, hstmp, 1, htmp, 1);
2182
2183 Vmath::Vadd(nq, htmp, 1, tmp, 1, tmp, 1);
2184 Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2185
2187}
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 2189 of file MMFSWE.cpp.
2192{
2194
2195 Array<OneD, NekDouble> hstartmp(nq);
2196 Array<OneD, NekDouble> tmp(nq);
2197
2199
2200 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2201 Array<OneD, Array<OneD, NekDouble>> Curlu(m_spacedim);
2203 {
2204 uvec[i] = Array<OneD, NekDouble>(nq);
2205 Curlu[i] = Array<OneD, NekDouble>(nq);
2206 }
2207
2208 ComputeVorticity(u, v, tmp);
2209
2211 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
2212 Vmath::Vdiv(nq, tmp, 1, hstartmp, 1, tmp, 1);
2213 Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2214
2216}
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 890 of file MMFSWE.cpp.
894{
896
897 NekDouble hC, huC, hvC, SL, SR, Sstar;
900
901 // Compute SL
902 if (hstar > hL)
903 {
904 SL = uL - cL * sqrt(0.5 * ((hstar * hstar + hstar * hL) / (hL * hL)));
905 }
906 else
907 {
908 SL = uL - cL;
909 }
910
911 // Compute SR
912 if (hstar > hR)
913 {
914 SR = uR + cR * sqrt(0.5 * ((hstar * hstar + hstar * hR) / (hR * hR)));
915 }
916 else
917 {
918 SR = uR + cR;
919 }
920
921 if (fabs(hR * (uR - SR) - hL * (uL - SL)) <= 1.0e-15)
922 {
923 Sstar = 0.0;
924 }
925 else
926 {
927 Sstar = (SL * hR * (uR - SR) - SR * hL * (uL - SL)) /
928 (hR * (uR - SR) - hL * (uL - SL));
929 }
930
931 if (SL >= 0)
932 {
933 hflux = hL * uL;
934 huflux = uL * uL * hL + 0.5 * g * hL * hL;
935 hvflux = hL * uL * vL;
936 }
937 else if (SR <= 0)
938 {
939 hflux = hR * uR;
940 huflux = uR * uR * hR + 0.5 * g * hR * hR;
941 hvflux = hR * uR * vR;
942 }
943 else
944 {
945 if ((SL < 0) && (Sstar >= 0))
946 {
947 hC = hL * ((SL - uL) / (SL - Sstar));
948 huC = hC * Sstar;
949 hvC = hC * vL;
950
951 hflux = hL * uL + SL * (hC - hL);
952 huflux = (uL * uL * hL + 0.5 * g * hL * hL) + SL * (huC - hL * uL);
953 hvflux = (uL * vL * hL) + SL * (hvC - hL * vL);
954 }
955 else
956 {
957 hC = hR * ((SR - uR) / (SR - Sstar));
958 huC = hC * Sstar;
959 hvC = hC * vR;
960
961 hflux = hR * uR + SR * (hC - hR);
962 huflux = (uR * uR * hR + 0.5 * g * hR * hR) + SR * (huC - hR * uR);
963 hvflux = (uR * vR * hR) + SR * (hvC - hR * vR);
964 }
965 }
966}
References m_g, and tinysimd::sqrt().
Referenced by RiemannSolverHLLC().
◆ ComputeMagAndDot()
|
protected |
Definition at line 1175 of file MMFSWE.cpp.
1180{
1181 NekDouble MF1x, MF1y, MF1z, MF2x, MF2y, MF2z;
1182 NekDouble MB1x, MB1y, MB1z, MB2x, MB2y, MB2z;
1183
1184 MF1x = m_MFtraceFwd[0][0][index];
1185 MF1y = m_MFtraceFwd[0][1][index];
1186 MF1z = m_MFtraceFwd[0][2][index];
1187
1188 MF2x = m_MFtraceFwd[1][0][index];
1189 MF2y = m_MFtraceFwd[1][1][index];
1190 MF2z = m_MFtraceFwd[1][2][index];
1191
1192 MB1x = m_MFtraceBwd[0][0][index];
1193 MB1y = m_MFtraceBwd[0][1][index];
1194 MB1z = m_MFtraceBwd[0][2][index];
1195
1196 MB2x = m_MFtraceBwd[1][0][index];
1197 MB2y = m_MFtraceBwd[1][1][index];
1198 MB2z = m_MFtraceBwd[1][2][index];
1199
1200 // MFtrace = MFtrace [ j*spacedim + k ], j = shape, k = sapce
1201 MageF1 = MF1x * MF1x + MF1y * MF1y + MF1z * MF1z;
1202 MageF2 = MF2x * MF2x + MF2y * MF2y + MF2z * MF2z;
1203 MageB1 = MB1x * MB1x + MB1y * MB1y + MB1z * MB1z;
1204 MageB2 = MB2x * MB2x + MB2y * MB2y + MB2z * MB2z;
1205
1206 eF1_cdot_eB1 = MF1x * MB1x + MF1y * MB1y + MF1z * MB1z;
1207 eF1_cdot_eB2 = MF1x * MB2x + MF1y * MB2y + MF1z * MB2z;
1208 eF2_cdot_eB1 = MF2x * MB1x + MF2y * MB1y + MF2z * MB1z;
1209 eF2_cdot_eB2 = MF2x * MB2x + MF2y * MB2y + MF2z * MB2z;
1210}
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 2150 of file MMFSWE.cpp.
2151{
2153
2154 Array<OneD, NekDouble> tmp(nq);
2156
2158}
References m_depth, Nektar::SolverUtils::EquationSystem::m_fields, and Vmath::Vadd().
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ ComputeNablaCdotVelocity()
Definition at line 2241 of file MMFSWE.cpp.
2242{
2244
2245 Array<OneD, NekDouble> velcoeff(nq, 0.0);
2246
2247 Array<OneD, NekDouble> Dtmp0(nq);
2248 Array<OneD, NekDouble> Dtmp1(nq);
2249 Array<OneD, NekDouble> Dtmp2(nq);
2250 Array<OneD, NekDouble> Drv(nq);
2251
2252 vellc = Array<OneD, NekDouble>(nq, 0.0);
2253
2254 // m_vellc = \nabla m_vel \cdot tan_i
2255 Array<OneD, NekDouble> tmp(nq);
2256 Array<OneD, NekDouble> vessel(nq);
2257
2259 {
2260 Vmath::Zero(nq, velcoeff, 1);
2262 {
2263 // a_j = tan_j cdot m_vel
2265 1, &velcoeff[0], 1, &velcoeff[0], 1);
2266 }
2267
2268 // d a_j / d x^k
2269 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
2270
2272 {
2273 // tan_j^k ( d a_j / d x^k )
2274 switch (k)
2275 {
2276 case (0):
2277 {
2279 1, &vellc[0], 1, &vellc[0], 1);
2280 }
2281 break;
2282
2283 case (1):
2284 {
2286 1, &vellc[0], 1, &vellc[0], 1);
2287 }
2288 break;
2289
2290 case (2):
2291 {
2293 1, &vellc[0], 1, &vellc[0], 1);
2294 }
2295 break;
2296 }
2297 }
2298 }
2299}
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 2744 of file MMFSWE.cpp.
2745{
2746 NekDouble uphi, f, dh;
2747
2748 uphi = ComputeUnstableJetuphi(theta);
2749 f = 2.0 * m_Omega * sin(theta);
2750
2751 dh = f * uphi + tan(theta) * uphi * uphi;
2752
2753 return dh;
2754}
NekDouble ComputeUnstableJetuphi(const NekDouble theta)
Definition: MMFSWE.cpp:2756
References ComputeUnstableJetuphi(), and m_Omega.
Referenced by UnstableJetFlow().
◆ ComputeUnstableJetuphi()
Definition at line 2756 of file MMFSWE.cpp.
2757{
2758 NekDouble uphi;
2759
2761 {
2764 }
2765
2766 else
2767 {
2768 uphi = 0.0;
2769 }
2770
2771 return uphi;
2772}
References m_en, m_theta0, m_theta1, and m_uthetamax.
Referenced by ComputeUnstableJetEta(), and UnstableJetFlow().
◆ ComputeVorticity()
|
protected |
Definition at line 2220 of file MMFSWE.cpp.
2223{
2225
2226 Array<OneD, NekDouble> tmp(nq);
2227
2228 Vorticity = Array<OneD, NekDouble>(nq, 0.0);
2229
2232 &Vorticity[0], 1);
2233
2235 Vmath::Neg(nq, tmp, 1);
2237
2238 Vmath::Vadd(nq, tmp, 1, Vorticity, 1, Vorticity, 1);
2239}
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 2919 of file MMFSWE.cpp.
2920{
2922
2923 // u = hu/h
2925 m_fields[1]->UpdatePhys(), 1);
2926
2927 // v = hv/ v
2929 m_fields[2]->UpdatePhys(), 1);
2930
2931 // \eta = h - d
2933 m_fields[0]->UpdatePhys(), 1);
2934}
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:764
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 2839 of file MMFSWE.cpp.
2842{
2844
2845 if (physin[0].get() == physout[0].get())
2846 {
2847 // copy indata and work with tmp array
2848 Array<OneD, Array<OneD, NekDouble>> tmp(3);
2849 for (int i = 0; i < 3; ++i)
2850 {
2851 // deep copy
2852 tmp[i] = Array<OneD, NekDouble>(nq);
2853 Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2854 }
2855
2856 // \eta = h - d
2858
2859 // u = hu/h
2860 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
2861
2862 // v = hv/h
2863 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
2864 }
2865 else
2866 {
2867 // \eta = h - d
2869
2870 // u = hu/h
2871 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
2872
2873 // v = hv/h
2874 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
2875 }
2876}
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vsub().
◆ create()
|
inlinestatic |
Creates an instance of this class.
Definition at line 67 of file MMFSWE.h.
70 {
72 MemoryManager<MMFSWE>::AllocateSharedPtr(pSession, pGraph);
73 p->InitObject();
75 }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp: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 1370 of file MMFSWE.cpp.
1373{
1375 {
1377 {
1378 ConservativeToPrimitive(inarray, outarray);
1379 SetBoundaryConditions(outarray, time);
1380 PrimitiveToConservative(outarray, outarray);
1381 }
1382 break;
1383 default:
1385 break;
1386 }
1387}
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &inarray, NekDouble time)
Definition: MMFSWE.cpp:1389
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
Definition: EquationSystem.h:418
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 425 of file MMFSWE.cpp.
428{
429 int i;
430 int nvariables = inarray.size();
433
434 // inarray in physical space
435 Array<OneD, Array<OneD, NekDouble>> physarray(nvariables);
436 Array<OneD, Array<OneD, NekDouble>> modarray(nvariables);
437
438 for (i = 0; i < nvariables; ++i)
439 {
440 physarray[i] = Array<OneD, NekDouble>(nq);
441 modarray[i] = Array<OneD, NekDouble>(ncoeffs);
442 }
443
444 // (h, hu, hv) -> (\eta, u, v)
445 ConservativeToPrimitive(inarray, physarray);
446
447 // Weak Directional Derivative
448 WeakDGSWEDirDeriv(physarray, modarray);
449 AddDivForGradient(physarray, modarray);
450
451 for (i = 0; i < nvariables; ++i)
452 {
453 Vmath::Neg(ncoeffs, modarray[i], 1);
454 }
455
456 // coriolis forcing
458 {
459 AddCoriolis(physarray, modarray);
460 }
461
462 // Bottom Elevation Effect
463 // Add g H \nabla m_depth
464 AddElevationEffect(physarray, modarray);
465
466 // Add terms concerning the rotation of the moving frame
468 {
469 AddRotation(physarray, modarray);
470 }
471
472 for (i = 0; i < nvariables; ++i)
473 {
474 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
475 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
476 }
477}
void AddDivForGradient(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:549
void AddCoriolis(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1213
void WeakDGSWEDirDeriv(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField)
Definition: MMFSWE.cpp:479
void AddRotation(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1287
void AddElevationEffect(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1258
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:706
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 1709 of file MMFSWE.cpp.
1710{
1712 {
1714 {
1716 }
1717 break;
1718
1720 {
1722 }
1723 break;
1724
1729 {
1731 }
1732 break;
1733
1734 default:
1735 break;
1736 }
1737}
void EvaluateStandardCoriolis(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1771
void EvaluateCoriolisForZonalFlow(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1739
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 1739 of file MMFSWE.cpp.
1740{
1742 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1743
1744 Array<OneD, NekDouble> x(nq);
1745 Array<OneD, NekDouble> y(nq);
1746 Array<OneD, NekDouble> z(nq);
1747
1749
1750 NekDouble x0j, x1j, x2j;
1751 NekDouble tmp;
1752
1753 outarray = Array<OneD, NekDouble>(nq, 0.0);
1754 for (int j = 0; j < nq; ++j)
1755 {
1756 x0j = x[j];
1757 x1j = y[j];
1758 x2j = z[j];
1759
1760 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1761 cos_theta);
1762
1763 // H = 2 \Omega *(- \cos \phi \cos \theta \sin \alpha + \sin \theta \cos
1764 // \alpha )
1765 tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
1766 sin_theta * cos(m_alpha);
1767 outarray[j] = 2.0 * m_Omega * tmp;
1768 }
1769}
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 1771 of file MMFSWE.cpp.
1772{
1774
1775 NekDouble x0j, x1j, x2j;
1776 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1777
1778 Array<OneD, NekDouble> x(nq);
1779 Array<OneD, NekDouble> y(nq);
1780 Array<OneD, NekDouble> z(nq);
1781
1783
1784 outarray = Array<OneD, NekDouble>(nq, 0.0);
1785 for (int j = 0; j < nq; ++j)
1786 {
1787 x0j = x[j];
1788 x1j = y[j];
1789 x2j = z[j];
1790
1791 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1792 cos_theta);
1793
1794 outarray[j] = 2.0 * m_Omega * sin_theta;
1795 }
1796}
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 1569 of file MMFSWE.cpp.
1570{
1572
1573 m_depth = Array<OneD, NekDouble>(nq, 0.0);
1574
1576 {
1578 {
1580 }
1581 break;
1582
1584 {
1585 // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1586 Array<OneD, NekDouble> x(nq);
1587 Array<OneD, NekDouble> y(nq);
1588 Array<OneD, NekDouble> z(nq);
1589
1591
1592 NekDouble x0j, x1j, x2j;
1593 NekDouble sin_varphi, cos_varphi, sin_theta, cos_theta;
1594
1595 for (int j = 0; j < nq; ++j)
1596 {
1597 x0j = x[j];
1598 x1j = y[j];
1599 x2j = z[j];
1600
1601 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1602 sin_theta, cos_theta);
1603 m_depth[j] =
1606 }
1607 }
1608 break;
1609
1611 {
1612 // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1613 Array<OneD, NekDouble> x(nq);
1614 Array<OneD, NekDouble> y(nq);
1615 Array<OneD, NekDouble> z(nq);
1616
1618
1619 int indx = 0;
1620 NekDouble x0j, x1j, x2j, dist2;
1621 NekDouble phi, theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
1622 NekDouble hRad, phic, thetac;
1623 NekDouble Tol = 0.000001;
1624
1625 hRad = m_pi / 9.0;
1626 phic = -m_pi / 2.0;
1627 thetac = m_pi / 6.0;
1628
1629 for (int j = 0; j < nq; ++j)
1630 {
1631 x0j = x[j];
1632 x1j = y[j];
1633 x2j = z[j];
1634
1635 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1636 sin_theta, cos_theta);
1637
1639 std::abs(sin(thetac) - sin_theta)) < Tol)
1640 {
1641 std::cout << "A point " << j
1642 << " is coincient with the singularity "
1643 << std::endl;
1644 indx = 1;
1645 }
1646
1647 phi = atan2(sin_varphi, cos_varphi);
1648 theta = atan2(sin_theta, cos_theta);
1649
1650 // Compute r
1651 dist2 = (phi - phic) * (phi - phic) +
1652 (theta - thetac) * (theta - thetac);
1653
1654 if (dist2 > hRad * hRad)
1655 {
1656 dist2 = hRad * hRad;
1657 }
1658
1660 }
1661
1662 if (!indx)
1663 {
1664 std::cout << "No point is coincident with the singularity point"
1665 << std::endl;
1666 }
1667 }
1668 break;
1669
1671 {
1672 for (int j = 0; j < nq; ++j)
1673 {
1675 }
1676 }
1677 break;
1678
1681 {
1683 }
1684 break;
1685
1686 default:
1687 {
1689 }
1690 break;
1691 }
1692
1693 // Comptue \nabla m_depth \cdot e^i
1695
1697 {
1698 m_Derivdepth[j] = Array<OneD, NekDouble>(nq);
1700 m_Derivdepth[j]);
1701 }
1702
1703 std::cout << "Water Depth (m_depth) was generated with mag = "
1707}
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 579 of file MMFSWE.cpp.
582{
584
585 switch (i)
586 {
587 // flux function for the h equation = [(\eta + d ) u, (\eta + d) v ]
588 case 0:
589 {
590 // h in flux 1
592
593 // hu in flux 0
594 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
595
596 // hv in flux 1
597 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
598 }
599 break;
600
601 // flux function for the hu equation = [Hu^2 + 0.5 g H^2, Huv]
602 case 1:
603 {
604 Array<OneD, NekDouble> tmp(nq);
605
606 // h in tmp
608
609 // hu in flux 1
610 Vmath::Vmul(nq, tmp, 1, physfield[1], 1, flux[1], 1);
611
612 // huu in flux 0
613 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
614
615 // hh in tmp
616 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
617
618 // huu + 0.5 g hh in flux 0
619 // Daxpy overwrites flux[0] on exit
621
622 // huv in flux 1
623 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
624 }
625 break;
626
627 // flux function for the hv equation = [Huv, Hv^2]
628 case 2:
629 {
630 Array<OneD, NekDouble> tmp(nq);
631
632 // h in tmp
634
635 // hv in flux 0
636 Vmath::Vmul(nq, tmp, 1, physfield[2], 1, flux[0], 1);
637
638 // hvv in flux 1
639 Vmath::Vmul(nq, flux[0], 1, physfield[2], 1, flux[1], 1);
640
641 // huv in flux 0
642 Vmath::Vmul(nq, flux[0], 1, physfield[1], 1, flux[0], 1);
643
644 // hh in tmp
645 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
646
647 // hvv + 0.5 g hh in flux 1
649 }
650 break;
651
652 // flux function 0.5 g h * h
653 case 3:
654 {
655 Array<OneD, NekDouble> h(nq);
656 flux[0] = Array<OneD, NekDouble>(nq, 0.0);
657
658 // h in tmp
660
661 // hh in tmp
662 Vmath::Vmul(nq, h, 1, h, 1, h, 1);
663
664 // 0.5 g hh in flux 0
666 }
667 break;
668
669 default:
670 break;
671 }
672}
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 2402 of file MMFSWE.cpp.
2405{
2407
2408 NekDouble uhat, vhat;
2409 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2410 NekDouble x0j, x1j, x2j;
2411
2412 Array<OneD, NekDouble> eta(nq, 0.0);
2413 Array<OneD, NekDouble> u(nq, 0.0);
2414 Array<OneD, NekDouble> v(nq, 0.0);
2415
2416 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2418 {
2419 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2420 }
2421
2422 Array<OneD, NekDouble> x(nq);
2423 Array<OneD, NekDouble> y(nq);
2424 Array<OneD, NekDouble> z(nq);
2425
2427
2428 for (int j = 0; j < nq; ++j)
2429 {
2430 x0j = x[j];
2431 x1j = y[j];
2432 x2j = z[j];
2433
2434 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2435 cos_theta);
2436
2437 // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2439 sin_theta * sin_theta;
2440
2441 uhat = m_u0 * cos_theta;
2442 vhat = 0.0;
2443
2444 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2445 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2446 uvec[2][j] = vhat * cos_theta;
2447 }
2448
2449 // Projection of u onto the tangent plane with conserving the mag. of the
2450 // velocity.
2451 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2452
2454 {
2455 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2456 }
2457
2458 // u is projected on the tangent plane with preserving its length
2459 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2460
2461 // Change it to the coordinate of moving frames
2462 // CartesianToMovingframes(0,uvecproj,u);
2463 // CartesianToMovingframes(1,uvecproj,v);
2464
2465 u = CartesianToMovingframes(uvec, 0);
2466 v = CartesianToMovingframes(uvec, 1);
2467
2468 switch (field)
2469 {
2470 case (0):
2471 {
2472 outfield = eta;
2473 }
2474 break;
2475
2476 case (1):
2477 {
2478 outfield = u;
2479 }
2480 break;
2481
2482 case (2):
2483 {
2484 outfield = v;
2485 }
2486 break;
2487 }
2488}
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(), 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 1022 of file MMFSWE.cpp.
1026{
1027 int nvariables = 3;
1028 NekDouble MageF1, MageF2, MageB1, MageB2;
1029 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1030 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1031
1033 NekDouble uRF, vRF, uLB, vLB;
1034 NekDouble velL, velR, lambdaF, lambdaB;
1035
1036 Array<OneD, NekDouble> EigF(nvariables);
1037 Array<OneD, NekDouble> EigB(nvariables);
1038
1039 // Compute Magnitude and Dot product of moving frames for the index
1040 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1041 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1042
1043 // Get the velocity in the normal to the edge
1046
1047 EigF[0] = velL - sqrt(g * hL);
1048 EigF[1] = velL;
1049 EigF[2] = velL + sqrt(g * hL);
1050
1051 EigB[0] = velR - sqrt(g * hR);
1052 EigB[1] = velR;
1053 EigB[2] = velR + sqrt(g * hR);
1054
1055 lambdaF = Vmath::Vamax(nvariables, EigF, 1);
1056 lambdaB = Vmath::Vamax(nvariables, EigB, 1);
1057
1058 // uRF = uR component in moving frames e^{Fwd}
1059 // vRF = vR component in moving frames e^{Fwd}
1060 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1061 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1062
1063 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1064 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1065 numfluxF[2] = 0.5 * lambdaF * (hL - hR);
1066
1067 numfluxF[3] = 0.5 * (hL * uL * uL * MageF1 + hR * uRF * uRF * MageB1 +
1068 0.5 * g * (hL * hL + hR * hR));
1069 numfluxF[4] = 0.5 * (hL * uL * vL * MageF1 + hR * uRF * vRF * MageB1);
1070 numfluxF[5] = 0.5 * lambdaF * (uL * hL - uRF * hR);
1071
1072 numfluxF[6] = 0.5 * (hL * uL * vL * MageF2 + hR * uRF * vRF * MageB2);
1073 numfluxF[7] = 0.5 * (hL * vL * vL * MageF2 + hR * vRF * vRF * MageB2 +
1074 0.5 * g * (hL * hL + hR * hR));
1075 numfluxF[8] = 0.5 * lambdaF * (vL * hL - vRF * hR);
1076
1077 // uLB = uL component in moving frames e^{Bwd}
1078 // vLB = vL component in moving frames e^{Bwd}
1079 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1080 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1081
1082 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1083 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1084 numfluxB[2] = 0.5 * lambdaB * (hL - hR);
1085
1086 numfluxB[3] = 0.5 * (hR * uR * uR * MageB1 + hR * uLB * uLB * MageF1 +
1087 0.5 * g * (hR * hR + hL * hL));
1088 numfluxB[4] = 0.5 * (hR * uR * vR * MageB1 + hR * uLB * vLB * MageF1);
1089 numfluxB[5] = 0.5 * lambdaB * (uLB * hL - uR * hR);
1090
1091 numfluxB[6] = 0.5 * (hR * uR * vR * MageB2 + hR * uLB * vLB * MageF2);
1092 numfluxB[7] = 0.5 * (hR * vR * vR * MageB2 + hR * vLB * vLB * MageF2 +
1093 0.5 * g * (hR * hR + hL * hL));
1094 numfluxB[8] = 0.5 * lambdaB * (vLB * hL - vR * hR);
1095}
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 674 of file MMFSWE.cpp.
677{
678
679 int i, k;
681 int nvariables = 3; // only the dependent variables
682
683 // get temporary arrays
684 Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
685 Array<OneD, Array<OneD, NekDouble>> Bwd(nvariables);
686 Array<OneD, NekDouble> DepthFwd(nTraceNumPoints);
687 Array<OneD, NekDouble> DepthBwd(nTraceNumPoints);
688
689 for (i = 0; i < nvariables; ++i)
690 {
691 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
692 Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
693 }
694
695 // get the physical values at the trace
696 for (i = 0; i < nvariables; ++i)
697 {
698 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
699
700 // Copy Fwd to Bwd at boundaries
702 }
705
706 // note that we are using the same depth - i.e. the depth is assumed
707 // continuous...
709 {
711 {
712 // rotate the values to the normal direction
713 NekDouble tmpX, tmpY;
714
715 // Fwd[1] = hu^+ = ( h^1^+ (e^1 \cdot n) + h^2^+ ( e^2 \cdot n) )
716 // Fwd[2] = hv^+ = ( h^1^+ (e^1 \cdot n^{\perp}) + h^2^+ ( e^2 \cdot
717 // n^{\perp}) )
718 // Bwd[1] = hu^- = ( h^1^- (e^1 \cdot n) + h^2^- ( e^2 \cdot n) )
719 // Bwd[2] = hv^- = ( h^1^- (e^1 \cdot n^{\perp}) + h^2^- ( e^2 \cdot
720 // n^{\perp}) )
721
722 for (k = 0; k < nTraceNumPoints; ++k)
723 {
724 tmpX = Fwd[1][k] * m_ncdotMFFwd[0][k] +
725 Fwd[2][k] * m_ncdotMFFwd[1][k];
726 tmpY = Fwd[1][k] * m_nperpcdotMFFwd[0][k] +
727 Fwd[2][k] * m_nperpcdotMFFwd[1][k];
728 Fwd[1][k] = tmpX;
729 Fwd[2][k] = tmpY;
730
731 tmpX = Bwd[1][k] * m_ncdotMFBwd[0][k] +
732 Bwd[2][k] * m_ncdotMFBwd[1][k];
733 tmpY = Bwd[1][k] * m_nperpcdotMFBwd[0][k] +
734 Bwd[2][k] * m_nperpcdotMFBwd[1][k];
735 Bwd[1][k] = tmpX;
736 Bwd[2][k] = tmpY;
737 }
738
739 // Solve the Riemann problem
740 NekDouble denomFwd, denomBwd;
741 Array<OneD, NekDouble> numfluxF(nvariables);
742 Array<OneD, NekDouble> numfluxB(nvariables);
743
744 NekDouble eF1n, eF2n, eB1n, eB2n;
745 NekDouble eF1t, eF2t, eB1t, eB2t;
746 for (k = 0; k < nTraceNumPoints; ++k)
747 {
748 RiemannSolverHLLC(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
749 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
750 Bwd[2][k], numfluxF, numfluxB);
751
752 // uflux = 1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
753 // \cdot n)(e^1 \cdot n^{\perp} ) ] ( e^2 \cdot n^{\perp} huflux
754 // - e^2 \cdot n hvflux
755 // vflux = -1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
756 // \cdot n)(e^1 \cdot n^{\perp} ) ] ( -e^1 \cdot n^{\perp}
757 // huflux + e^1 \cdot n hvflux
758 eF1n = m_ncdotMFFwd[0][k];
759 eF2n = m_ncdotMFFwd[1][k];
760 eB1n = m_ncdotMFBwd[0][k];
761 eB2n = m_ncdotMFBwd[1][k];
762
763 eF1t = m_nperpcdotMFFwd[0][k];
764 eF2t = m_nperpcdotMFFwd[1][k];
765 eB1t = m_nperpcdotMFBwd[0][k];
766 eB2t = m_nperpcdotMFBwd[1][k];
767
768 denomFwd = eF1n * eF2t - eF2n * eF1t;
769 denomBwd = eB1n * eB2t - eB2n * eB1t;
770
771 numfluxFwd[0][k] = numfluxF[0];
772 numfluxFwd[1][k] = (1.0 / denomFwd) *
773 (eF2t * numfluxF[1] - eF2n * numfluxF[2]);
774 numfluxFwd[2][k] =
775 (1.0 / denomFwd) *
776 (-1.0 * eF1t * numfluxF[1] + eF1n * numfluxF[2]);
777
778 numfluxBwd[0][k] = 1.0 * numfluxB[0];
779 numfluxBwd[1][k] = (1.0 / denomBwd) *
780 (eB2t * numfluxB[1] - eB2n * numfluxB[2]);
781 numfluxBwd[2][k] =
782 (1.0 / denomBwd) *
783 (-1.0 * eB1t * numfluxB[1] + eB1n * numfluxB[2]);
784 }
785 }
786 break;
787
791 {
792 Array<OneD, NekDouble> numfluxF(nvariables * (m_shapedim + 1));
793 Array<OneD, NekDouble> numfluxB(nvariables * (m_shapedim + 1));
794
795 for (k = 0; k < nTraceNumPoints; ++k)
796 {
798 {
799 AverageFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
800 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
801 Bwd[2][k], numfluxF, numfluxB);
802 }
803
805 {
806 LaxFriedrichFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
807 Fwd[2][k], Bwd[0][k] + DepthFwd[k],
808 Bwd[1][k], Bwd[2][k], numfluxF, numfluxB);
809 }
810
812 {
813 RusanovFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
814 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
815 Bwd[2][k], numfluxF, numfluxB);
816 }
817
818 int indx;
819 NekDouble tmpF0, tmpF1, tmpB0, tmpB1;
820 for (i = 0; i < nvariables; ++i)
821 {
822 indx = i * (m_shapedim + 1);
823
824 tmpF0 = numfluxF[indx] * m_ncdotMFFwd[0][k];
825 tmpF1 = numfluxF[indx + 1] * m_ncdotMFFwd[1][k];
826
827 tmpB0 = numfluxB[indx] * m_ncdotMFBwd[0][k];
828 tmpB1 = numfluxB[indx + 1] * m_ncdotMFBwd[1][k];
829
830 numfluxFwd[i][k] = tmpF0 + tmpF1 + numfluxF[indx + 2];
831 numfluxBwd[i][k] = tmpB0 + tmpB1 + numfluxB[indx + 2];
832 }
833 }
834 }
835 break;
836
837 default:
838 {
840 }
841 break;
842 }
843}
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:968
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:845
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:1022
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:1097
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:739
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 2936 of file MMFSWE.cpp.
2937{
2939
2940 // h = \eta + d
2942 m_fields[0]->UpdatePhys(), 1);
2943
2944 // hu = h * u
2946 m_fields[1]->UpdatePhys(), 1);
2947
2948 // hv = h * v
2950 m_fields[2]->UpdatePhys(), 1);
2951}
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 2878 of file MMFSWE.cpp.
2881{
2882
2884
2885 if (physin[0].get() == physout[0].get())
2886 {
2887 // copy indata and work with tmp array
2888 Array<OneD, Array<OneD, NekDouble>> tmp(3);
2889 for (int i = 0; i < 3; ++i)
2890 {
2891 // deep copy
2892 tmp[i] = Array<OneD, NekDouble>(nq);
2893 Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2894 }
2895
2896 // h = \eta + d
2898
2899 // hu = h * u
2900 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
2901
2902 // hv = h * v
2903 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
2904 }
2905 else
2906 {
2907 // h = \eta + d
2909
2910 // hu = h * u
2911 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
2912
2913 // hv = h * v
2914 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
2915 }
2916}
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vmul().
◆ RiemannSolverHLLC()
|
protected |
Definition at line 845 of file MMFSWE.cpp.
850{
852
855
856 NekDouble uRF, vRF, uLB, vLB;
857 NekDouble hstarF, hstarB;
858
859 // Temporary assignments
860 uRF = uR;
861 vRF = vR;
862 uLB = uL;
863 vLB = vL;
864
865 // the two-rarefaction wave assumption
866 hstarF = 0.5 * (cL + cR) + 0.25 * (uL - uRF);
867 hstarF *= hstarF;
868 hstarF *= (1.0 / g);
869
870 hstarB = 0.5 * (cL + cR) + 0.25 * (uLB - uR);
871 hstarB *= hstarB;
872 hstarB *= (1.0 / g);
873
874 NekDouble hfluxF, hufluxF, hvfluxF;
875 NekDouble hfluxB, hufluxB, hvfluxB;
876 Computehhuhvflux(hL, uL, vL, hR, uRF, vRF, hstarF, hfluxF, hufluxF,
877 hvfluxF);
878 Computehhuhvflux(hL, uLB, vLB, hR, uR, vR, hstarB, hfluxB, hufluxB,
879 hvfluxB);
880
881 numfluxF[0] = hfluxF;
882 numfluxF[1] = hufluxF;
883 numfluxF[2] = hvfluxF;
884
885 numfluxB[0] = hfluxB;
886 numfluxB[1] = hufluxB;
887 numfluxB[2] = hvfluxB;
888}
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:890
References Computehhuhvflux(), m_g, and tinysimd::sqrt().
Referenced by NumericalSWEFlux().
◆ RossbyWave()
|
protected |
Definition at line 2599 of file MMFSWE.cpp.
2600{
2602 NekDouble uhat, vhat;
2603 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2604 NekDouble x0j, x1j, x2j;
2605 NekDouble Ath, Bth, Cth, tmp;
2606
2607 Array<OneD, NekDouble> eta(nq, 0.0);
2608 Array<OneD, NekDouble> u(nq, 0.0);
2609 Array<OneD, NekDouble> v(nq, 0.0);
2610
2611 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2613 {
2614 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2615 }
2616
2617 Array<OneD, NekDouble> x(nq);
2618 Array<OneD, NekDouble> y(nq);
2619 Array<OneD, NekDouble> z(nq);
2620
2622
2623 NekDouble R = 4.0;
2624 NekDouble cos2theta, cosRtheta, cos6theta, cos2Rtheta, cosRm1theta;
2625 NekDouble cos2phi, cos4phi, sin4phi, cos8phi;
2626
2627 // disturbancees of Rossby-Haurwitz Wave
2628 NekDouble x0d, y0d, z0d, phi0, theta0;
2629 // NekDouble rad_earth = 6.37122 * 1000000;
2630
2631 phi0 = 40.0 * m_pi / 180.0;
2632 theta0 = 50.0 * m_pi / 180.0;
2633
2634 x0d = cos(phi0) * cos(theta0);
2635 y0d = sin(phi0) * cos(theta0);
2636 z0d = sin(theta0);
2637
2638 for (int j = 0; j < nq; ++j)
2639 {
2640 x0j = x[j];
2641 x1j = y[j];
2642 x2j = z[j];
2643
2644 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2645 cos_theta);
2646
2647 // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2648 // \sin \alpha + \sin \theta \cos \alpha )^2
2649
2650 // tmp = cos^{2R} \theta, R = 4;
2651 cos2theta = cos_theta * cos_theta;
2652 cosRm1theta = cos_theta * cos2theta;
2653 cosRtheta = cos2theta * cos2theta;
2654 cos6theta = cos2theta * cosRtheta;
2655 cos2Rtheta = cosRtheta * cosRtheta;
2656
2658 Ath = tmp * cos2theta +
2660 ((R + 1.0) * cosRtheta + (2 * R * R - R - 2.0) * cos2theta -
2661 2 * R * R);
2662
2664 Bth = tmp * cosRtheta *
2665 ((R * R + 2 * R + 2) - (R + 1.0) * (R + 1.0) * cos2theta);
2666
2667 Cth =
2669
2670 // cos (2 \phi) = 2 * cos \phi * cos \phi - 1.0
2671 cos2phi = 2.0 * cos_varphi * cos_varphi - 1.0;
2672 cos4phi = 2.0 * cos2phi * cos2phi - 1.0;
2673 cos8phi = 2.0 * cos4phi * cos4phi - 1.0;
2674
2675 // sin (2 \phi) = 2 * cos \phi * sin \phi
2676 sin4phi = 4.0 * sin_varphi * cos_varphi * cos2phi;
2677
2679
2680 // disturbances is added
2682 {
2683 eta[j] = eta[j] *
2684 (1.0 + (1.0 / 40.0) * (x0j * x0d + x1j * y0d + x2j * z0d));
2685 }
2686
2687 // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2688 // || e^2 ||^2
2689 uhat = m_angfreq * cos_theta +
2690 m_K * cosRm1theta * (R * sin_theta * sin_theta - cos2theta) *
2691 cos4phi;
2692 vhat = -1.0 * m_K * R * cosRm1theta * sin_theta * sin4phi;
2693
2694 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2695 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2696 uvec[2][j] = vhat * cos_theta;
2697 }
2698
2699 // NekDouble etamin, etaaver;
2700 // etamin = Vmath::Vmin(nq, eta, 1)*rad_earth;
2701 // etaaver = Vmath::Vsum(nq, eta, 1)/nq*rad_earth;
2702
2703 // Projection of u onto the tangent plane with conserving the mag. of the
2704 // velocity.
2705 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2706
2708 {
2709 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2710 }
2711
2712 // u is projected on the tangent plane with preserving its length
2713 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2714
2715 // Change it to the coordinate of moving frames
2716 // CartesianToMovingframes(0,uvecproj,u);
2717 // CartesianToMovingframes(1,uvecproj,v);
2718
2719 u = CartesianToMovingframes(uvec, 0);
2720 v = CartesianToMovingframes(uvec, 1);
2721
2722 switch (field)
2723 {
2724 case (0):
2725 {
2726 outfield = eta;
2727 }
2728 break;
2729
2730 case (1):
2731 {
2732 outfield = u;
2733 }
2734 break;
2735
2736 case (2):
2737 {
2738 outfield = v;
2739 }
2740 break;
2741 }
2742}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_angfreq, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, m_K, m_Omega, Nektar::SolverUtils::MMFSystem::m_pi, m_RossbyDisturbance, Nektar::SolverUtils::EquationSystem::m_spacedim, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ RusanovFlux()
|
protected |
Definition at line 1097 of file MMFSWE.cpp.
1101{
1102 int nvariables = 3;
1103 NekDouble MageF1, MageF2, MageB1, MageB2;
1104 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1105 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1106
1108 NekDouble uRF, vRF, uLB, vLB;
1109 NekDouble velL, velR;
1110
1111 Array<OneD, NekDouble> EigF(nvariables);
1112 Array<OneD, NekDouble> EigB(nvariables);
1113
1114 // Compute Magnitude and Dot product of moving frames for the index
1115 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1116 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1117
1118 // Get the velocity in the normal to the edge
1121
1122 NekDouble SL, SR;
1123 SL = fabs(velL) + sqrt(g * hL);
1124 SR = fabs(velR) + sqrt(g * hR);
1125
1126 NekDouble S;
1127 if (SL > SR)
1128 {
1129 S = SL;
1130 }
1131 else
1132 {
1133 S = SR;
1134 }
1135
1136 // uRF = uR component in moving frames e^{Fwd}
1137 // vRF = vR component in moving frames e^{Fwd}
1138 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1139 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1140
1141 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1142 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1143 numfluxF[2] = 0.5 * S * (hL - hR);
1144
1145 numfluxF[3] =
1146 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
1147 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1148 numfluxF[5] = 0.5 * S * (uL * hL - uRF * hR);
1149
1150 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1151 numfluxF[7] =
1152 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1153 numfluxF[8] = 0.5 * S * (vL * hL - vRF * hR);
1154
1155 // uLB = uL component in moving frames e^{Bwd}
1156 // vLB = vL component in moving frames e^{Bwd}
1157 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1158 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1159
1160 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1161 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1162 numfluxB[2] = 0.5 * S * (hL - hR);
1163
1164 numfluxB[3] =
1165 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1166 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1167 numfluxB[5] = 0.5 * S * (uLB * hL - uR * hR);
1168
1169 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1170 numfluxB[7] =
1171 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1172 numfluxB[8] = 0.5 * S * (vLB * hL - vR * hR);
1173}
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 1389 of file MMFSWE.cpp.
1391{
1392
1394 int cnt = 0;
1395
1396 // loop over Boundary Regions
1398 {
1399
1400 // Zonal Boundary Condition
1402 {
1404 {
1406 }
1408 {
1409 // ZonalBoundary2D(n,cnt,inarray);
1410 }
1411 }
1412
1413 // Wall Boundary Condition
1415 {
1417 {
1419 }
1421 {
1422 WallBoundary2D(n, cnt, inarray);
1423 }
1424 }
1425
1426 // Time Dependent Boundary Condition (specified in meshfile)
1428 "eTimeDependent")
1429 {
1430 for (int i = 0; i < nvariables; ++i)
1431 {
1432 m_fields[i]->EvaluateBoundaryConditions(time);
1433 }
1434 }
1435 cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
1436 }
1437}
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Definition: MMFSWE.cpp:1439
References ASSERTL0, Nektar::SolverUtils::EquationSystem::m_expdim, Nektar::SolverUtils::EquationSystem::m_fields, and WallBoundary2D().
Referenced by DoOdeProjection().
◆ SteadyZonalFlow()
|
protected |
Definition at line 2059 of file MMFSWE.cpp.
2061{
2063 NekDouble uhat, vhat;
2064 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2065 NekDouble x0j, x1j, x2j, tmp;
2066
2067 Array<OneD, NekDouble> eta(nq, 0.0);
2068 Array<OneD, NekDouble> u(nq, 0.0);
2069 Array<OneD, NekDouble> v(nq, 0.0);
2070
2071 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2073 {
2074 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2075 }
2076
2077 Array<OneD, NekDouble> x(nq);
2078 Array<OneD, NekDouble> y(nq);
2079 Array<OneD, NekDouble> z(nq);
2080
2082
2083 for (int j = 0; j < nq; ++j)
2084 {
2085 x0j = x[j];
2086 x1j = y[j];
2087 x2j = z[j];
2088
2089 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2090 cos_theta);
2091
2092 // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2093 // \sin \alpha + \sin \theta \cos \alpha )^2
2094 tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
2095 sin_theta * cos(m_alpha);
2097
2098 // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2099 // || e^2 ||^2
2101 sin_theta * cos_varphi * sin(m_alpha));
2103
2104 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2105 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2106 uvec[2][j] = vhat * cos_theta;
2107 }
2108
2109 // Projection of u onto the tangent plane with conserving the mag. of the
2110 // velocity.
2111 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2112
2114 {
2115 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2116 }
2117
2118 // u is projected on the tangent plane with preserving its length
2119 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2120
2121 // Change it to the coordinate of moving frames
2122 // CartesianToMovingframes(0,uvecproj,u);
2123 // CartesianToMovingframes(1,uvecproj,v);
2124
2125 u = CartesianToMovingframes(uvec, 0);
2126 v = CartesianToMovingframes(uvec, 1);
2127
2128 switch (field)
2129 {
2130 case (0):
2131 {
2132 outfield = eta;
2133 }
2134 break;
2135
2136 case (1):
2137 {
2138 outfield = u;
2139 }
2140 break;
2141
2142 case (2):
2143 {
2144 outfield = v;
2145 }
2146 break;
2147 }
2148}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_H0, m_Hvar, Nektar::SolverUtils::EquationSystem::m_spacedim, m_u0, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ TestSWE2Dproblem()
|
private |
Definition at line 1995 of file MMFSWE.cpp.
1998{
2000
2001 Array<OneD, NekDouble> x0(nq);
2002 Array<OneD, NekDouble> x1(nq);
2003 Array<OneD, NekDouble> x2(nq);
2004
2005 m_fields[0]->GetCoords(x0, x1, x2);
2006
2007 Array<OneD, NekDouble> eta0(nq);
2008 Array<OneD, NekDouble> u0(nq);
2009 Array<OneD, NekDouble> v0(nq);
2010
2011 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2012
2014 {
2015 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2016 }
2017
2018 for (int i = 0; i < nq; ++i)
2019 {
2020 eta0[i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2021 (1.0 / cosh(0.395 * x0[i]))) *
2022 (3.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2023 exp(-0.5 * x1[i] * x1[i]);
2024 uvec[0][i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2025 (1.0 / cosh(0.395 * x0[i]))) *
2026 (-9.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2027 exp(-0.5 * x1[i] * x1[i]);
2028 uvec[1][i] = (-2.0 * 0.395 * tanh(0.395 * x0[i])) *
2029 (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2030 (1.0 / cosh(0.395 * x0[i]))) *
2031 (2.0 * x1[i]) * exp(-0.5 * x1[i] * x1[i]);
2032 }
2033
2034 u0 = CartesianToMovingframes(uvec, 0);
2035 v0 = CartesianToMovingframes(uvec, 1);
2036
2037 switch (field)
2038 {
2039 case (0):
2040 {
2041 outfield = eta0;
2042 }
2043 break;
2044
2045 case (1):
2046 {
2047 outfield = u0;
2048 }
2049 break;
2050
2051 case (2):
2052 {
2053 outfield = v0;
2054 }
2055 break;
2056 }
2057}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::EquationSystem::m_fields, and Nektar::SolverUtils::EquationSystem::m_spacedim.
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ TestVorticityComputation()
Definition at line 2953 of file MMFSWE.cpp.
2954{
2955 // Construct beta
2956 int i, k;
2957 int n = 1, m = 1;
2959
2960 NekDouble alpha, beta_theta, beta_phi;
2961
2962 NekDouble xp, yp, zp, Re;
2963 NekDouble theta, phi, sin_theta, cos_theta, sin_varphi, cos_varphi;
2964 NekDouble cosntheta3;
2965
2966 NekDouble thetax, thetay, thetaz;
2967 NekDouble phix, phiy, phiz;
2968
2969 Array<OneD, NekDouble> x(nq);
2970 Array<OneD, NekDouble> y(nq);
2971 Array<OneD, NekDouble> z(nq);
2972
2973 Array<OneD, NekDouble> u(nq);
2974 Array<OneD, NekDouble> v(nq);
2975
2976 Array<OneD, NekDouble> vorticitycompt(nq);
2977 Array<OneD, NekDouble> vorticityexact(nq);
2978
2980
2981 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2983 {
2984 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2985 }
2986
2987 // Generate SphericalCoords
2988 for (k = 0; k < nq; ++k)
2989 {
2990 xp = x[k];
2991 yp = y[k];
2992 zp = z[k];
2993
2994 Re = sqrt(xp * xp + yp * yp + zp * zp);
2995
2996 CartesianToSpherical(xp, yp, zp, sin_varphi, cos_varphi, sin_theta,
2997 cos_theta);
2998
2999 alpha = sin_varphi;
3000
3001 theta = atan2(sin_theta, cos_theta);
3002 phi = atan2(sin_varphi, cos_varphi);
3003
3004 cosntheta3 = cos(n * theta) * cos(n * theta) * cos(n * theta);
3005
3006 beta_theta = -4.0 * n * cosntheta3 * cos(m * phi) * sin(n * theta) / Re;
3007 beta_phi = -m * cosntheta3 * sin(m * phi) / Re;
3008
3009 thetax = -1.0 * cos_varphi * sin_theta;
3010 thetay = -1.0 * sin_varphi * sin_theta;
3011 thetaz = cos_theta;
3012
3013 phix = -1.0 * sin_varphi;
3014 phiy = cos_varphi;
3015 phiz = 0.0;
3016
3017 uvec[0][k] = alpha * (beta_theta * thetax + beta_phi * phix);
3018 uvec[1][k] = alpha * (beta_theta * thetay + beta_phi * phiy);
3019 uvec[2][k] = alpha * (beta_theta * thetaz + beta_phi * phiz);
3020
3021 vorticityexact[k] = -4.0 * n / Re / Re * cos_theta * cos_theta *
3022 cos_varphi * cos(m * phi) * sin(n * theta);
3023 }
3024
3025 u = CartesianToMovingframes(uvec, 0);
3026 v = CartesianToMovingframes(uvec, 1);
3027
3028 std::cout << "chi migi1" << std::endl;
3029
3030 ComputeVorticity(u, v, vorticitycompt);
3031
3032 Vmath::Vsub(nq, vorticityexact, 1, vorticitycompt, 1, vorticitycompt, 1);
3033
3036 << std::endl;
3037}
References Nektar::SolverUtils::MMFSystem::AvgAbsInt(), Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), ComputeVorticity(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_spacedim, tinysimd::sqrt(), Vmath::Vamax(), Vmath::Vsub(), and Nektar::UnitTests::z().
Referenced by v_InitObject().
◆ UnstableJetFlow()
|
protected |
Definition at line 2490 of file MMFSWE.cpp.
2493{
2495
2496 NekDouble uhat, vhat;
2497 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2498 NekDouble x0j, x1j, x2j;
2499 NekDouble Ttheta, Tphi;
2500
2501 Array<OneD, NekDouble> eta(nq, 0.0);
2502 Array<OneD, NekDouble> u(nq, 0.0);
2503 Array<OneD, NekDouble> v(nq, 0.0);
2504
2505 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2507 {
2508 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2509 }
2510
2511 Array<OneD, NekDouble> x(nq);
2512 Array<OneD, NekDouble> y(nq);
2513 Array<OneD, NekDouble> z(nq);
2514
2516
2517 int Nint = 1000;
2518 NekDouble dth, intj;
2519 for (int j = 0; j < nq; ++j)
2520 {
2521 x0j = x[j];
2522 x1j = y[j];
2523 x2j = z[j];
2524
2525 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2526 cos_theta);
2527
2528 Ttheta = atan2(sin_theta, cos_theta);
2529 Tphi = atan2(sin_varphi, cos_varphi);
2530
2531 uhat = ComputeUnstableJetuphi(Ttheta);
2532 vhat = 0.0;
2533
2534 // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2535 dth = Ttheta / Nint;
2536 eta[j] = dth * 0.5 *
2538 for (int i = 1; i < Nint - 1; i++)
2539 {
2540 intj = i * dth;
2541 eta[j] = eta[j] + dth * ComputeUnstableJetEta(intj);
2542 }
2543 eta[j] = (-1.0 / m_g) * eta[j];
2544
2545 // Add perturbation
2547 {
2548 eta[j] = eta[j] + m_hbar * cos_theta * exp(-9.0 * Tphi * Tphi) *
2549 exp(-225.0 * (m_pi / 4.0 - Ttheta) *
2550 (m_pi / 4.0 - Ttheta));
2551 }
2552
2553 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2554 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2555 uvec[2][j] = vhat * cos_theta;
2556 }
2557
2558 // Projection of u onto the tangent plane with conserving the mag. of the
2559 // velocity.
2560 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2561
2563 {
2564 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2565 }
2566
2567 // u is projected on the tangent plane with preserving its length
2568 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2569
2570 // Change it to the coordinate of moving frames
2571 // CartesianToMovingframes(0,uvecproj,u);
2572 // CartesianToMovingframes(1,uvecproj,v);
2573
2574 u = CartesianToMovingframes(uvec, 0);
2575 v = CartesianToMovingframes(uvec, 1);
2576
2577 switch (field)
2578 {
2579 case (0):
2580 {
2581 outfield = eta;
2582 }
2583 break;
2584
2585 case (1):
2586 {
2587 outfield = u;
2588 }
2589 break;
2590
2591 case (2):
2592 {
2593 outfield = v;
2594 }
2595 break;
2596 }
2597}
NekDouble ComputeUnstableJetEta(const NekDouble theta)
Definition: MMFSWE.cpp:2744
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), ComputeUnstableJetEta(), ComputeUnstableJetuphi(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_hbar, Nektar::SolverUtils::MMFSystem::m_pi, m_PurturbedJet, Nektar::SolverUtils::EquationSystem::m_spacedim, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ UnsteadyZonalFlow()
|
protected |
Definition at line 2301 of file MMFSWE.cpp.
2303{
2305 NekDouble uhat, vhat;
2306 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2307 NekDouble x0j, x1j, x2j, tmp;
2308
2309 NekDouble TR, Ttheta;
2310
2311 Array<OneD, NekDouble> eta(nq, 0.0);
2312 Array<OneD, NekDouble> u(nq, 0.0);
2313 Array<OneD, NekDouble> v(nq, 0.0);
2314
2315 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2317 {
2318 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2319 }
2320
2321 Array<OneD, NekDouble> x(nq);
2322 Array<OneD, NekDouble> y(nq);
2323 Array<OneD, NekDouble> z(nq);
2324
2326
2327 for (int j = 0; j < nq; ++j)
2328 {
2329 x0j = x[j];
2330 x1j = y[j];
2331 x2j = z[j];
2332
2333 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2334 cos_theta);
2335
2336 // \eta = ( - ( u_0 ( - T_R sin \alpha \cos \theta + \cos \alpha \sin
2337 // \theta ) + \Omega \sin \theta )^2 + \Omega \sin \theta )^2 + (\Omega
2338 // \sin \theta )^2
2339 TR =
2341 Ttheta =
2344
2348 eta[j] = 0.5 * eta[j] / m_g;
2349
2350 // u = u_0*(TR*\sin \alpha * \sin \theta + \cos \alpha * \cos \theta
2351 // v = - u_0 Ttheta * \sin \alpha
2352 uhat =
2355
2356 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2357 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2358 uvec[2][j] = vhat * cos_theta;
2359 }
2360
2361 // Projection of u onto the tangent plane with conserving the mag. of the
2362 // velocity.
2363 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2364
2366 {
2367 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2368 }
2369
2370 // u is projected on the tangent plane with preserving its length
2371 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2372
2373 // Change it to the coordinate of moving frames
2374 // CartesianToMovingframes(0,uvecproj,u);
2375 // CartesianToMovingframes(1,uvecproj,v);
2376
2377 u = CartesianToMovingframes(uvec, 0);
2378 v = CartesianToMovingframes(uvec, 1);
2379
2380 switch (field)
2381 {
2382 case (0):
2383 {
2384 outfield = eta;
2385 }
2386 break;
2387
2388 case (1):
2389 {
2390 outfield = u;
2391 }
2392 break;
2393
2394 case (2):
2395 {
2396 outfield = v;
2397 }
2398 break;
2399 }
2400}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_Omega, Nektar::SolverUtils::EquationSystem::m_spacedim, m_u0, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ v_DoInitialise()
|
overrideprotectedvirtual |
Sets up initial conditions.
Sets the initial conditions.
Reimplemented from Nektar::SolverUtils::UnsteadySystem.
Definition at line 1552 of file MMFSWE.cpp.
1553{
1554 // Compute m_depth and m_Derivdepth
1555 EvaluateWaterDepth();
1556 EvaluateCoriolis();
1557 SetInitialConditions(0.0, dumpInitialConditions);
1558 PrimitiveToConservative();
1559
1560 // transfer the initial conditions to modal values
1562 {
1565 m_fields[i]->UpdateCoeffs());
1566 }
1567}
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:674
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 248 of file MMFSWE.cpp.
249{
251
252 int i, nchk = 1;
253 int nvariables = 0;
256
258 {
259 for (i = 0; i < nfields; ++i)
260 {
261 m_intVariables.push_back(i);
262 }
263 nvariables = nfields;
264 }
265 else
266 {
267 nvariables = m_intVariables.size();
268 }
269
270 // Set up wrapper to fields data storage.
271 Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
272 Array<OneD, Array<OneD, NekDouble>> tmp(nvariables);
273
274 // Order storage to list time-integrated fields first.
275 for (i = 0; i < nvariables; ++i)
276 {
279 }
280
281 // Initialise time integration scheme
283
284 // Check uniqueness of checkpoint output
288 "Only one of IO_CheckTime and IO_CheckSteps "
289 "should be set!");
290
291 LibUtilities::Timer timer;
292 bool doCheckTime = false;
293 int step = 0;
294 NekDouble intTime = 0.0;
295 NekDouble cpuTime = 0.0;
296 NekDouble elapsed = 0.0;
297
298 NekDouble Mass = 0.0, Energy = 0.0, Enstrophy = 0.0, Vorticity = 0.0;
299 Array<OneD, NekDouble> zeta(nq);
300 Array<OneD, Array<OneD, NekDouble>> fieldsprimitive(nvariables);
301 for (int i = 0; i < nvariables; ++i)
302 {
303 fieldsprimitive[i] = Array<OneD, NekDouble>(nq);
304 }
305
307 {
308 timer.Start();
310 timer.Stop();
311
313 elapsed = timer.TimePerTest(1);
314 intTime += elapsed;
315 cpuTime += elapsed;
316
317 // Write out status information
319 m_session->GetComm()->GetRank() == 0)
320 {
321 std::cout << "Steps: " << std::setw(8) << std::left << step + 1
322 << " "
324
325 std::stringstream ss;
326 ss << cpuTime << "s";
327 std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
328 << std::endl;
329
330 // Printout Mass, Energy, Enstrophy
331 ConservativeToPrimitive(fields, fieldsprimitive);
332
333 // Vorticity zeta
334 ComputeVorticity(fieldsprimitive[1], fieldsprimitive[2], zeta);
336
337 // Masss = h^*
339
340 // Energy = 0.5*( h^*(u^2 + v^2) + g ( h^2 - h_s^s ) )
341 Energy = (ComputeEnergy(fieldsprimitive[0], fieldsprimitive[1],
342 fieldsprimitive[2]) -
343 m_Energy0) /
344 m_Energy0;
345
346 // Enstrophy = 0.5/h^* ( \mathbf{k} \cdot (\nabla \times \mathbf{v}
347 // ) + f )^2
348 Enstrophy =
349 (ComputeEnstrophy(fieldsprimitive[0], fieldsprimitive[1],
350 fieldsprimitive[2]) -
351 m_Enstrophy0) /
352 m_Enstrophy0;
353
354 std::cout << "dMass = " << std::setw(8) << std::left << Mass << " "
355 << ", dEnergy = " << std::setw(8) << std::left << Energy
356 << " "
357 << ", dEnstrophy = " << std::setw(8) << std::left
358 << Enstrophy << " "
359 << ", dVorticity = " << std::setw(8) << std::left
360 << Vorticity << std::endl
361 << std::endl;
362
363 cpuTime = 0.0;
364 }
365
366 // (h, hu, hv) -> (\eta, u, v)
367 ConservativeToPrimitive();
368
369 // Transform data into coefficient space
370 for (i = 0; i < nvariables; ++i)
371 {
376 }
377
378 // Write out checkpoint files
380 doCheckTime)
381 {
382 Checkpoint_Output(nchk++);
383 doCheckTime = false;
384
385 // (\eta, u, v) -> (h, hu, hv)
386 PrimitiveToConservative();
387 }
388
389 // Step advance
390 ++step;
391 }
392
393 // Print out summary statistics
395 {
397 {
399 << std::endl
401 }
402
404 {
405 std::cout << "Time-integration : " << intTime << "s" << std::endl;
406 }
407 }
408
409 for (i = 0; i < nvariables; ++i)
410 {
413 }
414
415 // (h, hu, hv) -> (\eta, u, v)
416 ConservativeToPrimitive();
417
418 for (i = 0; i < nvariables; ++i)
419 {
421 m_fields[i]->UpdateCoeffs());
422 }
423}
NekDouble ComputeEnstrophy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2189
NekDouble ComputeMass(const Array< OneD, const NekDouble > &eta)
Definition: MMFSWE.cpp:2150
NekDouble ComputeEnergy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2160
int m_infosteps
Number of time steps between outputting status information.
Definition: EquationSystem.h:390
SOLVER_UTILS_EXPORT void Checkpoint_Output(const int n)
Write checkpoint file of m_fields.
Definition: EquationSystem.cpp:1176
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:88
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
Definition: UnsteadySystem.h:96
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
Definition: UnsteadySystem.h:85
std::vector< int > m_intVariables
Definition: UnsteadySystem.h:93
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 3225 of file MMFSWE.cpp.
3228{
3230 {
3232 {
3233 SteadyZonalFlow(field, outfield);
3234 }
3235 break;
3236
3238 {
3239 UnsteadyZonalFlow(field, time, outfield);
3240 }
3241 break;
3242
3244 {
3245 IsolatedMountainFlow(field, time, outfield);
3246 }
3247 break;
3248
3250 {
3251 UnstableJetFlow(field, time, outfield);
3252 }
3253 break;
3254
3256 {
3257 RossbyWave(field, outfield);
3258 }
3259 break;
3260
3262 {
3263 TestSWE2Dproblem(time, field, outfield);
3264 }
3265 break;
3266
3267 default:
3268 break;
3269 }
3270}
void UnstableJetFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2490
void RossbyWave(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2599
void IsolatedMountainFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2402
void UnsteadyZonalFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2301
void TestSWE2Dproblem(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:1995
void SteadyZonalFlow(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2059
References Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, IsolatedMountainFlow(), m_TestType, RossbyWave(), SteadyZonalFlow(), TestSWE2Dproblem(), UnstableJetFlow(), and UnsteadyZonalFlow().
Referenced by v_L2Error().
◆ v_GenerateSummary()
|
overrideprotectedvirtual |
Print Summary.
Reimplemented from Nektar::SolverUtils::MMFSystem.
Definition at line 3272 of file MMFSWE.cpp.
3273{
3274 MMFSystem::v_GenerateSummary(s);
3276
3278 {
3280 {
3282 }
3283 break;
3284
3286 {
3288 m_RossbyDisturbance);
3289 }
3290 break;
3291
3293 {
3295 }
3296 break;
3297
3298 default:
3299 break;
3300 }
3301}
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 61 of file MMFSWE.cpp.
62{
63 // Call to the initialisation object
64 UnsteadySystem::v_InitObject(DeclareFields);
65
68 Array<OneD, Array<OneD, NekDouble>> Anisotropy(shapedim);
69 for (int j = 0; j < shapedim; ++j)
70 {
71 Anisotropy[j] = Array<OneD, NekDouble>(nq, 1.0);
72 }
73
74 MMFSystem::MMFInitObject(Anisotropy);
75
76 // Load acceleration of gravity
78
79 // Add Coriolois effects
81
82 // Add Rotation of the sphere along the pole
84
87
88 // Define TestType
90 "No TESTTYPE defined in session.");
93 {
95 {
97 break;
98 }
99 }
100
101 // Variable Setting for each test problem
102 NekDouble gms = 9.80616;
103
105 {
107 {
108 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
109 NekDouble rad_earth = 6.37122 * 1000000;
110 NekDouble Omegams;
111
112 // Nondimensionalized coeffs.
113 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
114
118 m_Omega = Omegams * SecondToDay;
119
122
124 }
125 break;
126
128 {
129 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
130 NekDouble rad_earth = 6.37122 * 1000000;
131 NekDouble Omegams;
132
133 // Nondimensionalized coeffs.
134 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
135
139 m_Omega = Omegams * SecondToDay;
140
143 }
144 break;
145
147 {
148 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
149 NekDouble rad_earth = 6.37122 * 1000000;
150 NekDouble Omegams;
151
152 // Nondimensionalized coeffs.
153 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
154
156
159
161 m_Omega = Omegams * SecondToDay;
162
163 m_H0 = 5960.0 / rad_earth;
164 m_hs0 = 2000.0 / rad_earth;
165 }
166 break;
167
169 {
170 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
171 NekDouble rad_earth = 6.37122 * 1000000;
172 NekDouble Omegams;
173
174 // Nondimensionalized coeffs.
175 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
176
179
181 m_Omega = Omegams * SecondToDay;
182
185
186 m_uthetamax = 80 * SecondToDay / rad_earth;
190 m_hbar = 120.0 / rad_earth;
191
195 }
196 break;
197
199 {
200 NekDouble SecondToDay = 60.0 * 60.0 * 24.0;
201 NekDouble rad_earth = 6.37122 * 1000000;
202 NekDouble Omegams;
203
204 // Nondimensionalized coeffs.
205 m_g = (gms * SecondToDay * SecondToDay) / rad_earth;
206
208 m_Omega = Omegams * SecondToDay;
209
212
213 m_angfreq = 7.848 * 0.000001 * SecondToDay;
214 m_K = 7.848 * 0.000001 * SecondToDay;
215 }
216 break;
217
218 default:
219 break;
220 }
221
222 // TestVorticityComputation
224 {
225 TestVorticityComputation();
226 }
227
228 // If explicit it computes RHS and PROJECTION for the time integration
230 {
233 }
234 // Otherwise it gives an error (no implicit integration)
235 else
236 {
238 }
239}
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:425
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:1370
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:108
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 3039 of file MMFSWE.cpp.
3042{
3044 NekDouble L2error = -1.0;
3045
3047 {
3049 {
3051 m_fields[field]->UpdatePhys());
3052 }
3053
3054 switch (field)
3055 {
3056 case (0):
3057 {
3058 // I(h) = I (h - h_exact ) / I (h_exact)
3059 Array<OneD, NekDouble> exactsolution(nq);
3061
3062 // exactsoln = u - u_T so that L2 compute u_T
3064
3066 &exactsolution[0], 1, &exactsolution[0], 1);
3067 Vmath::Vabs(nq, exactsolution, 1, exactsolution, 1);
3068
3069 L2error = (m_fields[0]->Integral(exactsolution)) / L2exact;
3070 }
3071 break;
3072
3073 case (1):
3074 {
3075 // I2 (u) = I( (u - u_ext)^2 + (v - v_ext)^2 )^{1/2} / I(
3076 // u_ext^2 + v_ext^2 )^{1/2}
3077 Array<OneD, NekDouble> exactu(nq);
3078 Array<OneD, NekDouble> exactv(nq);
3079 Array<OneD, NekDouble> tmp(nq);
3080
3081 // L2exact = \int (\sqrt{exactu*exactu+exactv*exactv})
3082 NekDouble L2exact;
3085 Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3086 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3087 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3088
3089 L2exact = m_fields[1]->Integral(tmp);
3090
3091 // L2exact = \int
3092 // (\sqrt{(u-exactu)*(u-exactu)+(v-exactv)*(v-exactv)})
3094 &exactu[0], 1);
3096 &exactv[0], 1);
3097 Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3098 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3099 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3100
3101 L2error = (m_fields[1]->Integral(tmp)) / L2exact;
3102 }
3103 break;
3104
3105 case (2):
3106 {
3107 L2error = 0.0;
3108 }
3109 break;
3110
3111 default:
3112 break;
3113 }
3114
3115 if (Normalised == true)
3116 {
3118
3121
3122 L2error = sqrt(L2error * L2error / Vol);
3123 }
3124 }
3125 else
3126 {
3127 Array<OneD, NekDouble> L2INF(2);
3128 L2INF = ErrorExtraPoints(field);
3129 L2error = L2INF[0];
3130 }
3131
3132 return L2error;
3133}
void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
Definition: MMFSWE.cpp:3225
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:899
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:774
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Definition: EquationSystem.h:442
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(), 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 3135 of file MMFSWE.cpp.
3138{
3140
3142 {
3144 m_fields[field]->UpdatePhys());
3145 }
3146
3148
3149 // Obtain \vec{u} in cartesian coordinate
3150 Array<OneD, NekDouble> x(nq);
3151 Array<OneD, NekDouble> y(nq);
3152 Array<OneD, NekDouble> z(nq);
3153
3155
3156 switch (field)
3157 {
3158 case (0):
3159 {
3160 NekDouble Letaint;
3161
3162 Array<OneD, NekDouble> exactsolution(nq);
3163
3166 exactsolution);
3167
3168 Letaint = Vmath::Vamax(nq, exactsolution, 1);
3169
3171 1, &exactsolution[0], 1);
3172
3174 }
3175 break;
3176
3177 case (1):
3178 {
3179 Array<OneD, NekDouble> exactu(nq);
3180 Array<OneD, NekDouble> exactv(nq);
3181 Array<OneD, NekDouble> tmpu(nq);
3182 Array<OneD, NekDouble> tmpv(nq);
3183 Array<OneD, NekDouble> Lerr(nq);
3184 Array<OneD, NekDouble> uT(nq);
3185
3188
3189 // Compute max[sqrt{(u-uex)^2 + (v-vex)^2}]
3192
3193 Vmath::Vsub(nq, &exactu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3194 Vmath::Vsub(nq, &exactv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3195
3196 Vmath::Vmul(nq, &tmpu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3197 Vmath::Vmul(nq, &tmpv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3198
3199 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &Lerr[0], 1);
3200 Vmath::Vsqrt(nq, &Lerr[0], 1, &Lerr[0], 1);
3201
3202 // uT = max[sqrt( u_T^2 + v_T^2 ) ]
3203 Vmath::Vmul(nq, &exactu[0], 1, &exactu[0], 1, &tmpu[0], 1);
3204 Vmath::Vmul(nq, &exactv[0], 1, &exactv[0], 1, &tmpv[0], 1);
3205 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &uT[0], 1);
3206 Vmath::Vsqrt(nq, &uT[0], 1, &uT[0], 1);
3207
3209 }
3210 break;
3211
3212 case (2):
3213 {
3214 LinfError = 0.0;
3215 }
3216 break;
3217
3218 default:
3219 break;
3220 }
3221
3223}
SOLVER_UTILS_EXPORT NekDouble LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Linf error computation.
Definition: EquationSystem.h:599
SOLVER_UTILS_EXPORT void EvaluateExactSolution(int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Evaluates an exact solution.
Definition: EquationSystem.h:682
References Nektar::SolverUtils::EquationSystem::EvaluateExactSolution(), Nektar::SolverUtils::EquationSystem::LinfError(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_time, Vmath::Vadd(), Vmath::Vamax(), Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Nektar::UnitTests::z().
◆ v_SetInitialConditions()
|
overrideprotectedvirtual |
Set the physical fields based on a restart file, or a function describing the initial condition given in the session.
- Parameters
-
initialtime Time at which to evaluate the function. dumpInitialConditions Write the initial condition to file?
Reimplemented from Nektar::SolverUtils::EquationSystem.
Definition at line 1798 of file MMFSWE.cpp.
1801{
1803
1805 {
1807 {
1808 Array<OneD, NekDouble> eta0(nq);
1809 Array<OneD, NekDouble> u0(nq);
1810 Array<OneD, NekDouble> v0(nq);
1811
1812 TestSWE2Dproblem(initialtime, 0, eta0);
1813 m_fields[0]->SetPhys(eta0);
1814
1815 TestSWE2Dproblem(initialtime, 1, u0);
1816 m_fields[1]->SetPhys(u0);
1817
1818 TestSWE2Dproblem(initialtime, 2, v0);
1819 m_fields[2]->SetPhys(v0);
1820
1821 // forward transform to fill the modal coeffs
1823 {
1826 m_fields[i]->UpdateCoeffs());
1827 }
1828 }
1829 break;
1830
1832 {
1833 Array<OneD, NekDouble> eta0(nq);
1834 Array<OneD, NekDouble> u0(nq);
1835 Array<OneD, NekDouble> v0(nq);
1836 Array<OneD, NekDouble> zeta0(nq);
1837
1838 SteadyZonalFlow(0, eta0);
1839 m_fields[0]->SetPhys(eta0);
1840
1841 SteadyZonalFlow(1, u0);
1842 m_fields[1]->SetPhys(u0);
1843
1844 SteadyZonalFlow(2, v0);
1845 m_fields[2]->SetPhys(v0);
1846
1847 // ComputeVorticity(u0, v0, zeta0);
1849
1853
1854 // forward transform to fill the modal coeffs
1856 {
1859 m_fields[i]->UpdateCoeffs());
1860 }
1861 }
1862 break;
1863
1865 {
1866 Array<OneD, NekDouble> eta0(nq);
1867 Array<OneD, NekDouble> u0(nq);
1868 Array<OneD, NekDouble> v0(nq);
1869
1870 UnsteadyZonalFlow(0, initialtime, eta0);
1871 m_fields[0]->SetPhys(eta0);
1872
1873 UnsteadyZonalFlow(1, initialtime, u0);
1874 m_fields[1]->SetPhys(u0);
1875
1876 UnsteadyZonalFlow(2, initialtime, v0);
1877 m_fields[2]->SetPhys(v0);
1878
1882
1883 // forward transform to fill the modal coeffs
1885 {
1888 m_fields[i]->UpdateCoeffs());
1889 }
1890 }
1891 break;
1892
1894 {
1895 Array<OneD, NekDouble> eta0(nq);
1896 Array<OneD, NekDouble> u0(nq);
1897 Array<OneD, NekDouble> v0(nq);
1898
1899 IsolatedMountainFlow(0, initialtime, eta0);
1900 m_fields[0]->SetPhys(eta0);
1901
1902 IsolatedMountainFlow(1, initialtime, u0);
1903 m_fields[1]->SetPhys(u0);
1904
1905 IsolatedMountainFlow(2, initialtime, v0);
1906 m_fields[2]->SetPhys(v0);
1907
1911
1912 // forward transform to fill the modal coeffs
1914 {
1917 m_fields[i]->UpdateCoeffs());
1918 }
1919 }
1920 break;
1921
1923 {
1924 Array<OneD, NekDouble> eta0(nq);
1925 Array<OneD, NekDouble> u0(nq);
1926 Array<OneD, NekDouble> v0(nq);
1927
1928 UnstableJetFlow(0, initialtime, eta0);
1929 m_fields[0]->SetPhys(eta0);
1930
1931 UnstableJetFlow(1, initialtime, u0);
1932 m_fields[1]->SetPhys(u0);
1933
1934 UnstableJetFlow(2, initialtime, v0);
1935 m_fields[2]->SetPhys(v0);
1936
1940
1941 // forward transform to fill the modal coeffs
1943 {
1946 m_fields[i]->UpdateCoeffs());
1947 }
1948 }
1949 break;
1950
1952 {
1953 Array<OneD, NekDouble> eta0(nq);
1954 Array<OneD, NekDouble> u0(nq);
1955 Array<OneD, NekDouble> v0(nq);
1956
1957 RossbyWave(0, eta0);
1958 m_fields[0]->SetPhys(eta0);
1959
1960 RossbyWave(1, u0);
1961 m_fields[1]->SetPhys(u0);
1962
1963 RossbyWave(2, v0);
1964 m_fields[2]->SetPhys(v0);
1965
1969
1970 // forward transform to fill the modal coeffs
1972 {
1975 m_fields[i]->UpdateCoeffs());
1976 }
1977 }
1978 break;
1979
1980 default:
1981 break;
1982 }
1983
1984 if (dumpInitialConditions)
1985 {
1986 // dump initial conditions to file
1988 WriteFld(outname);
1989
1991 Checkpoint_Output_Cartesian(outname);
1992 }
1993}
void Checkpoint_Output_Cartesian(std::string outname)
Definition: MMFSWE.cpp:2774
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 1439 of file MMFSWE.cpp.
1441{
1442
1443 int i;
1445 int nvariables = physarray.size();
1446
1447 // get physical values of the forward trace
1448 Array<OneD, Array<OneD, NekDouble>> Fwd0(nvariables);
1449 for (i = 0; i < nvariables; ++i)
1450 {
1451 Fwd0[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1452 m_fields[i]->ExtractTracePhys(physarray[i], Fwd0[i]);
1453 }
1454
1455 Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
1456 for (i = 0; i < nvariables; ++i)
1457 {
1458 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1459 m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1460 }
1461
1462 // Adjust the physical values of the trace to take
1463 // user defined boundaries into account
1464 int e, id1, id2, npts;
1465
1467 ++e)
1468 {
1469 npts = m_fields[0]
1470 ->GetBndCondExpansions()[bcRegion]
1471 ->GetExp(e)
1472 ->GetTotPoints();
1473 id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
1474 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
1475 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1476
1478 {
1479 case 1:
1480 {
1481 // negate the forward flux
1482 Vmath::Neg(npts, &Fwd[1][id2], 1);
1483 }
1484 break;
1485 case 2:
1486 {
1487 Array<OneD, NekDouble> tmp_n(npts);
1488 Array<OneD, NekDouble> tmp_t(npts);
1489
1491 &tmp_n[0], 1);
1493 &tmp_n[0], 1, &tmp_n[0], 1);
1494
1496 &tmp_t[0], 1);
1498 1, &tmp_t[0], 1, &tmp_t[0], 1);
1499
1500 // negate the normal flux
1501 Vmath::Neg(npts, tmp_n, 1);
1502
1503 Array<OneD, NekDouble> denom(npts);
1504 Array<OneD, NekDouble> tmp_u(npts);
1505 Array<OneD, NekDouble> tmp_v(npts);
1506
1507 // denom = (e^1 \cdot n ) (e^2 \cdot t) - (e^2 \cdot n ) (e^1
1508 // \cdot t)
1510 &m_nperpcdotMFFwd[0][id2], 1, &denom[0], 1);
1512 &m_nperpcdotMFFwd[1][id2], 1, &denom[0], 1,
1513 &denom[0], 1);
1514
1516 &tmp_u[0], 1);
1518 &tmp_u[0], 1, &tmp_u[0], 1);
1519 Vmath::Vdiv(npts, &tmp_u[0], 1, &denom[0], 1, &tmp_u[0], 1);
1520
1521 Vmath::Vcopy(npts, &tmp_u[0], 1, &Fwd[1][id2], 1);
1522
1524 &tmp_v[0], 1);
1526 &tmp_v[0], 1, &tmp_v[0], 1);
1527 Vmath::Vdiv(npts, &tmp_v[0], 1, &denom[0], 1, &tmp_v[0], 1);
1528
1529 Vmath::Vcopy(npts, &tmp_v[0], 1, &Fwd[2][id2], 1);
1530 }
1531 break;
1532
1533 default:
1535 }
1536
1537 // copy boundary adjusted values into the boundary expansion
1538 for (i = 0; i < nvariables; ++i)
1539 {
1540 Vmath::Vcopy(npts, &Fwd[i][id2], 1,
1541 &(m_fields[i]
1542 ->GetBndCondExpansions()[bcRegion]
1543 ->UpdatePhys())[id1],
1544 1);
1545 }
1546 }
1547}
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:749
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 479 of file MMFSWE.cpp.
482{
483 int i;
488
489 Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
490 Array<OneD, Array<OneD, NekDouble>> physfield(nvariables);
491
493 {
494 fluxvector[i] = Array<OneD, NekDouble>(nq);
495 }
496
497 // InField is Primitive
498 for (i = 0; i < nvariables; ++i)
499 {
500 physfield[i] = InField[i];
501 }
502
503 // Get the ith component of the flux vector in (physical space)
504 // fluxvector[0] = component for e^1 cdot \nabla \varphi
505 // fluxvector[1] = component for e^2 cdot \nabla \varphi
506 Array<OneD, NekDouble> fluxtmp(nq);
507 Array<OneD, NekDouble> tmp(ncoeffs);
508
509 // Compute Divergence Components
510 for (i = 0; i < nvariables; ++i)
511 {
512 GetSWEFluxVector(i, physfield, fluxvector);
513
514 OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
516 {
517 // Directional derivation with respect to the j'th moving frame
518 // tmp_j = ( \nabla \phi, fluxvector[j] \mathbf{e}^j )
520 fluxvector[j], tmp);
521 Vmath::Vadd(ncoeffs, &tmp[0], 1, &OutField[i][0], 1,
522 &OutField[i][0], 1);
523 }
524 }
525
526 // Numerical Flux
527 Array<OneD, Array<OneD, NekDouble>> numfluxFwd(nvariables);
528 Array<OneD, Array<OneD, NekDouble>> numfluxBwd(nvariables);
529
530 for (i = 0; i < nvariables; ++i)
531 {
532 numfluxFwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
533 numfluxBwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
534 }
535
536 NumericalSWEFlux(physfield, numfluxFwd, numfluxBwd);
537
538 // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
539 for (i = 0; i < nvariables; ++i)
540 {
541 Vmath::Neg(ncoeffs, OutField[i], 1);
542 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
543 OutField[i]);
545 }
546}
void NumericalSWEFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd)
Definition: MMFSWE.cpp:674
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:744
References Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetNpoints(), GetSWEFluxVector(), Nektar::SolverUtils::EquationSystem::GetTraceNpoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Neg(), NumericalSWEFlux(), and Vmath::Vadd().
Referenced by DoOdeRhs().
Friends And Related Function Documentation
◆ MemoryManager< MMFSWE >
|
friend |
Member Data Documentation
◆ className
|
static |
Initial value:
=
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:197
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: MMFSWE.h:67
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:99
Name of class.
◆ m_AddCoriolis
|
protected |
Definition at line 92 of file MMFSWE.h.
Referenced by DoOdeRhs(), and v_InitObject().
◆ m_AddRotation
|
protected |
Definition at line 92 of file MMFSWE.h.
Referenced by DoOdeRhs(), and v_InitObject().
◆ m_alpha
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by EvaluateCoriolisForZonalFlow(), SteadyZonalFlow(), UnsteadyZonalFlow(), v_GenerateSummary(), and v_InitObject().
◆ m_angfreq
|
protected |
Definition at line 99 of file MMFSWE.h.
Referenced by RossbyWave(), and v_InitObject().
◆ m_coriolis
Coriolis force.
Definition at line 90 of file MMFSWE.h.
Referenced by AddCoriolis(), ComputeEnstrophy(), and EvaluateCoriolis().
◆ m_depth
Still water depth.
Definition at line 86 of file MMFSWE.h.
Referenced by AddCoriolis(), AddElevationEffect(), AddRotation(), Checkpoint_Output_Cartesian(), ComputeEnergy(), ComputeEnstrophy(), ComputeMass(), ConservativeToPrimitive(), EvaluateWaterDepth(), GetSWEFluxVector(), NumericalSWEFlux(), and PrimitiveToConservative().
◆ m_Derivdepth
Definition at line 87 of file MMFSWE.h.
Referenced by AddElevationEffect(), and EvaluateWaterDepth().
◆ m_en
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_Energy0
|
protected |
Definition at line 101 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_Enstrophy0
|
protected |
Definition at line 101 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_g
|
protected |
Definition at line 94 of file MMFSWE.h.
Referenced by AddElevationEffect(), AverageFlux(), ComputeEnergy(), Computehhuhvflux(), EvaluateWaterDepth(), GetSWEFluxVector(), IsolatedMountainFlow(), LaxFriedrichFlux(), RiemannSolverHLLC(), RossbyWave(), RusanovFlux(), UnstableJetFlow(), UnsteadyZonalFlow(), and v_InitObject().
◆ m_H0
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by Checkpoint_Output_Cartesian(), ComputeEnergy(), EvaluateWaterDepth(), RossbyWave(), SteadyZonalFlow(), and v_InitObject().
◆ m_hbar
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by UnstableJetFlow(), and v_InitObject().
◆ m_hs0
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by EvaluateWaterDepth(), and v_InitObject().
◆ m_Hvar
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by SteadyZonalFlow(), and v_InitObject().
◆ m_K
|
protected |
Definition at line 99 of file MMFSWE.h.
Referenced by RossbyWave(), and v_InitObject().
◆ m_k2
|
protected |
Definition at line 96 of file MMFSWE.h.
Referenced by EvaluateWaterDepth(), and v_InitObject().
◆ m_Mass0
|
protected |
Definition at line 101 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_Omega
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by ComputeUnstableJetEta(), EvaluateCoriolisForZonalFlow(), EvaluateStandardCoriolis(), EvaluateWaterDepth(), IsolatedMountainFlow(), RossbyWave(), UnsteadyZonalFlow(), and v_InitObject().
◆ m_planeNumber
|
protected |
◆ m_primitive
|
protected |
◆ m_PurturbedJet
|
private |
Definition at line 271 of file MMFSWE.h.
Referenced by UnstableJetFlow(), v_GenerateSummary(), and v_InitObject().
◆ m_RossbyDisturbance
|
private |
Definition at line 270 of file MMFSWE.h.
Referenced by RossbyWave(), v_GenerateSummary(), and v_InitObject().
◆ m_TestType
| TestType Nektar::MMFSWE::m_TestType |
Definition at line 79 of file MMFSWE.h.
Referenced by EvaluateCoriolis(), EvaluateWaterDepth(), v_EvaluateExactSolution(), v_GenerateSummary(), v_InitObject(), and v_SetInitialConditions().
◆ m_theta0
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_theta1
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_u0
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by IsolatedMountainFlow(), SteadyZonalFlow(), UnsteadyZonalFlow(), and v_InitObject().
◆ m_uthetamax
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_veldotMF
◆ m_vellc
◆ m_velocity
Advection velocity.
Definition at line 107 of file MMFSWE.h.
Referenced by ComputeNablaCdotVelocity().
◆ m_Vorticity0
|
protected |
Definition at line 101 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
