#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_movingFrameData |
| Moving reference frame status in the inertial frame X, Y, Z, Theta_x, Theta_y, Theta_z, U, V, W, Omega_x, Omega_y, Omega_z, A_x, A_y, A_z, DOmega_x, DOmega_y, DOmega_z, pivot_x, pivot_y, pivot_z. More... | |
| std::vector< std::string > | m_strFrameData |
| variable name in m_movingFrameData More... | |
| int | m_NumQuadPointsError |
| Number of Quadrature points used to work out the error. More... | |
| enum HomogeneousType | m_HomogeneousType |
| NekDouble | m_LhomX |
| physical length in X direction (if homogeneous) More... | |
| NekDouble | m_LhomY |
| physical length in Y direction (if homogeneous) More... | |
| NekDouble | m_LhomZ |
| physical length in Z direction (if homogeneous) More... | |
| int | m_npointsX |
| number of points in X direction (if homogeneous) More... | |
| int | m_npointsY |
| number of points in Y direction (if homogeneous) More... | |
| int | m_npointsZ |
| number of points in Z direction (if homogeneous) More... | |
| int | m_HomoDirec |
| number of homogenous directions More... | |
| Protected Attributes inherited from Nektar::SolverUtils::ALEHelper | |
| Array< OneD, MultiRegions::ExpListSharedPtr > | m_fieldsALE |
| Array< OneD, Array< OneD, NekDouble > > | m_gridVelocity |
| Array< OneD, Array< OneD, NekDouble > > | m_gridVelocityTrace |
| std::vector< ALEBaseShPtr > | m_ALEs |
| bool | m_ALESolver = false |
| bool | m_ImplicitALESolver = false |
| NekDouble | m_prevStageTime = 0.0 |
| int | m_spaceDim |
Private Member Functions | |
| void | TestSWE2Dproblem (const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield) |
| void | Checkpoint_Output_Cartesian (std::string outname) |
Private Attributes | |
| int | m_RossbyDisturbance |
| int | m_PurturbedJet |
Friends | |
| class | MemoryManager< MMFSWE > |
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()
|
overridedefault |
Destructor.
◆ 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 1206 of file MMFSWE.cpp.
1208{
1209 int ncoeffs = outarray[0].size();
1210 int nq = physarray[0].size();
1211
1212 Array<OneD, NekDouble> h(nq);
1213 Array<OneD, NekDouble> tmp(nq);
1214 Array<OneD, NekDouble> tmpc(ncoeffs);
1215
1216 // physarray is primitive
1217 // conservative formulation compute h
1218 // h = \eta + d
1220
1221 int indx = 0;
1223 {
1224 if (j == 0)
1225 {
1226 indx = 2;
1227 }
1228
1229 else if (j == 1)
1230 {
1231 indx = 1;
1232 }
1233
1234 // add to hu equation
1236 Vmath::Vmul(nq, h, 1, tmp, 1, tmp, 1);
1237
1238 if (j == 1)
1239 {
1240 Vmath::Neg(nq, tmp, 1);
1241 }
1242
1243 // N \cdot (e^1 \times e^2 )
1244 // Vmath::Vmul(nq, &m_MF1crossMF2dotSN[0], 1, &tmp[0], 1, &tmp[0], 1);
1245 m_fields[0]->IProductWRTBase(tmp, tmpc);
1246 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1247 }
1248}
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 542 of file MMFSWE.cpp.
544{
545 // routine works for both primitive and conservative formulations
546 int ncoeffs = outarray[0].size();
547 int nq = physarray[0].size();
548
549 Array<OneD, NekDouble> h(nq);
550 Array<OneD, NekDouble> tmp(nq);
551 Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
553 {
554 fluxvector[i] = Array<OneD, NekDouble>(nq);
555 }
556
557 // Get 0.5 g H*H / || e^m ||^2
558 GetSWEFluxVector(3, physarray, fluxvector);
559
560 Array<OneD, NekDouble> tmpc(ncoeffs);
562 {
564
565 Vmath::Neg(nq, &tmp[0], 1);
566 m_fields[0]->IProductWRTBase(tmp, tmpc);
567
568 Vmath::Vadd(ncoeffs, outarray[j + 1], 1, tmpc, 1, outarray[j + 1], 1);
569 }
570}
void GetSWEFluxVector(const int i, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Definition: MMFSWE.cpp:572
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 1251 of file MMFSWE.cpp.
1253{
1254 int ncoeffs = outarray[0].size();
1255 int nq = physarray[0].size();
1256
1257 Array<OneD, NekDouble> h(nq);
1258 Array<OneD, NekDouble> tmp(nq);
1259 Array<OneD, NekDouble> tmpc(ncoeffs);
1260
1261 // physarray is primitive
1262 // conservative formulation compute h
1263 // h = \eta + d
1265
1267 {
1270
1271 m_fields[0]->IProductWRTBase(tmp, tmpc);
1272
1273 Vmath::Vadd(ncoeffs, tmpc, 1, outarray[j + 1], 1, outarray[j + 1], 1);
1274 }
1275}
Array< OneD, Array< OneD, NekDouble > > m_Derivdepth
Definition: MMFSWE.h:86
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 1280 of file MMFSWE.cpp.
1282{
1283 // routine works for both primitive and conservative formulations
1284 int ncoeffs = outarray[0].size();
1285 int nq = physarray[0].size();
1286
1287 // Compute h
1288 Array<OneD, NekDouble> h(nq);
1290
1291 Array<OneD, NekDouble> de0dt_cdot_e0;
1292 Array<OneD, NekDouble> de0dt_cdot_e1;
1293 Array<OneD, NekDouble> de1dt_cdot_e0;
1294 Array<OneD, NekDouble> de1dt_cdot_e1;
1295 Compute_demdt_cdot_ek(0, 0, physarray, de0dt_cdot_e0);
1296 Compute_demdt_cdot_ek(1, 0, physarray, de1dt_cdot_e0);
1297 Compute_demdt_cdot_ek(0, 1, physarray, de0dt_cdot_e1);
1298 Compute_demdt_cdot_ek(1, 1, physarray, de1dt_cdot_e1);
1299
1300 Array<OneD, NekDouble> Rott1(nq);
1301 Array<OneD, NekDouble> Rott2(nq);
1302 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e0, 1, Rott1, 1);
1303 Vmath::Vmul(nq, physarray[1], 1, de0dt_cdot_e1, 1, Rott2, 1);
1304 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e0, 1, Rott1, 1, Rott1, 1);
1305 Vmath::Vvtvp(nq, physarray[2], 1, de1dt_cdot_e1, 1, Rott2, 1, Rott2, 1);
1306
1307 // Multiply H and \partial \phi / \partial t which is assumed to be u_{\phi}
1308 Vmath::Vmul(nq, &h[0], 1, &Rott1[0], 1, &Rott1[0], 1);
1309 Vmath::Vmul(nq, &h[0], 1, &Rott2[0], 1, &Rott2[0], 1);
1310
1311 Vmath::Neg(nq, Rott1, 1);
1312 Vmath::Neg(nq, Rott2, 1);
1313
1314 Array<OneD, NekDouble> tmpc1(ncoeffs);
1315 Array<OneD, NekDouble> tmpc2(ncoeffs);
1316 m_fields[0]->IProductWRTBase(Rott1, tmpc1);
1317 m_fields[0]->IProductWRTBase(Rott2, tmpc2);
1318
1319 Vmath::Vadd(ncoeffs, tmpc1, 1, outarray[1], 1, outarray[1], 1);
1320 Vmath::Vadd(ncoeffs, tmpc2, 1, outarray[2], 1, outarray[2], 1);
1321}
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:1327
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 961 of file MMFSWE.cpp.
965{
966 NekDouble MageF1, MageF2, MageB1, MageB2;
967 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
968 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
969
971 NekDouble uRF, vRF, uLB, vLB;
972
973 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
974 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
975
976 // uRF = uR component in moving frames e^{Fwd}
977 // vRF = vR component in moving frames e^{Fwd}
978 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
979 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
980
981 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
982 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
983 numfluxF[2] = 0.0;
984
985 numfluxF[3] =
986 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
987 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
988 numfluxF[5] = 0.0;
989
990 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
991 numfluxF[7] =
992 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
993 numfluxF[8] = 0.0;
994
995 // uLB = uL component in moving frames e^{Bwd}
996 // vLB = vL component in moving frames e^{Bwd}
997 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
998 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
999
1000 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1001 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1002 numfluxB[2] = 0.0;
1003
1004 numfluxB[3] =
1005 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1006 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1007 numfluxB[5] = 0.0;
1008
1009 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1010 numfluxB[7] =
1011 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1012 numfluxB[8] = 0.0;
1013}
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:1168
References ComputeMagAndDot(), and m_g.
Referenced by NumericalSWEFlux().
◆ Checkpoint_Output_Cartesian()
|
private |
Definition at line 2767 of file MMFSWE.cpp.
2768{
2771
2772 NekDouble rad_earth = 6.37122 * 1000000;
2773
2774 int nvariables = 7;
2775
2776 // vector u in Cartesian coordinates
2777 std::vector<std::string> variables(nvariables);
2778
2779 variables[0] = "eta";
2780 variables[1] = "hstar";
2781 variables[2] = "vorticity";
2782 variables[3] = "ux";
2783 variables[4] = "uy";
2784 variables[5] = "uz";
2785 variables[6] = "null";
2786
2787 // Obtain \vec{u} in cartesian coordinate
2788 Array<OneD, Array<OneD, NekDouble>> fieldphys(nvariables);
2789 std::vector<Array<OneD, NekDouble>> fieldcoeffs(nvariables);
2790 for (int i = 0; i < nvariables; ++i)
2791 {
2792 fieldphys[i] = Array<OneD, NekDouble>(nq, 0.0);
2793 fieldcoeffs[i] = Array<OneD, NekDouble>(ncoeffs);
2794 }
2795
2797 &fieldphys[0][0], 1);
2798 // Vmath::Vadd(nq, &(m_fields[0]->GetPhys())[0], 1, &m_depth[0], 1,
2799 // &fieldphys[1][0], 1);
2800
2802 Vmath::Neg(nq, &fieldphys[1][0], 1);
2804 Vmath::Smul(nq, rad_earth, &fieldphys[1][0], 1, &fieldphys[1][0], 1);
2805
2806 Array<OneD, NekDouble> utmp(nq);
2807 Array<OneD, NekDouble> vtmp(nq);
2808
2811
2812 ComputeVorticity(utmp, vtmp, fieldphys[2]);
2813 // u_x = u e^1_x + v e^2_x
2814 int indx = 3;
2816 {
2818 &fieldphys[k + indx][0], 1);
2820 &fieldphys[k + indx][0], 1, &fieldphys[k + indx][0], 1);
2821 }
2822
2823 for (int i = 0; i < nvariables; ++i)
2824 {
2825 m_fields[0]->FwdTrans(fieldphys[i], fieldcoeffs[i]);
2826 }
2827
2829}
void ComputeVorticity(const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v, Array< OneD, NekDouble > &Vorticity)
Definition: MMFSWE.cpp:2213
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1257
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 1327 of file MMFSWE.cpp.
1331{
1332 int j, k;
1334
1335 Array<OneD, NekDouble> tmp(nq, 0.0);
1336 Array<OneD, NekDouble> tmpderiv(nq);
1337
1338 outarray = Array<OneD, NekDouble>(nq, 0.0);
1340 {
1342 {
1343 // Compute d e^m / d \xi_1 and d e^m / d \xi_2
1346
1347 Vmath::Vmul(nq, &physarray[j + 1][0], 1, &tmpderiv[0], 1,
1348 &tmpderiv[0], 1);
1349
1351 &outarray[0], 1, &outarray[0], 1);
1352 }
1353 }
1354}
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 2153 of file MMFSWE.cpp.
2156{
2158
2159 Array<OneD, NekDouble> tmp(nq);
2160 Array<OneD, NekDouble> htmp(nq);
2161 Array<OneD, NekDouble> hstmp(nq);
2162
2163 Vmath::Vmul(nq, u, 1, u, 1, tmp, 1);
2164 Vmath::Vvtvp(nq, v, 1, v, 1, tmp, 1, tmp, 1);
2165 Vmath::Vmul(nq, eta, 1, tmp, 1, tmp, 1);
2166
2168 Vmath::Vmul(nq, htmp, 1, htmp, 1, htmp, 1);
2169
2171 Vmath::Vmul(nq, hstmp, 1, hstmp, 1, hstmp, 1);
2172
2173 Vmath::Vsub(nq, htmp, 1, hstmp, 1, htmp, 1);
2175
2176 Vmath::Vadd(nq, htmp, 1, tmp, 1, tmp, 1);
2177 Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2178
2180}
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 2182 of file MMFSWE.cpp.
2185{
2187
2188 Array<OneD, NekDouble> hstartmp(nq);
2189 Array<OneD, NekDouble> tmp(nq);
2190
2192
2193 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2194 Array<OneD, Array<OneD, NekDouble>> Curlu(m_spacedim);
2196 {
2197 uvec[i] = Array<OneD, NekDouble>(nq);
2198 Curlu[i] = Array<OneD, NekDouble>(nq);
2199 }
2200
2201 ComputeVorticity(u, v, tmp);
2202
2204 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
2205 Vmath::Vdiv(nq, tmp, 1, hstartmp, 1, tmp, 1);
2206 Vmath::Smul(nq, 0.5, tmp, 1, tmp, 1);
2207
2209}
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 883 of file MMFSWE.cpp.
887{
889
890 NekDouble hC, huC, hvC, SL, SR, Sstar;
893
894 // Compute SL
895 if (hstar > hL)
896 {
897 SL = uL - cL * sqrt(0.5 * ((hstar * hstar + hstar * hL) / (hL * hL)));
898 }
899 else
900 {
901 SL = uL - cL;
902 }
903
904 // Compute SR
905 if (hstar > hR)
906 {
907 SR = uR + cR * sqrt(0.5 * ((hstar * hstar + hstar * hR) / (hR * hR)));
908 }
909 else
910 {
911 SR = uR + cR;
912 }
913
914 if (fabs(hR * (uR - SR) - hL * (uL - SL)) <= 1.0e-15)
915 {
916 Sstar = 0.0;
917 }
918 else
919 {
920 Sstar = (SL * hR * (uR - SR) - SR * hL * (uL - SL)) /
921 (hR * (uR - SR) - hL * (uL - SL));
922 }
923
924 if (SL >= 0)
925 {
926 hflux = hL * uL;
927 huflux = uL * uL * hL + 0.5 * g * hL * hL;
928 hvflux = hL * uL * vL;
929 }
930 else if (SR <= 0)
931 {
932 hflux = hR * uR;
933 huflux = uR * uR * hR + 0.5 * g * hR * hR;
934 hvflux = hR * uR * vR;
935 }
936 else
937 {
938 if ((SL < 0) && (Sstar >= 0))
939 {
940 hC = hL * ((SL - uL) / (SL - Sstar));
941 huC = hC * Sstar;
942 hvC = hC * vL;
943
944 hflux = hL * uL + SL * (hC - hL);
945 huflux = (uL * uL * hL + 0.5 * g * hL * hL) + SL * (huC - hL * uL);
946 hvflux = (uL * vL * hL) + SL * (hvC - hL * vL);
947 }
948 else
949 {
950 hC = hR * ((SR - uR) / (SR - Sstar));
951 huC = hC * Sstar;
952 hvC = hC * vR;
953
954 hflux = hR * uR + SR * (hC - hR);
955 huflux = (uR * uR * hR + 0.5 * g * hR * hR) + SR * (huC - hR * uR);
956 hvflux = (uR * vR * hR) + SR * (hvC - hR * vR);
957 }
958 }
959}
References m_g, and tinysimd::sqrt().
Referenced by RiemannSolverHLLC().
◆ ComputeMagAndDot()
|
protected |
Definition at line 1168 of file MMFSWE.cpp.
1173{
1174 NekDouble MF1x, MF1y, MF1z, MF2x, MF2y, MF2z;
1175 NekDouble MB1x, MB1y, MB1z, MB2x, MB2y, MB2z;
1176
1177 MF1x = m_MFtraceFwd[0][0][index];
1178 MF1y = m_MFtraceFwd[0][1][index];
1179 MF1z = m_MFtraceFwd[0][2][index];
1180
1181 MF2x = m_MFtraceFwd[1][0][index];
1182 MF2y = m_MFtraceFwd[1][1][index];
1183 MF2z = m_MFtraceFwd[1][2][index];
1184
1185 MB1x = m_MFtraceBwd[0][0][index];
1186 MB1y = m_MFtraceBwd[0][1][index];
1187 MB1z = m_MFtraceBwd[0][2][index];
1188
1189 MB2x = m_MFtraceBwd[1][0][index];
1190 MB2y = m_MFtraceBwd[1][1][index];
1191 MB2z = m_MFtraceBwd[1][2][index];
1192
1193 // MFtrace = MFtrace [ j*spacedim + k ], j = shape, k = sapce
1194 MageF1 = MF1x * MF1x + MF1y * MF1y + MF1z * MF1z;
1195 MageF2 = MF2x * MF2x + MF2y * MF2y + MF2z * MF2z;
1196 MageB1 = MB1x * MB1x + MB1y * MB1y + MB1z * MB1z;
1197 MageB2 = MB2x * MB2x + MB2y * MB2y + MB2z * MB2z;
1198
1199 eF1_cdot_eB1 = MF1x * MB1x + MF1y * MB1y + MF1z * MB1z;
1200 eF1_cdot_eB2 = MF1x * MB2x + MF1y * MB2y + MF1z * MB2z;
1201 eF2_cdot_eB1 = MF2x * MB1x + MF2y * MB1y + MF2z * MB1z;
1202 eF2_cdot_eB2 = MF2x * MB2x + MF2y * MB2y + MF2z * MB2z;
1203}
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 2143 of file MMFSWE.cpp.
2144{
2146
2147 Array<OneD, NekDouble> tmp(nq);
2149
2151}
References m_depth, Nektar::SolverUtils::EquationSystem::m_fields, and Vmath::Vadd().
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ ComputeNablaCdotVelocity()
Definition at line 2234 of file MMFSWE.cpp.
2235{
2237
2238 Array<OneD, NekDouble> velcoeff(nq, 0.0);
2239
2240 Array<OneD, NekDouble> Dtmp0(nq);
2241 Array<OneD, NekDouble> Dtmp1(nq);
2242 Array<OneD, NekDouble> Dtmp2(nq);
2243 Array<OneD, NekDouble> Drv(nq);
2244
2245 vellc = Array<OneD, NekDouble>(nq, 0.0);
2246
2247 // m_vellc = \nabla m_vel \cdot tan_i
2248 Array<OneD, NekDouble> tmp(nq);
2249 Array<OneD, NekDouble> vessel(nq);
2250
2252 {
2253 Vmath::Zero(nq, velcoeff, 1);
2255 {
2256 // a_j = tan_j cdot m_vel
2258 1, &velcoeff[0], 1, &velcoeff[0], 1);
2259 }
2260
2261 // d a_j / d x^k
2262 m_fields[0]->PhysDeriv(velcoeff, Dtmp0, Dtmp1, Dtmp2);
2263
2265 {
2266 // tan_j^k ( d a_j / d x^k )
2267 switch (k)
2268 {
2269 case (0):
2270 {
2272 1, &vellc[0], 1, &vellc[0], 1);
2273 }
2274 break;
2275
2276 case (1):
2277 {
2279 1, &vellc[0], 1, &vellc[0], 1);
2280 }
2281 break;
2282
2283 case (2):
2284 {
2286 1, &vellc[0], 1, &vellc[0], 1);
2287 }
2288 break;
2289 }
2290 }
2291 }
2292}
References Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_shapedim, Nektar::SolverUtils::EquationSystem::m_spacedim, m_velocity, Vmath::Vvtvp(), and Vmath::Zero().
◆ ComputeUnstableJetEta()
Definition at line 2737 of file MMFSWE.cpp.
2738{
2739 NekDouble uphi, f, dh;
2740
2741 uphi = ComputeUnstableJetuphi(theta);
2742 f = 2.0 * m_Omega * sin(theta);
2743
2744 dh = f * uphi + tan(theta) * uphi * uphi;
2745
2746 return dh;
2747}
NekDouble ComputeUnstableJetuphi(const NekDouble theta)
Definition: MMFSWE.cpp:2749
References ComputeUnstableJetuphi(), and m_Omega.
Referenced by UnstableJetFlow().
◆ ComputeUnstableJetuphi()
Definition at line 2749 of file MMFSWE.cpp.
2750{
2751 NekDouble uphi;
2752
2754 {
2757 }
2758
2759 else
2760 {
2761 uphi = 0.0;
2762 }
2763
2764 return uphi;
2765}
References m_en, m_theta0, m_theta1, and m_uthetamax.
Referenced by ComputeUnstableJetEta(), and UnstableJetFlow().
◆ ComputeVorticity()
|
protected |
Definition at line 2213 of file MMFSWE.cpp.
2216{
2218
2219 Array<OneD, NekDouble> tmp(nq);
2220
2221 Vorticity = Array<OneD, NekDouble>(nq, 0.0);
2222
2225 &Vorticity[0], 1);
2226
2228 Vmath::Neg(nq, tmp, 1);
2230
2231 Vmath::Vadd(nq, tmp, 1, Vorticity, 1, Vorticity, 1);
2232}
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 2912 of file MMFSWE.cpp.
2913{
2915
2916 // u = hu/h
2918 m_fields[1]->UpdatePhys(), 1);
2919
2920 // v = hv/ v
2922 m_fields[2]->UpdatePhys(), 1);
2923
2924 // \eta = h - d
2926 m_fields[0]->UpdatePhys(), 1);
2927}
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:759
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Vdiv(), and Vmath::Vsub().
Referenced by DoOdeProjection(), DoOdeRhs(), and v_DoSolve().
◆ ConservativeToPrimitive() [2/2]
|
protected |
Definition at line 2832 of file MMFSWE.cpp.
2835{
2837
2838 if (physin[0].get() == physout[0].get())
2839 {
2840 // copy indata and work with tmp array
2841 Array<OneD, Array<OneD, NekDouble>> tmp(3);
2842 for (int i = 0; i < 3; ++i)
2843 {
2844 // deep copy
2845 tmp[i] = Array<OneD, NekDouble>(nq);
2846 Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2847 }
2848
2849 // \eta = h - d
2851
2852 // u = hu/h
2853 Vmath::Vdiv(nq, tmp[1], 1, tmp[0], 1, physout[1], 1);
2854
2855 // v = hv/h
2856 Vmath::Vdiv(nq, tmp[2], 1, tmp[0], 1, physout[2], 1);
2857 }
2858 else
2859 {
2860 // \eta = h - d
2862
2863 // u = hu/h
2864 Vmath::Vdiv(nq, physin[1], 1, physin[0], 1, physout[1], 1);
2865
2866 // v = hv/h
2867 Vmath::Vdiv(nq, physin[2], 1, physin[0], 1, physout[2], 1);
2868 }
2869}
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Vmath::Vcopy(), Vmath::Vdiv(), and Vmath::Vsub().
◆ create()
|
inlinestatic |
Creates an instance of this class.
Definition at line 66 of file MMFSWE.h.
69 {
71 MemoryManager<MMFSWE>::AllocateSharedPtr(pSession, pGraph);
72 p->InitObject();
74 }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:166
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition: EquationSystem.h:65
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.
◆ DoOdeProjection()
|
protected |
Compute the projection.
Compute the projection for the linear advection equation.
- Parameters
-
inarray Given fields. outarray Calculated solution. time Time.
Definition at line 1363 of file MMFSWE.cpp.
1366{
1368 {
1370 {
1371 ConservativeToPrimitive(inarray, outarray);
1372 SetBoundaryConditions(outarray, time);
1373 PrimitiveToConservative(outarray, outarray);
1374 }
1375 break;
1376 default:
1378 break;
1379 }
1380}
void SetBoundaryConditions(Array< OneD, Array< OneD, NekDouble > > &inarray, NekDouble time)
Definition: MMFSWE.cpp:1382
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 418 of file MMFSWE.cpp.
421{
422 int i;
423 int nvariables = inarray.size();
426
427 // inarray in physical space
428 Array<OneD, Array<OneD, NekDouble>> physarray(nvariables);
429 Array<OneD, Array<OneD, NekDouble>> modarray(nvariables);
430
431 for (i = 0; i < nvariables; ++i)
432 {
433 physarray[i] = Array<OneD, NekDouble>(nq);
434 modarray[i] = Array<OneD, NekDouble>(ncoeffs);
435 }
436
437 // (h, hu, hv) -> (\eta, u, v)
438 ConservativeToPrimitive(inarray, physarray);
439
440 // Weak Directional Derivative
441 WeakDGSWEDirDeriv(physarray, modarray);
442 AddDivForGradient(physarray, modarray);
443
444 for (i = 0; i < nvariables; ++i)
445 {
446 Vmath::Neg(ncoeffs, modarray[i], 1);
447 }
448
449 // coriolis forcing
451 {
452 AddCoriolis(physarray, modarray);
453 }
454
455 // Bottom Elevation Effect
456 // Add g H \nabla m_depth
457 AddElevationEffect(physarray, modarray);
458
459 // Add terms concerning the rotation of the moving frame
461 {
462 AddRotation(physarray, modarray);
463 }
464
465 for (i = 0; i < nvariables; ++i)
466 {
467 m_fields[i]->MultiplyByElmtInvMass(modarray[i], modarray[i]);
468 m_fields[i]->BwdTrans(modarray[i], outarray[i]);
469 }
470}
void AddDivForGradient(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:542
void AddCoriolis(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1206
void WeakDGSWEDirDeriv(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField)
Definition: MMFSWE.cpp:472
void AddRotation(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1280
void AddElevationEffect(Array< OneD, Array< OneD, NekDouble > > &physarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: MMFSWE.cpp:1251
SOLVER_UTILS_EXPORT int GetNcoeffs()
Definition: EquationSystem.h:701
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 1702 of file MMFSWE.cpp.
1703{
1705 {
1707 {
1709 }
1710 break;
1711
1713 {
1715 }
1716 break;
1717
1722 {
1724 }
1725 break;
1726
1727 default:
1728 break;
1729 }
1730}
void EvaluateStandardCoriolis(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1764
void EvaluateCoriolisForZonalFlow(Array< OneD, NekDouble > &outarray)
Definition: MMFSWE.cpp:1732
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 1732 of file MMFSWE.cpp.
1733{
1735 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1736
1737 Array<OneD, NekDouble> x(nq);
1738 Array<OneD, NekDouble> y(nq);
1739 Array<OneD, NekDouble> z(nq);
1740
1742
1743 NekDouble x0j, x1j, x2j;
1744 NekDouble tmp;
1745
1746 outarray = Array<OneD, NekDouble>(nq, 0.0);
1747 for (int j = 0; j < nq; ++j)
1748 {
1749 x0j = x[j];
1750 x1j = y[j];
1751 x2j = z[j];
1752
1753 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1754 cos_theta);
1755
1756 // H = 2 \Omega *(- \cos \phi \cos \theta \sin \alpha + \sin \theta \cos
1757 // \alpha )
1758 tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
1759 sin_theta * cos(m_alpha);
1760 outarray[j] = 2.0 * m_Omega * tmp;
1761 }
1762}
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 1764 of file MMFSWE.cpp.
1765{
1767
1768 NekDouble x0j, x1j, x2j;
1769 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
1770
1771 Array<OneD, NekDouble> x(nq);
1772 Array<OneD, NekDouble> y(nq);
1773 Array<OneD, NekDouble> z(nq);
1774
1776
1777 outarray = Array<OneD, NekDouble>(nq, 0.0);
1778 for (int j = 0; j < nq; ++j)
1779 {
1780 x0j = x[j];
1781 x1j = y[j];
1782 x2j = z[j];
1783
1784 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
1785 cos_theta);
1786
1787 outarray[j] = 2.0 * m_Omega * sin_theta;
1788 }
1789}
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 1562 of file MMFSWE.cpp.
1563{
1565
1566 m_depth = Array<OneD, NekDouble>(nq, 0.0);
1567
1569 {
1571 {
1573 }
1574 break;
1575
1577 {
1578 // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1579 Array<OneD, NekDouble> x(nq);
1580 Array<OneD, NekDouble> y(nq);
1581 Array<OneD, NekDouble> z(nq);
1582
1584
1585 NekDouble x0j, x1j, x2j;
1586 NekDouble sin_varphi, cos_varphi, sin_theta, cos_theta;
1587
1588 for (int j = 0; j < nq; ++j)
1589 {
1590 x0j = x[j];
1591 x1j = y[j];
1592 x2j = z[j];
1593
1594 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1595 sin_theta, cos_theta);
1596 m_depth[j] =
1599 }
1600 }
1601 break;
1602
1604 {
1605 // H_0 = k_1 - k_2 - (1/2/g) (Omega sin \phi)
1606 Array<OneD, NekDouble> x(nq);
1607 Array<OneD, NekDouble> y(nq);
1608 Array<OneD, NekDouble> z(nq);
1609
1611
1612 int indx = 0;
1613 NekDouble x0j, x1j, x2j, dist2;
1614 NekDouble phi, theta, sin_varphi, cos_varphi, sin_theta, cos_theta;
1615 NekDouble hRad, phic, thetac;
1616 NekDouble Tol = 0.000001;
1617
1618 hRad = m_pi / 9.0;
1619 phic = -m_pi / 2.0;
1620 thetac = m_pi / 6.0;
1621
1622 for (int j = 0; j < nq; ++j)
1623 {
1624 x0j = x[j];
1625 x1j = y[j];
1626 x2j = z[j];
1627
1628 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi,
1629 sin_theta, cos_theta);
1630
1632 std::abs(sin(thetac) - sin_theta)) < Tol)
1633 {
1634 std::cout << "A point " << j
1635 << " is coincient with the singularity "
1636 << std::endl;
1637 indx = 1;
1638 }
1639
1640 phi = atan2(sin_varphi, cos_varphi);
1641 theta = atan2(sin_theta, cos_theta);
1642
1643 // Compute r
1644 dist2 = (phi - phic) * (phi - phic) +
1645 (theta - thetac) * (theta - thetac);
1646
1647 if (dist2 > hRad * hRad)
1648 {
1649 dist2 = hRad * hRad;
1650 }
1651
1653 }
1654
1655 if (!indx)
1656 {
1657 std::cout << "No point is coincident with the singularity point"
1658 << std::endl;
1659 }
1660 }
1661 break;
1662
1664 {
1665 for (int j = 0; j < nq; ++j)
1666 {
1668 }
1669 }
1670 break;
1671
1674 {
1676 }
1677 break;
1678
1679 default:
1680 {
1682 }
1683 break;
1684 }
1685
1686 // Comptue \nabla m_depth \cdot e^i
1688
1690 {
1691 m_Derivdepth[j] = Array<OneD, NekDouble>(nq);
1693 m_Derivdepth[j]);
1694 }
1695
1696 std::cout << "Water Depth (m_depth) was generated with mag = "
1700}
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 572 of file MMFSWE.cpp.
575{
577
578 switch (i)
579 {
580 // flux function for the h equation = [(\eta + d ) u, (\eta + d) v ]
581 case 0:
582 {
583 // h in flux 1
585
586 // hu in flux 0
587 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
588
589 // hv in flux 1
590 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
591 }
592 break;
593
594 // flux function for the hu equation = [Hu^2 + 0.5 g H^2, Huv]
595 case 1:
596 {
597 Array<OneD, NekDouble> tmp(nq);
598
599 // h in tmp
601
602 // hu in flux 1
603 Vmath::Vmul(nq, tmp, 1, physfield[1], 1, flux[1], 1);
604
605 // huu in flux 0
606 Vmath::Vmul(nq, flux[1], 1, physfield[1], 1, flux[0], 1);
607
608 // hh in tmp
609 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
610
611 // huu + 0.5 g hh in flux 0
612 // Daxpy overwrites flux[0] on exit
614
615 // huv in flux 1
616 Vmath::Vmul(nq, flux[1], 1, physfield[2], 1, flux[1], 1);
617 }
618 break;
619
620 // flux function for the hv equation = [Huv, Hv^2]
621 case 2:
622 {
623 Array<OneD, NekDouble> tmp(nq);
624
625 // h in tmp
627
628 // hv in flux 0
629 Vmath::Vmul(nq, tmp, 1, physfield[2], 1, flux[0], 1);
630
631 // hvv in flux 1
632 Vmath::Vmul(nq, flux[0], 1, physfield[2], 1, flux[1], 1);
633
634 // huv in flux 0
635 Vmath::Vmul(nq, flux[0], 1, physfield[1], 1, flux[0], 1);
636
637 // hh in tmp
638 Vmath::Vmul(nq, tmp, 1, tmp, 1, tmp, 1);
639
640 // hvv + 0.5 g hh in flux 1
642 }
643 break;
644
645 // flux function 0.5 g h * h
646 case 3:
647 {
648 Array<OneD, NekDouble> h(nq);
649 flux[0] = Array<OneD, NekDouble>(nq, 0.0);
650
651 // h in tmp
653
654 // hh in tmp
655 Vmath::Vmul(nq, h, 1, h, 1, h, 1);
656
657 // 0.5 g hh in flux 0
659 }
660 break;
661
662 default:
663 break;
664 }
665}
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 2395 of file MMFSWE.cpp.
2398{
2400
2401 NekDouble uhat, vhat;
2402 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2403 NekDouble x0j, x1j, x2j;
2404
2405 Array<OneD, NekDouble> eta(nq, 0.0);
2406 Array<OneD, NekDouble> u(nq, 0.0);
2407 Array<OneD, NekDouble> v(nq, 0.0);
2408
2409 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2411 {
2412 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2413 }
2414
2415 Array<OneD, NekDouble> x(nq);
2416 Array<OneD, NekDouble> y(nq);
2417 Array<OneD, NekDouble> z(nq);
2418
2420
2421 for (int j = 0; j < nq; ++j)
2422 {
2423 x0j = x[j];
2424 x1j = y[j];
2425 x2j = z[j];
2426
2427 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2428 cos_theta);
2429
2430 // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2432 sin_theta * sin_theta;
2433
2434 uhat = m_u0 * cos_theta;
2435 vhat = 0.0;
2436
2437 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2438 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2439 uvec[2][j] = vhat * cos_theta;
2440 }
2441
2442 // Projection of u onto the tangent plane with conserving the mag. of the
2443 // velocity.
2444 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2445
2447 {
2448 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2449 }
2450
2451 // u is projected on the tangent plane with preserving its length
2452 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2453
2454 // Change it to the coordinate of moving frames
2455 // CartesianToMovingframes(0,uvecproj,u);
2456 // CartesianToMovingframes(1,uvecproj,v);
2457
2458 u = CartesianToMovingframes(uvec, 0);
2459 v = CartesianToMovingframes(uvec, 1);
2460
2462 {
2463 case (0):
2464 {
2465 outfield = eta;
2466 }
2467 break;
2468
2469 case (1):
2470 {
2471 outfield = u;
2472 }
2473 break;
2474
2475 case (2):
2476 {
2477 outfield = v;
2478 }
2479 break;
2480 }
2481}
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble > > &uvec, unsigned int field)
Definition: MMFSystem.cpp:771
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_Omega, Nektar::SolverUtils::EquationSystem::m_spacedim, m_u0, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ LaxFriedrichFlux()
|
protected |
Definition at line 1015 of file MMFSWE.cpp.
1019{
1020 int nvariables = 3;
1021 NekDouble MageF1, MageF2, MageB1, MageB2;
1022 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1023 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1024
1026 NekDouble uRF, vRF, uLB, vLB;
1027 NekDouble velL, velR, lambdaF, lambdaB;
1028
1029 Array<OneD, NekDouble> EigF(nvariables);
1030 Array<OneD, NekDouble> EigB(nvariables);
1031
1032 // Compute Magnitude and Dot product of moving frames for the index
1033 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1034 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1035
1036 // Get the velocity in the normal to the edge
1039
1040 EigF[0] = velL - sqrt(g * hL);
1041 EigF[1] = velL;
1042 EigF[2] = velL + sqrt(g * hL);
1043
1044 EigB[0] = velR - sqrt(g * hR);
1045 EigB[1] = velR;
1046 EigB[2] = velR + sqrt(g * hR);
1047
1048 lambdaF = Vmath::Vamax(nvariables, EigF, 1);
1049 lambdaB = Vmath::Vamax(nvariables, EigB, 1);
1050
1051 // uRF = uR component in moving frames e^{Fwd}
1052 // vRF = vR component in moving frames e^{Fwd}
1053 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1054 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1055
1056 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1057 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1058 numfluxF[2] = 0.5 * lambdaF * (hL - hR);
1059
1060 numfluxF[3] = 0.5 * (hL * uL * uL * MageF1 + hR * uRF * uRF * MageB1 +
1061 0.5 * g * (hL * hL + hR * hR));
1062 numfluxF[4] = 0.5 * (hL * uL * vL * MageF1 + hR * uRF * vRF * MageB1);
1063 numfluxF[5] = 0.5 * lambdaF * (uL * hL - uRF * hR);
1064
1065 numfluxF[6] = 0.5 * (hL * uL * vL * MageF2 + hR * uRF * vRF * MageB2);
1066 numfluxF[7] = 0.5 * (hL * vL * vL * MageF2 + hR * vRF * vRF * MageB2 +
1067 0.5 * g * (hL * hL + hR * hR));
1068 numfluxF[8] = 0.5 * lambdaF * (vL * hL - vRF * hR);
1069
1070 // uLB = uL component in moving frames e^{Bwd}
1071 // vLB = vL component in moving frames e^{Bwd}
1072 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1073 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1074
1075 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1076 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1077 numfluxB[2] = 0.5 * lambdaB * (hL - hR);
1078
1079 numfluxB[3] = 0.5 * (hR * uR * uR * MageB1 + hR * uLB * uLB * MageF1 +
1080 0.5 * g * (hR * hR + hL * hL));
1081 numfluxB[4] = 0.5 * (hR * uR * vR * MageB1 + hR * uLB * vLB * MageF1);
1082 numfluxB[5] = 0.5 * lambdaB * (uLB * hL - uR * hR);
1083
1084 numfluxB[6] = 0.5 * (hR * uR * vR * MageB2 + hR * uLB * vLB * MageF2);
1085 numfluxB[7] = 0.5 * (hR * vR * vR * MageB2 + hR * vLB * vLB * MageF2 +
1086 0.5 * g * (hR * hR + hL * hL));
1087 numfluxB[8] = 0.5 * lambdaB * (vLB * hL - vR * hR);
1088}
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 667 of file MMFSWE.cpp.
670{
671
672 int i, k;
674 int nvariables = 3; // only the dependent variables
675
676 // get temporary arrays
677 Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
678 Array<OneD, Array<OneD, NekDouble>> Bwd(nvariables);
679 Array<OneD, NekDouble> DepthFwd(nTraceNumPoints);
680 Array<OneD, NekDouble> DepthBwd(nTraceNumPoints);
681
682 for (i = 0; i < nvariables; ++i)
683 {
684 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
685 Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
686 }
687
688 // get the physical values at the trace
689 for (i = 0; i < nvariables; ++i)
690 {
691 m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
692
693 // Copy Fwd to Bwd at boundaries
695 }
698
699 // note that we are using the same depth - i.e. the depth is assumed
700 // continuous...
702 {
704 {
705 // rotate the values to the normal direction
706 NekDouble tmpX, tmpY;
707
708 // Fwd[1] = hu^+ = ( h^1^+ (e^1 \cdot n) + h^2^+ ( e^2 \cdot n) )
709 // Fwd[2] = hv^+ = ( h^1^+ (e^1 \cdot n^{\perp}) + h^2^+ ( e^2 \cdot
710 // n^{\perp}) )
711 // Bwd[1] = hu^- = ( h^1^- (e^1 \cdot n) + h^2^- ( e^2 \cdot n) )
712 // Bwd[2] = hv^- = ( h^1^- (e^1 \cdot n^{\perp}) + h^2^- ( e^2 \cdot
713 // n^{\perp}) )
714
715 for (k = 0; k < nTraceNumPoints; ++k)
716 {
717 tmpX = Fwd[1][k] * m_ncdotMFFwd[0][k] +
718 Fwd[2][k] * m_ncdotMFFwd[1][k];
719 tmpY = Fwd[1][k] * m_nperpcdotMFFwd[0][k] +
720 Fwd[2][k] * m_nperpcdotMFFwd[1][k];
721 Fwd[1][k] = tmpX;
722 Fwd[2][k] = tmpY;
723
724 tmpX = Bwd[1][k] * m_ncdotMFBwd[0][k] +
725 Bwd[2][k] * m_ncdotMFBwd[1][k];
726 tmpY = Bwd[1][k] * m_nperpcdotMFBwd[0][k] +
727 Bwd[2][k] * m_nperpcdotMFBwd[1][k];
728 Bwd[1][k] = tmpX;
729 Bwd[2][k] = tmpY;
730 }
731
732 // Solve the Riemann problem
733 NekDouble denomFwd, denomBwd;
734 Array<OneD, NekDouble> numfluxF(nvariables);
735 Array<OneD, NekDouble> numfluxB(nvariables);
736
737 NekDouble eF1n, eF2n, eB1n, eB2n;
738 NekDouble eF1t, eF2t, eB1t, eB2t;
739 for (k = 0; k < nTraceNumPoints; ++k)
740 {
741 RiemannSolverHLLC(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
742 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
743 Bwd[2][k], numfluxF, numfluxB);
744
745 // uflux = 1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
746 // \cdot n)(e^1 \cdot n^{\perp} ) ] ( e^2 \cdot n^{\perp} huflux
747 // - e^2 \cdot n hvflux
748 // vflux = -1/[ ( e^1 \cdot n) ( e^2 \cdot n^{\perp} ) - (e^2
749 // \cdot n)(e^1 \cdot n^{\perp} ) ] ( -e^1 \cdot n^{\perp}
750 // huflux + e^1 \cdot n hvflux
751 eF1n = m_ncdotMFFwd[0][k];
752 eF2n = m_ncdotMFFwd[1][k];
753 eB1n = m_ncdotMFBwd[0][k];
754 eB2n = m_ncdotMFBwd[1][k];
755
756 eF1t = m_nperpcdotMFFwd[0][k];
757 eF2t = m_nperpcdotMFFwd[1][k];
758 eB1t = m_nperpcdotMFBwd[0][k];
759 eB2t = m_nperpcdotMFBwd[1][k];
760
761 denomFwd = eF1n * eF2t - eF2n * eF1t;
762 denomBwd = eB1n * eB2t - eB2n * eB1t;
763
764 numfluxFwd[0][k] = numfluxF[0];
765 numfluxFwd[1][k] = (1.0 / denomFwd) *
766 (eF2t * numfluxF[1] - eF2n * numfluxF[2]);
767 numfluxFwd[2][k] =
768 (1.0 / denomFwd) *
769 (-1.0 * eF1t * numfluxF[1] + eF1n * numfluxF[2]);
770
771 numfluxBwd[0][k] = 1.0 * numfluxB[0];
772 numfluxBwd[1][k] = (1.0 / denomBwd) *
773 (eB2t * numfluxB[1] - eB2n * numfluxB[2]);
774 numfluxBwd[2][k] =
775 (1.0 / denomBwd) *
776 (-1.0 * eB1t * numfluxB[1] + eB1n * numfluxB[2]);
777 }
778 }
779 break;
780
784 {
785 Array<OneD, NekDouble> numfluxF(nvariables * (m_shapedim + 1));
786 Array<OneD, NekDouble> numfluxB(nvariables * (m_shapedim + 1));
787
788 for (k = 0; k < nTraceNumPoints; ++k)
789 {
791 {
792 AverageFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
793 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
794 Bwd[2][k], numfluxF, numfluxB);
795 }
796
798 {
799 LaxFriedrichFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
800 Fwd[2][k], Bwd[0][k] + DepthFwd[k],
801 Bwd[1][k], Bwd[2][k], numfluxF, numfluxB);
802 }
803
805 {
806 RusanovFlux(k, Fwd[0][k] + DepthFwd[k], Fwd[1][k],
807 Fwd[2][k], Bwd[0][k] + DepthFwd[k], Bwd[1][k],
808 Bwd[2][k], numfluxF, numfluxB);
809 }
810
811 int indx;
812 NekDouble tmpF0, tmpF1, tmpB0, tmpB1;
813 for (i = 0; i < nvariables; ++i)
814 {
815 indx = i * (m_shapedim + 1);
816
817 tmpF0 = numfluxF[indx] * m_ncdotMFFwd[0][k];
818 tmpF1 = numfluxF[indx + 1] * m_ncdotMFFwd[1][k];
819
820 tmpB0 = numfluxB[indx] * m_ncdotMFBwd[0][k];
821 tmpB1 = numfluxB[indx + 1] * m_ncdotMFBwd[1][k];
822
823 numfluxFwd[i][k] = tmpF0 + tmpF1 + numfluxF[indx + 2];
824 numfluxBwd[i][k] = tmpB0 + tmpB1 + numfluxB[indx + 2];
825 }
826 }
827 }
828 break;
829
830 default:
831 {
833 }
834 break;
835 }
836}
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:961
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:838
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:1015
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:1090
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:734
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 2929 of file MMFSWE.cpp.
2930{
2932
2933 // h = \eta + d
2935 m_fields[0]->UpdatePhys(), 1);
2936
2937 // hu = h * u
2939 m_fields[1]->UpdatePhys(), 1);
2940
2941 // hv = h * v
2943 m_fields[2]->UpdatePhys(), 1);
2944}
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Nektar::SolverUtils::EquationSystem::m_fields, Vmath::Vadd(), and Vmath::Vmul().
Referenced by DoOdeProjection(), v_DoInitialise(), and v_DoSolve().
◆ PrimitiveToConservative() [2/2]
|
protected |
Definition at line 2871 of file MMFSWE.cpp.
2874{
2875
2877
2878 if (physin[0].get() == physout[0].get())
2879 {
2880 // copy indata and work with tmp array
2881 Array<OneD, Array<OneD, NekDouble>> tmp(3);
2882 for (int i = 0; i < 3; ++i)
2883 {
2884 // deep copy
2885 tmp[i] = Array<OneD, NekDouble>(nq);
2886 Vmath::Vcopy(nq, physin[i], 1, tmp[i], 1);
2887 }
2888
2889 // h = \eta + d
2891
2892 // hu = h * u
2893 Vmath::Vmul(nq, physout[0], 1, tmp[1], 1, physout[1], 1);
2894
2895 // hv = h * v
2896 Vmath::Vmul(nq, physout[0], 1, tmp[2], 1, physout[2], 1);
2897 }
2898 else
2899 {
2900 // h = \eta + d
2902
2903 // hu = h * u
2904 Vmath::Vmul(nq, physout[0], 1, physin[1], 1, physout[1], 1);
2905
2906 // hv = h * v
2907 Vmath::Vmul(nq, physout[0], 1, physin[2], 1, physout[2], 1);
2908 }
2909}
References Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_depth, Vmath::Vadd(), Vmath::Vcopy(), and Vmath::Vmul().
◆ RiemannSolverHLLC()
|
protected |
Definition at line 838 of file MMFSWE.cpp.
843{
845
848
849 NekDouble uRF, vRF, uLB, vLB;
850 NekDouble hstarF, hstarB;
851
852 // Temporary assignments
853 uRF = uR;
854 vRF = vR;
855 uLB = uL;
856 vLB = vL;
857
858 // the two-rarefaction wave assumption
859 hstarF = 0.5 * (cL + cR) + 0.25 * (uL - uRF);
860 hstarF *= hstarF;
861 hstarF *= (1.0 / g);
862
863 hstarB = 0.5 * (cL + cR) + 0.25 * (uLB - uR);
864 hstarB *= hstarB;
865 hstarB *= (1.0 / g);
866
867 NekDouble hfluxF, hufluxF, hvfluxF;
868 NekDouble hfluxB, hufluxB, hvfluxB;
869 Computehhuhvflux(hL, uL, vL, hR, uRF, vRF, hstarF, hfluxF, hufluxF,
870 hvfluxF);
871 Computehhuhvflux(hL, uLB, vLB, hR, uR, vR, hstarB, hfluxB, hufluxB,
872 hvfluxB);
873
874 numfluxF[0] = hfluxF;
875 numfluxF[1] = hufluxF;
876 numfluxF[2] = hvfluxF;
877
878 numfluxB[0] = hfluxB;
879 numfluxB[1] = hufluxB;
880 numfluxB[2] = hvfluxB;
881}
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:883
References Computehhuhvflux(), m_g, and tinysimd::sqrt().
Referenced by NumericalSWEFlux().
◆ RossbyWave()
|
protected |
Definition at line 2592 of file MMFSWE.cpp.
2593{
2595 NekDouble uhat, vhat;
2596 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2597 NekDouble x0j, x1j, x2j;
2598 NekDouble Ath, Bth, Cth, tmp;
2599
2600 Array<OneD, NekDouble> eta(nq, 0.0);
2601 Array<OneD, NekDouble> u(nq, 0.0);
2602 Array<OneD, NekDouble> v(nq, 0.0);
2603
2604 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2606 {
2607 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2608 }
2609
2610 Array<OneD, NekDouble> x(nq);
2611 Array<OneD, NekDouble> y(nq);
2612 Array<OneD, NekDouble> z(nq);
2613
2615
2616 NekDouble R = 4.0;
2617 NekDouble cos2theta, cosRtheta, cos6theta, cos2Rtheta, cosRm1theta;
2618 NekDouble cos2phi, cos4phi, sin4phi, cos8phi;
2619
2620 // disturbancees of Rossby-Haurwitz Wave
2621 NekDouble x0d, y0d, z0d, phi0, theta0;
2622 // NekDouble rad_earth = 6.37122 * 1000000;
2623
2624 phi0 = 40.0 * m_pi / 180.0;
2625 theta0 = 50.0 * m_pi / 180.0;
2626
2627 x0d = cos(phi0) * cos(theta0);
2628 y0d = sin(phi0) * cos(theta0);
2629 z0d = sin(theta0);
2630
2631 for (int j = 0; j < nq; ++j)
2632 {
2633 x0j = x[j];
2634 x1j = y[j];
2635 x2j = z[j];
2636
2637 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2638 cos_theta);
2639
2640 // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2641 // \sin \alpha + \sin \theta \cos \alpha )^2
2642
2643 // tmp = cos^{2R} \theta, R = 4;
2644 cos2theta = cos_theta * cos_theta;
2645 cosRm1theta = cos_theta * cos2theta;
2646 cosRtheta = cos2theta * cos2theta;
2647 cos6theta = cos2theta * cosRtheta;
2648 cos2Rtheta = cosRtheta * cosRtheta;
2649
2651 Ath = tmp * cos2theta +
2653 ((R + 1.0) * cosRtheta + (2 * R * R - R - 2.0) * cos2theta -
2654 2 * R * R);
2655
2657 Bth = tmp * cosRtheta *
2658 ((R * R + 2 * R + 2) - (R + 1.0) * (R + 1.0) * cos2theta);
2659
2660 Cth =
2662
2663 // cos (2 \phi) = 2 * cos \phi * cos \phi - 1.0
2664 cos2phi = 2.0 * cos_varphi * cos_varphi - 1.0;
2665 cos4phi = 2.0 * cos2phi * cos2phi - 1.0;
2666 cos8phi = 2.0 * cos4phi * cos4phi - 1.0;
2667
2668 // sin (2 \phi) = 2 * cos \phi * sin \phi
2669 sin4phi = 4.0 * sin_varphi * cos_varphi * cos2phi;
2670
2672
2673 // disturbances is added
2675 {
2676 eta[j] = eta[j] *
2677 (1.0 + (1.0 / 40.0) * (x0j * x0d + x1j * y0d + x2j * z0d));
2678 }
2679
2680 // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2681 // || e^2 ||^2
2682 uhat = m_angfreq * cos_theta +
2683 m_K * cosRm1theta * (R * sin_theta * sin_theta - cos2theta) *
2684 cos4phi;
2685 vhat = -1.0 * m_K * R * cosRm1theta * sin_theta * sin4phi;
2686
2687 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2688 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2689 uvec[2][j] = vhat * cos_theta;
2690 }
2691
2692 // NekDouble etamin, etaaver;
2693 // etamin = Vmath::Vmin(nq, eta, 1)*rad_earth;
2694 // etaaver = Vmath::Vsum(nq, eta, 1)/nq*rad_earth;
2695
2696 // Projection of u onto the tangent plane with conserving the mag. of the
2697 // velocity.
2698 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2699
2701 {
2702 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2703 }
2704
2705 // u is projected on the tangent plane with preserving its length
2706 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2707
2708 // Change it to the coordinate of moving frames
2709 // CartesianToMovingframes(0,uvecproj,u);
2710 // CartesianToMovingframes(1,uvecproj,v);
2711
2712 u = CartesianToMovingframes(uvec, 0);
2713 v = CartesianToMovingframes(uvec, 1);
2714
2716 {
2717 case (0):
2718 {
2719 outfield = eta;
2720 }
2721 break;
2722
2723 case (1):
2724 {
2725 outfield = u;
2726 }
2727 break;
2728
2729 case (2):
2730 {
2731 outfield = v;
2732 }
2733 break;
2734 }
2735}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_angfreq, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, m_K, m_Omega, Nektar::SolverUtils::MMFSystem::m_pi, m_RossbyDisturbance, Nektar::SolverUtils::EquationSystem::m_spacedim, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ RusanovFlux()
|
protected |
Definition at line 1090 of file MMFSWE.cpp.
1094{
1095 int nvariables = 3;
1096 NekDouble MageF1, MageF2, MageB1, MageB2;
1097 NekDouble eF1_cdot_eB1, eF1_cdot_eB2;
1098 NekDouble eF2_cdot_eB1, eF2_cdot_eB2;
1099
1101 NekDouble uRF, vRF, uLB, vLB;
1102 NekDouble velL, velR;
1103
1104 Array<OneD, NekDouble> EigF(nvariables);
1105 Array<OneD, NekDouble> EigB(nvariables);
1106
1107 // Compute Magnitude and Dot product of moving frames for the index
1108 ComputeMagAndDot(index, MageF1, MageF2, MageB1, MageB2, eF1_cdot_eB1,
1109 eF1_cdot_eB2, eF2_cdot_eB1, eF2_cdot_eB2);
1110
1111 // Get the velocity in the normal to the edge
1114
1115 NekDouble SL, SR;
1116 SL = fabs(velL) + sqrt(g * hL);
1117 SR = fabs(velR) + sqrt(g * hR);
1118
1119 NekDouble S;
1120 if (SL > SR)
1121 {
1122 S = SL;
1123 }
1124 else
1125 {
1126 S = SR;
1127 }
1128
1129 // uRF = uR component in moving frames e^{Fwd}
1130 // vRF = vR component in moving frames e^{Fwd}
1131 uRF = (uR * eF1_cdot_eB1 + vR * eF1_cdot_eB2) / MageF1;
1132 vRF = (uR * eF2_cdot_eB1 + vR * eF2_cdot_eB2) / MageF2;
1133
1134 numfluxF[0] = 0.5 * (hL * uL + hR * uRF);
1135 numfluxF[1] = 0.5 * (hL * vL + hR * vRF);
1136 numfluxF[2] = 0.5 * S * (hL - hR);
1137
1138 numfluxF[3] =
1139 0.5 * (hL * uL * uL + hR * uRF * uRF + 0.5 * g * (hL * hL + hR * hR));
1140 numfluxF[4] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1141 numfluxF[5] = 0.5 * S * (uL * hL - uRF * hR);
1142
1143 numfluxF[6] = 0.5 * (hL * uL * vL + hR * uRF * vRF);
1144 numfluxF[7] =
1145 0.5 * (hL * vL * vL + hR * vRF * vRF + 0.5 * g * (hL * hL + hR * hR));
1146 numfluxF[8] = 0.5 * S * (vL * hL - vRF * hR);
1147
1148 // uLB = uL component in moving frames e^{Bwd}
1149 // vLB = vL component in moving frames e^{Bwd}
1150 uLB = (uL * eF1_cdot_eB1 + vL * eF2_cdot_eB1) / MageB1;
1151 vLB = (uL * eF1_cdot_eB2 + vL * eF2_cdot_eB2) / MageB2;
1152
1153 numfluxB[0] = 0.5 * (hR * uR + hR * uLB);
1154 numfluxB[1] = 0.5 * (hR * vR + hR * vLB);
1155 numfluxB[2] = 0.5 * S * (hL - hR);
1156
1157 numfluxB[3] =
1158 0.5 * (hR * uR * uR + hR * uLB * uLB + 0.5 * g * (hR * hR + hL * hL));
1159 numfluxB[4] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1160 numfluxB[5] = 0.5 * S * (uLB * hL - uR * hR);
1161
1162 numfluxB[6] = 0.5 * (hR * uR * vR + hR * uLB * vLB);
1163 numfluxB[7] =
1164 0.5 * (hR * vR * vR + hR * vLB * vLB + 0.5 * g * (hR * hR + hL * hL));
1165 numfluxB[8] = 0.5 * S * (vLB * hL - vR * hR);
1166}
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 1382 of file MMFSWE.cpp.
1384{
1385
1387 int cnt = 0;
1388
1389 // loop over Boundary Regions
1391 {
1392
1393 // Zonal Boundary Condition
1395 {
1397 {
1399 }
1401 {
1402 // ZonalBoundary2D(n,cnt,inarray);
1403 }
1404 }
1405
1406 // Wall Boundary Condition
1408 {
1410 {
1412 }
1414 {
1415 WallBoundary2D(n, cnt, inarray);
1416 }
1417 }
1418
1419 // Time Dependent Boundary Condition (specified in meshfile)
1421 "eTimeDependent")
1422 {
1423 for (int i = 0; i < nvariables; ++i)
1424 {
1425 m_fields[i]->EvaluateBoundaryConditions(time);
1426 }
1427 }
1428 cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
1429 }
1430}
void WallBoundary2D(int bcRegion, int cnt, Array< OneD, Array< OneD, NekDouble > > &physarray)
Definition: MMFSWE.cpp:1432
References ASSERTL0, Nektar::SolverUtils::EquationSystem::m_expdim, Nektar::SolverUtils::EquationSystem::m_fields, and WallBoundary2D().
Referenced by DoOdeProjection().
◆ SteadyZonalFlow()
|
protected |
Definition at line 2052 of file MMFSWE.cpp.
2054{
2056 NekDouble uhat, vhat;
2057 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2058 NekDouble x0j, x1j, x2j, tmp;
2059
2060 Array<OneD, NekDouble> eta(nq, 0.0);
2061 Array<OneD, NekDouble> u(nq, 0.0);
2062 Array<OneD, NekDouble> v(nq, 0.0);
2063
2064 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2066 {
2067 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2068 }
2069
2070 Array<OneD, NekDouble> x(nq);
2071 Array<OneD, NekDouble> y(nq);
2072 Array<OneD, NekDouble> z(nq);
2073
2075
2076 for (int j = 0; j < nq; ++j)
2077 {
2078 x0j = x[j];
2079 x1j = y[j];
2080 x2j = z[j];
2081
2082 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2083 cos_theta);
2084
2085 // H = H_0 - (1/g)*(a \Omega u_0 + 0.5*u_0^2 )*(- \cos \phi \cos \theta
2086 // \sin \alpha + \sin \theta \cos \alpha )^2
2087 tmp = -1.0 * cos_varphi * cos_theta * sin(m_alpha) +
2088 sin_theta * cos(m_alpha);
2090
2091 // u = (\vec{u} \cdot e^1 )/ || e^1 ||^2 , v = (\vec{u} \cdot e^2 )/
2092 // || e^2 ||^2
2094 sin_theta * cos_varphi * sin(m_alpha));
2096
2097 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2098 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2099 uvec[2][j] = vhat * cos_theta;
2100 }
2101
2102 // Projection of u onto the tangent plane with conserving the mag. of the
2103 // velocity.
2104 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2105
2107 {
2108 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2109 }
2110
2111 // u is projected on the tangent plane with preserving its length
2112 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2113
2114 // Change it to the coordinate of moving frames
2115 // CartesianToMovingframes(0,uvecproj,u);
2116 // CartesianToMovingframes(1,uvecproj,v);
2117
2118 u = CartesianToMovingframes(uvec, 0);
2119 v = CartesianToMovingframes(uvec, 1);
2120
2122 {
2123 case (0):
2124 {
2125 outfield = eta;
2126 }
2127 break;
2128
2129 case (1):
2130 {
2131 outfield = u;
2132 }
2133 break;
2134
2135 case (2):
2136 {
2137 outfield = v;
2138 }
2139 break;
2140 }
2141}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_H0, m_Hvar, Nektar::SolverUtils::EquationSystem::m_spacedim, m_u0, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ TestSWE2Dproblem()
|
private |
Definition at line 1988 of file MMFSWE.cpp.
1991{
1993
1994 Array<OneD, NekDouble> x0(nq);
1995 Array<OneD, NekDouble> x1(nq);
1996 Array<OneD, NekDouble> x2(nq);
1997
1998 m_fields[0]->GetCoords(x0, x1, x2);
1999
2000 Array<OneD, NekDouble> eta0(nq);
2001 Array<OneD, NekDouble> u0(nq);
2002 Array<OneD, NekDouble> v0(nq);
2003
2004 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2005
2007 {
2008 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2009 }
2010
2011 for (int i = 0; i < nq; ++i)
2012 {
2013 eta0[i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2014 (1.0 / cosh(0.395 * x0[i]))) *
2015 (3.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2016 exp(-0.5 * x1[i] * x1[i]);
2017 uvec[0][i] = (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2018 (1.0 / cosh(0.395 * x0[i]))) *
2019 (-9.0 + 6.0 * x1[i] * x1[i]) / (4.0) *
2020 exp(-0.5 * x1[i] * x1[i]);
2021 uvec[1][i] = (-2.0 * 0.395 * tanh(0.395 * x0[i])) *
2022 (0.771 * 0.395 * 0.395 * (1.0 / cosh(0.395 * x0[i])) *
2023 (1.0 / cosh(0.395 * x0[i]))) *
2024 (2.0 * x1[i]) * exp(-0.5 * x1[i] * x1[i]);
2025 }
2026
2027 u0 = CartesianToMovingframes(uvec, 0);
2028 v0 = CartesianToMovingframes(uvec, 1);
2029
2031 {
2032 case (0):
2033 {
2034 outfield = eta0;
2035 }
2036 break;
2037
2038 case (1):
2039 {
2040 outfield = u0;
2041 }
2042 break;
2043
2044 case (2):
2045 {
2046 outfield = v0;
2047 }
2048 break;
2049 }
2050}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::m_fields, and Nektar::SolverUtils::EquationSystem::m_spacedim.
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ TestVorticityComputation()
Definition at line 2946 of file MMFSWE.cpp.
2947{
2948 // Construct beta
2949 int i, k;
2950 int n = 1, m = 1;
2952
2953 NekDouble alpha, beta_theta, beta_phi;
2954
2955 NekDouble xp, yp, zp, Re;
2956 NekDouble theta, phi, sin_theta, cos_theta, sin_varphi, cos_varphi;
2957 NekDouble cosntheta3;
2958
2959 NekDouble thetax, thetay, thetaz;
2960 NekDouble phix, phiy, phiz;
2961
2962 Array<OneD, NekDouble> x(nq);
2963 Array<OneD, NekDouble> y(nq);
2964 Array<OneD, NekDouble> z(nq);
2965
2966 Array<OneD, NekDouble> u(nq);
2967 Array<OneD, NekDouble> v(nq);
2968
2969 Array<OneD, NekDouble> vorticitycompt(nq);
2970 Array<OneD, NekDouble> vorticityexact(nq);
2971
2973
2974 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2976 {
2977 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2978 }
2979
2980 // Generate SphericalCoords
2981 for (k = 0; k < nq; ++k)
2982 {
2983 xp = x[k];
2984 yp = y[k];
2985 zp = z[k];
2986
2987 Re = sqrt(xp * xp + yp * yp + zp * zp);
2988
2989 CartesianToSpherical(xp, yp, zp, sin_varphi, cos_varphi, sin_theta,
2990 cos_theta);
2991
2992 alpha = sin_varphi;
2993
2994 theta = atan2(sin_theta, cos_theta);
2995 phi = atan2(sin_varphi, cos_varphi);
2996
2997 cosntheta3 = cos(n * theta) * cos(n * theta) * cos(n * theta);
2998
2999 beta_theta = -4.0 * n * cosntheta3 * cos(m * phi) * sin(n * theta) / Re;
3000 beta_phi = -m * cosntheta3 * sin(m * phi) / Re;
3001
3002 thetax = -1.0 * cos_varphi * sin_theta;
3003 thetay = -1.0 * sin_varphi * sin_theta;
3004 thetaz = cos_theta;
3005
3006 phix = -1.0 * sin_varphi;
3007 phiy = cos_varphi;
3008 phiz = 0.0;
3009
3010 uvec[0][k] = alpha * (beta_theta * thetax + beta_phi * phix);
3011 uvec[1][k] = alpha * (beta_theta * thetay + beta_phi * phiy);
3012 uvec[2][k] = alpha * (beta_theta * thetaz + beta_phi * phiz);
3013
3014 vorticityexact[k] = -4.0 * n / Re / Re * cos_theta * cos_theta *
3015 cos_varphi * cos(m * phi) * sin(n * theta);
3016 }
3017
3018 u = CartesianToMovingframes(uvec, 0);
3019 v = CartesianToMovingframes(uvec, 1);
3020
3021 std::cout << "chi migi1" << std::endl;
3022
3023 ComputeVorticity(u, v, vorticitycompt);
3024
3025 Vmath::Vsub(nq, vorticityexact, 1, vorticitycompt, 1, vorticitycompt, 1);
3026
3029 << std::endl;
3030}
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 2483 of file MMFSWE.cpp.
2486{
2488
2489 NekDouble uhat, vhat;
2490 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2491 NekDouble x0j, x1j, x2j;
2492 NekDouble Ttheta, Tphi;
2493
2494 Array<OneD, NekDouble> eta(nq, 0.0);
2495 Array<OneD, NekDouble> u(nq, 0.0);
2496 Array<OneD, NekDouble> v(nq, 0.0);
2497
2498 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2500 {
2501 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2502 }
2503
2504 Array<OneD, NekDouble> x(nq);
2505 Array<OneD, NekDouble> y(nq);
2506 Array<OneD, NekDouble> z(nq);
2507
2509
2510 int Nint = 1000;
2511 NekDouble dth, intj;
2512 for (int j = 0; j < nq; ++j)
2513 {
2514 x0j = x[j];
2515 x1j = y[j];
2516 x2j = z[j];
2517
2518 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2519 cos_theta);
2520
2521 Ttheta = atan2(sin_theta, cos_theta);
2522 Tphi = atan2(sin_varphi, cos_varphi);
2523
2524 uhat = ComputeUnstableJetuphi(Ttheta);
2525 vhat = 0.0;
2526
2527 // eta = - (1/g) (\Omega u_0 + 0.4 u^2 ) \sin^2 \theta
2528 dth = Ttheta / Nint;
2529 eta[j] = dth * 0.5 *
2531 for (int i = 1; i < Nint - 1; i++)
2532 {
2533 intj = i * dth;
2534 eta[j] = eta[j] + dth * ComputeUnstableJetEta(intj);
2535 }
2536 eta[j] = (-1.0 / m_g) * eta[j];
2537
2538 // Add perturbation
2540 {
2541 eta[j] = eta[j] + m_hbar * cos_theta * exp(-9.0 * Tphi * Tphi) *
2542 exp(-225.0 * (m_pi / 4.0 - Ttheta) *
2543 (m_pi / 4.0 - Ttheta));
2544 }
2545
2546 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2547 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2548 uvec[2][j] = vhat * cos_theta;
2549 }
2550
2551 // Projection of u onto the tangent plane with conserving the mag. of the
2552 // velocity.
2553 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2554
2556 {
2557 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2558 }
2559
2560 // u is projected on the tangent plane with preserving its length
2561 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2562
2563 // Change it to the coordinate of moving frames
2564 // CartesianToMovingframes(0,uvecproj,u);
2565 // CartesianToMovingframes(1,uvecproj,v);
2566
2567 u = CartesianToMovingframes(uvec, 0);
2568 v = CartesianToMovingframes(uvec, 1);
2569
2571 {
2572 case (0):
2573 {
2574 outfield = eta;
2575 }
2576 break;
2577
2578 case (1):
2579 {
2580 outfield = u;
2581 }
2582 break;
2583
2584 case (2):
2585 {
2586 outfield = v;
2587 }
2588 break;
2589 }
2590}
NekDouble ComputeUnstableJetEta(const NekDouble theta)
Definition: MMFSWE.cpp:2737
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), ComputeUnstableJetEta(), ComputeUnstableJetuphi(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::GetTotPoints(), Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_hbar, Nektar::SolverUtils::MMFSystem::m_pi, m_PurturbedJet, Nektar::SolverUtils::EquationSystem::m_spacedim, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ UnsteadyZonalFlow()
|
protected |
Definition at line 2294 of file MMFSWE.cpp.
2296{
2298 NekDouble uhat, vhat;
2299 NekDouble sin_theta, cos_theta, sin_varphi, cos_varphi;
2300 NekDouble x0j, x1j, x2j, tmp;
2301
2302 NekDouble TR, Ttheta;
2303
2304 Array<OneD, NekDouble> eta(nq, 0.0);
2305 Array<OneD, NekDouble> u(nq, 0.0);
2306 Array<OneD, NekDouble> v(nq, 0.0);
2307
2308 Array<OneD, Array<OneD, NekDouble>> uvec(m_spacedim);
2310 {
2311 uvec[i] = Array<OneD, NekDouble>(nq, 0.0);
2312 }
2313
2314 Array<OneD, NekDouble> x(nq);
2315 Array<OneD, NekDouble> y(nq);
2316 Array<OneD, NekDouble> z(nq);
2317
2319
2320 for (int j = 0; j < nq; ++j)
2321 {
2322 x0j = x[j];
2323 x1j = y[j];
2324 x2j = z[j];
2325
2326 CartesianToSpherical(x0j, x1j, x2j, sin_varphi, cos_varphi, sin_theta,
2327 cos_theta);
2328
2329 // \eta = ( - ( u_0 ( - T_R sin \alpha \cos \theta + \cos \alpha \sin
2330 // \theta ) + \Omega \sin \theta )^2 + \Omega \sin \theta )^2 + (\Omega
2331 // \sin \theta )^2
2332 TR =
2334 Ttheta =
2337
2341 eta[j] = 0.5 * eta[j] / m_g;
2342
2343 // u = u_0*(TR*\sin \alpha * \sin \theta + \cos \alpha * \cos \theta
2344 // v = - u_0 Ttheta * \sin \alpha
2345 uhat =
2348
2349 uvec[0][j] = -1.0 * uhat * sin_varphi - vhat * sin_theta * cos_varphi;
2350 uvec[1][j] = uhat * cos_varphi - vhat * sin_theta * sin_varphi;
2351 uvec[2][j] = vhat * cos_theta;
2352 }
2353
2354 // Projection of u onto the tangent plane with conserving the mag. of the
2355 // velocity.
2356 Array<OneD, Array<OneD, NekDouble>> uvecproj(m_spacedim);
2357
2359 {
2360 uvecproj[i] = Array<OneD, NekDouble>(nq, 0.0);
2361 }
2362
2363 // u is projected on the tangent plane with preserving its length
2364 // GramSchumitz(m_surfaceNormal, uvec, uvecproj, true);
2365
2366 // Change it to the coordinate of moving frames
2367 // CartesianToMovingframes(0,uvecproj,u);
2368 // CartesianToMovingframes(1,uvecproj,v);
2369
2370 u = CartesianToMovingframes(uvec, 0);
2371 v = CartesianToMovingframes(uvec, 1);
2372
2374 {
2375 case (0):
2376 {
2377 outfield = eta;
2378 }
2379 break;
2380
2381 case (1):
2382 {
2383 outfield = u;
2384 }
2385 break;
2386
2387 case (2):
2388 {
2389 outfield = v;
2390 }
2391 break;
2392 }
2393}
References Nektar::SolverUtils::MMFSystem::CartesianToMovingframes(), Nektar::SolverUtils::MMFSystem::CartesianToSpherical(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::GetTotPoints(), m_alpha, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_Omega, Nektar::SolverUtils::EquationSystem::m_spacedim, m_u0, and Nektar::UnitTests::z().
Referenced by v_EvaluateExactSolution(), and v_SetInitialConditions().
◆ v_DoInitialise()
|
overrideprotectedvirtual |
Sets up initial conditions.
Sets the initial conditions.
Reimplemented from Nektar::SolverUtils::UnsteadySystem.
Definition at line 1545 of file MMFSWE.cpp.
1546{
1547 // Compute m_depth and m_Derivdepth
1548 EvaluateWaterDepth();
1549 EvaluateCoriolis();
1550 SetInitialConditions(0.0, dumpInitialConditions);
1551 PrimitiveToConservative();
1552
1553 // transfer the initial conditions to modal values
1555 {
1558 m_fields[i]->UpdateCoeffs());
1559 }
1560}
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:669
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 241 of file MMFSWE.cpp.
242{
244
245 int i, nchk = 1;
246 int nvariables = 0;
249
251 {
252 for (i = 0; i < nfields; ++i)
253 {
254 m_intVariables.push_back(i);
255 }
256 nvariables = nfields;
257 }
258 else
259 {
260 nvariables = m_intVariables.size();
261 }
262
263 // Set up wrapper to fields data storage.
264 Array<OneD, Array<OneD, NekDouble>> fields(nvariables);
265 Array<OneD, Array<OneD, NekDouble>> tmp(nvariables);
266
267 // Order storage to list time-integrated fields first.
268 for (i = 0; i < nvariables; ++i)
269 {
272 }
273
274 // Initialise time integration scheme
276
277 // Check uniqueness of checkpoint output
281 "Only one of IO_CheckTime and IO_CheckSteps "
282 "should be set!");
283
284 LibUtilities::Timer timer;
285 bool doCheckTime = false;
286 int step = 0;
287 NekDouble intTime = 0.0;
288 NekDouble cpuTime = 0.0;
289 NekDouble elapsed = 0.0;
290
291 NekDouble Mass = 0.0, Energy = 0.0, Enstrophy = 0.0, Vorticity = 0.0;
292 Array<OneD, NekDouble> zeta(nq);
293 Array<OneD, Array<OneD, NekDouble>> fieldsprimitive(nvariables);
294 for (int i = 0; i < nvariables; ++i)
295 {
296 fieldsprimitive[i] = Array<OneD, NekDouble>(nq);
297 }
298
300 {
301 timer.Start();
303 timer.Stop();
304
306 elapsed = timer.TimePerTest(1);
307 intTime += elapsed;
308 cpuTime += elapsed;
309
310 // Write out status information
312 m_session->GetComm()->GetRank() == 0)
313 {
314 std::cout << "Steps: " << std::setw(8) << std::left << step + 1
315 << " "
317
318 std::stringstream ss;
319 ss << cpuTime << "s";
320 std::cout << " CPU Time: " << std::setw(8) << std::left << ss.str()
321 << std::endl;
322
323 // Printout Mass, Energy, Enstrophy
324 ConservativeToPrimitive(fields, fieldsprimitive);
325
326 // Vorticity zeta
327 ComputeVorticity(fieldsprimitive[1], fieldsprimitive[2], zeta);
329
330 // Masss = h^*
332
333 // Energy = 0.5*( h^*(u^2 + v^2) + g ( h^2 - h_s^s ) )
334 Energy = (ComputeEnergy(fieldsprimitive[0], fieldsprimitive[1],
335 fieldsprimitive[2]) -
336 m_Energy0) /
337 m_Energy0;
338
339 // Enstrophy = 0.5/h^* ( \mathbf{k} \cdot (\nabla \times \mathbf{v}
340 // ) + f )^2
341 Enstrophy =
342 (ComputeEnstrophy(fieldsprimitive[0], fieldsprimitive[1],
343 fieldsprimitive[2]) -
344 m_Enstrophy0) /
345 m_Enstrophy0;
346
347 std::cout << "dMass = " << std::setw(8) << std::left << Mass << " "
348 << ", dEnergy = " << std::setw(8) << std::left << Energy
349 << " "
350 << ", dEnstrophy = " << std::setw(8) << std::left
351 << Enstrophy << " "
352 << ", dVorticity = " << std::setw(8) << std::left
353 << Vorticity << std::endl
354 << std::endl;
355
356 cpuTime = 0.0;
357 }
358
359 // (h, hu, hv) -> (\eta, u, v)
360 ConservativeToPrimitive();
361
362 // Transform data into coefficient space
363 for (i = 0; i < nvariables; ++i)
364 {
369 }
370
371 // Write out checkpoint files
373 doCheckTime)
374 {
375 Checkpoint_Output(nchk++);
376 doCheckTime = false;
377
378 // (\eta, u, v) -> (h, hu, hv)
379 PrimitiveToConservative();
380 }
381
382 // Step advance
383 ++step;
384 }
385
386 // Print out summary statistics
388 {
390 {
392 << std::endl
394 }
395
397 {
398 std::cout << "Time-integration : " << intTime << "s" << std::endl;
399 }
400 }
401
402 for (i = 0; i < nvariables; ++i)
403 {
406 }
407
408 // (h, hu, hv) -> (\eta, u, v)
409 ConservativeToPrimitive();
410
411 for (i = 0; i < nvariables; ++i)
412 {
414 m_fields[i]->UpdateCoeffs());
415 }
416}
NekDouble ComputeEnstrophy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2182
NekDouble ComputeMass(const Array< OneD, const NekDouble > &eta)
Definition: MMFSWE.cpp:2143
NekDouble ComputeEnergy(const Array< OneD, const NekDouble > &eta, const Array< OneD, const NekDouble > &u, const Array< OneD, const NekDouble > &v)
Definition: MMFSWE.cpp:2153
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:1195
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:89
NekDouble m_cflSafetyFactor
CFL safety factor (comprise between 0 to 1).
Definition: UnsteadySystem.h:97
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
Definition: UnsteadySystem.h:86
std::vector< int > m_intVariables
Definition: UnsteadySystem.h:94
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:49
References tinysimd::abs(), ASSERTL0, Nektar::SolverUtils::EquationSystem::Checkpoint_Output(), ComputeEnergy(), ComputeEnstrophy(), ComputeMass(), ComputeVorticity(), ConservativeToPrimitive(), Nektar::NekConstants::kNekZeroTol, Nektar::SolverUtils::UnsteadySystem::m_cflSafetyFactor, Nektar::SolverUtils::EquationSystem::m_checksteps, Nektar::SolverUtils::EquationSystem::m_checktime, m_Energy0, m_Enstrophy0, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_fintime, Nektar::SolverUtils::EquationSystem::m_infosteps, Nektar::SolverUtils::UnsteadySystem::m_intScheme, Nektar::SolverUtils::UnsteadySystem::m_intVariables, m_Mass0, Nektar::SolverUtils::UnsteadySystem::m_ode, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::EquationSystem::m_steps, Nektar::SolverUtils::EquationSystem::m_time, Nektar::SolverUtils::EquationSystem::m_timestep, m_Vorticity0, PrimitiveToConservative(), Nektar::LibUtilities::Timer::Start(), Nektar::LibUtilities::Timer::Stop(), and Nektar::LibUtilities::Timer::TimePerTest().
◆ v_EvaluateExactSolution()
|
overrideprotectedvirtual |
Reimplemented from Nektar::SolverUtils::EquationSystem.
Definition at line 3218 of file MMFSWE.cpp.
3221{
3223 {
3225 {
3227 }
3228 break;
3229
3231 {
3233 }
3234 break;
3235
3237 {
3239 }
3240 break;
3241
3243 {
3245 }
3246 break;
3247
3249 {
3251 }
3252 break;
3253
3255 {
3257 }
3258 break;
3259
3260 default:
3261 break;
3262 }
3263}
void UnstableJetFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2483
void RossbyWave(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2592
void IsolatedMountainFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2395
void UnsteadyZonalFlow(unsigned int field, const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2294
void TestSWE2Dproblem(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:1988
void SteadyZonalFlow(unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFSWE.cpp:2052
References Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, FilterPython_Function::field, IsolatedMountainFlow(), m_TestType, RossbyWave(), SteadyZonalFlow(), TestSWE2Dproblem(), UnstableJetFlow(), and UnsteadyZonalFlow().
Referenced by v_L2Error().
◆ v_GenerateSummary()
|
overrideprotectedvirtual |
Print Summary.
Reimplemented from Nektar::SolverUtils::MMFSystem.
Definition at line 3265 of file MMFSWE.cpp.
3266{
3267 MMFSystem::v_GenerateSummary(s);
3269
3271 {
3273 {
3275 }
3276 break;
3277
3279 {
3281 m_RossbyDisturbance);
3282 }
3283 break;
3284
3286 {
3288 }
3289 break;
3290
3291 default:
3292 break;
3293 }
3294}
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:418
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:1363
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:51
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
Definition: UnsteadySystem.h:109
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:85
References ASSERTL0, Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs(), Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineProjection(), DoOdeProjection(), DoOdeRhs(), Nektar::SolverUtils::eSphere, Nektar::eTestIsolatedMountain, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, m_AddCoriolis, m_AddRotation, m_alpha, m_angfreq, m_en, Nektar::SolverUtils::UnsteadySystem::m_explicitAdvection, Nektar::SolverUtils::EquationSystem::m_fields, m_g, m_H0, m_hbar, m_hs0, m_Hvar, m_K, m_k2, Nektar::SolverUtils::UnsteadySystem::m_ode, m_Omega, Nektar::SolverUtils::MMFSystem::m_pi, m_PurturbedJet, m_RossbyDisturbance, Nektar::SolverUtils::EquationSystem::m_session, Nektar::SolverUtils::MMFSystem::m_surfaceType, m_TestType, m_theta0, m_theta1, m_u0, m_uthetamax, Nektar::SolverUtils::MMFSystem::MMFInitObject(), Nektar::SIZE_TestType, Nektar::TestTypeMap, TestVorticityComputation(), and Nektar::SolverUtils::UnsteadySystem::v_InitObject().
◆ v_L2Error()
|
overrideprotectedvirtual |
Virtual function for the L_2 error computation between fields and a given exact solution.
Compute the error in the L2-norm.
- Parameters
-
field The field to compare. exactsoln The exact solution to compare with. Normalised Normalise L2-error.
- Returns
- Error in the L2-norm.
Reimplemented from Nektar::SolverUtils::EquationSystem.
Definition at line 3032 of file MMFSWE.cpp.
3035{
3037 NekDouble L2error = -1.0;
3038
3040 {
3042 {
3045 }
3046
3048 {
3049 case (0):
3050 {
3051 // I(h) = I (h - h_exact ) / I (h_exact)
3052 Array<OneD, NekDouble> exactsolution(nq);
3054
3055 // exactsoln = u - u_T so that L2 compute u_T
3057
3059 &exactsolution[0], 1, &exactsolution[0], 1);
3060 Vmath::Vabs(nq, exactsolution, 1, exactsolution, 1);
3061
3062 L2error = (m_fields[0]->Integral(exactsolution)) / L2exact;
3063 }
3064 break;
3065
3066 case (1):
3067 {
3068 // I2 (u) = I( (u - u_ext)^2 + (v - v_ext)^2 )^{1/2} / I(
3069 // u_ext^2 + v_ext^2 )^{1/2}
3070 Array<OneD, NekDouble> exactu(nq);
3071 Array<OneD, NekDouble> exactv(nq);
3072 Array<OneD, NekDouble> tmp(nq);
3073
3074 // L2exact = \int (\sqrt{exactu*exactu+exactv*exactv})
3075 NekDouble L2exact;
3078 Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3079 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3080 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3081
3082 L2exact = m_fields[1]->Integral(tmp);
3083
3084 // L2exact = \int
3085 // (\sqrt{(u-exactu)*(u-exactu)+(v-exactv)*(v-exactv)})
3087 &exactu[0], 1);
3089 &exactv[0], 1);
3090 Vmath::Vmul(nq, exactu, 1, exactu, 1, tmp, 1);
3091 Vmath::Vvtvp(nq, exactv, 1, exactv, 1, tmp, 1, tmp, 1);
3092 Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
3093
3094 L2error = (m_fields[1]->Integral(tmp)) / L2exact;
3095 }
3096 break;
3097
3098 case (2):
3099 {
3100 L2error = 0.0;
3101 }
3102 break;
3103
3104 default:
3105 break;
3106 }
3107
3108 if (Normalised == true)
3109 {
3111
3114
3115 L2error = sqrt(L2error * L2error / Vol);
3116 }
3117 }
3118 else
3119 {
3120 Array<OneD, NekDouble> L2INF(2);
3122 L2error = L2INF[0];
3123 }
3124
3125 return L2error;
3126}
void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
Definition: MMFSWE.cpp:3218
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:918
SOLVER_UTILS_EXPORT int GetNpoints()
Definition: EquationSystem.h:769
int m_NumQuadPointsError
Number of Quadrature points used to work out the error.
Definition: EquationSystem.h:437
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.hpp:352
References Nektar::SolverUtils::EquationSystem::ErrorExtraPoints(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::GetNpoints(), Nektar::SolverUtils::EquationSystem::m_comm, Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_NumQuadPointsError, Nektar::SolverUtils::EquationSystem::m_time, Nektar::LibUtilities::ReduceSum, tinysimd::sqrt(), v_EvaluateExactSolution(), Vmath::Vabs(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().
◆ v_LinfError()
|
overrideprotectedvirtual |
Virtual function for the L_inf error computation between fields and a given exact solution.
Compute the error in the L_inf-norm
- Parameters
-
field The field to compare. exactsoln The exact solution to compare with.
- Returns
- Error in the L_inft-norm.
Reimplemented from Nektar::SolverUtils::EquationSystem.
Definition at line 3128 of file MMFSWE.cpp.
3131{
3133
3135 {
3138 }
3139
3141
3142 // Obtain \vec{u} in cartesian coordinate
3143 Array<OneD, NekDouble> x(nq);
3144 Array<OneD, NekDouble> y(nq);
3145 Array<OneD, NekDouble> z(nq);
3146
3148
3150 {
3151 case (0):
3152 {
3153 NekDouble Letaint;
3154
3155 Array<OneD, NekDouble> exactsolution(nq);
3156
3159 exactsolution);
3160
3161 Letaint = Vmath::Vamax(nq, exactsolution, 1);
3162
3164 1, &exactsolution[0], 1);
3165
3167 }
3168 break;
3169
3170 case (1):
3171 {
3172 Array<OneD, NekDouble> exactu(nq);
3173 Array<OneD, NekDouble> exactv(nq);
3174 Array<OneD, NekDouble> tmpu(nq);
3175 Array<OneD, NekDouble> tmpv(nq);
3176 Array<OneD, NekDouble> Lerr(nq);
3177 Array<OneD, NekDouble> uT(nq);
3178
3181
3182 // Compute max[sqrt{(u-uex)^2 + (v-vex)^2}]
3185
3186 Vmath::Vsub(nq, &exactu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3187 Vmath::Vsub(nq, &exactv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3188
3189 Vmath::Vmul(nq, &tmpu[0], 1, &tmpu[0], 1, &tmpu[0], 1);
3190 Vmath::Vmul(nq, &tmpv[0], 1, &tmpv[0], 1, &tmpv[0], 1);
3191
3192 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &Lerr[0], 1);
3193 Vmath::Vsqrt(nq, &Lerr[0], 1, &Lerr[0], 1);
3194
3195 // uT = max[sqrt( u_T^2 + v_T^2 ) ]
3196 Vmath::Vmul(nq, &exactu[0], 1, &exactu[0], 1, &tmpu[0], 1);
3197 Vmath::Vmul(nq, &exactv[0], 1, &exactv[0], 1, &tmpv[0], 1);
3198 Vmath::Vadd(nq, &tmpu[0], 1, &tmpv[0], 1, &uT[0], 1);
3199 Vmath::Vsqrt(nq, &uT[0], 1, &uT[0], 1);
3200
3202 }
3203 break;
3204
3205 case (2):
3206 {
3207 LinfError = 0.0;
3208 }
3209 break;
3210
3211 default:
3212 break;
3213 }
3214
3216}
SOLVER_UTILS_EXPORT NekDouble LinfError(unsigned int field, const Array< OneD, NekDouble > &exactsoln=NullNekDouble1DArray)
Linf error computation.
Definition: EquationSystem.h:594
SOLVER_UTILS_EXPORT void EvaluateExactSolution(int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Evaluates an exact solution.
Definition: EquationSystem.h:677
References Nektar::SolverUtils::EquationSystem::EvaluateExactSolution(), FilterPython_Function::field, Nektar::SolverUtils::EquationSystem::LinfError(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::EquationSystem::m_time, Vmath::Vadd(), Vmath::Vamax(), Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Nektar::UnitTests::z().
◆ v_SetInitialConditions()
|
overrideprotectedvirtual |
Set the physical fields based on a restart file, or a function describing the initial condition given in the session.
- Parameters
-
initialtime Time at which to evaluate the function. dumpInitialConditions Write the initial condition to file?
Reimplemented from Nektar::SolverUtils::EquationSystem.
Definition at line 1791 of file MMFSWE.cpp.
1794{
1796
1798 {
1800 {
1801 Array<OneD, NekDouble> eta0(nq);
1802 Array<OneD, NekDouble> u0(nq);
1803 Array<OneD, NekDouble> v0(nq);
1804
1805 TestSWE2Dproblem(initialtime, 0, eta0);
1806 m_fields[0]->SetPhys(eta0);
1807
1808 TestSWE2Dproblem(initialtime, 1, u0);
1809 m_fields[1]->SetPhys(u0);
1810
1811 TestSWE2Dproblem(initialtime, 2, v0);
1812 m_fields[2]->SetPhys(v0);
1813
1814 // forward transform to fill the modal coeffs
1816 {
1819 m_fields[i]->UpdateCoeffs());
1820 }
1821 }
1822 break;
1823
1825 {
1826 Array<OneD, NekDouble> eta0(nq);
1827 Array<OneD, NekDouble> u0(nq);
1828 Array<OneD, NekDouble> v0(nq);
1829 Array<OneD, NekDouble> zeta0(nq);
1830
1831 SteadyZonalFlow(0, eta0);
1832 m_fields[0]->SetPhys(eta0);
1833
1834 SteadyZonalFlow(1, u0);
1835 m_fields[1]->SetPhys(u0);
1836
1837 SteadyZonalFlow(2, v0);
1838 m_fields[2]->SetPhys(v0);
1839
1840 // ComputeVorticity(u0, v0, zeta0);
1842
1846
1847 // forward transform to fill the modal coeffs
1849 {
1852 m_fields[i]->UpdateCoeffs());
1853 }
1854 }
1855 break;
1856
1858 {
1859 Array<OneD, NekDouble> eta0(nq);
1860 Array<OneD, NekDouble> u0(nq);
1861 Array<OneD, NekDouble> v0(nq);
1862
1863 UnsteadyZonalFlow(0, initialtime, eta0);
1864 m_fields[0]->SetPhys(eta0);
1865
1866 UnsteadyZonalFlow(1, initialtime, u0);
1867 m_fields[1]->SetPhys(u0);
1868
1869 UnsteadyZonalFlow(2, initialtime, v0);
1870 m_fields[2]->SetPhys(v0);
1871
1875
1876 // forward transform to fill the modal coeffs
1878 {
1881 m_fields[i]->UpdateCoeffs());
1882 }
1883 }
1884 break;
1885
1887 {
1888 Array<OneD, NekDouble> eta0(nq);
1889 Array<OneD, NekDouble> u0(nq);
1890 Array<OneD, NekDouble> v0(nq);
1891
1892 IsolatedMountainFlow(0, initialtime, eta0);
1893 m_fields[0]->SetPhys(eta0);
1894
1895 IsolatedMountainFlow(1, initialtime, u0);
1896 m_fields[1]->SetPhys(u0);
1897
1898 IsolatedMountainFlow(2, initialtime, v0);
1899 m_fields[2]->SetPhys(v0);
1900
1904
1905 // forward transform to fill the modal coeffs
1907 {
1910 m_fields[i]->UpdateCoeffs());
1911 }
1912 }
1913 break;
1914
1916 {
1917 Array<OneD, NekDouble> eta0(nq);
1918 Array<OneD, NekDouble> u0(nq);
1919 Array<OneD, NekDouble> v0(nq);
1920
1921 UnstableJetFlow(0, initialtime, eta0);
1922 m_fields[0]->SetPhys(eta0);
1923
1924 UnstableJetFlow(1, initialtime, u0);
1925 m_fields[1]->SetPhys(u0);
1926
1927 UnstableJetFlow(2, initialtime, v0);
1928 m_fields[2]->SetPhys(v0);
1929
1933
1934 // forward transform to fill the modal coeffs
1936 {
1939 m_fields[i]->UpdateCoeffs());
1940 }
1941 }
1942 break;
1943
1945 {
1946 Array<OneD, NekDouble> eta0(nq);
1947 Array<OneD, NekDouble> u0(nq);
1948 Array<OneD, NekDouble> v0(nq);
1949
1950 RossbyWave(0, eta0);
1951 m_fields[0]->SetPhys(eta0);
1952
1953 RossbyWave(1, u0);
1954 m_fields[1]->SetPhys(u0);
1955
1956 RossbyWave(2, v0);
1957 m_fields[2]->SetPhys(v0);
1958
1962
1963 // forward transform to fill the modal coeffs
1965 {
1968 m_fields[i]->UpdateCoeffs());
1969 }
1970 }
1971 break;
1972
1973 default:
1974 break;
1975 }
1976
1977 if (dumpInitialConditions)
1978 {
1979 // dump initial conditions to file
1981 WriteFld(outname);
1982
1984 Checkpoint_Output_Cartesian(outname);
1985 }
1986}
void Checkpoint_Output_Cartesian(std::string outname)
Definition: MMFSWE.cpp:2767
References Checkpoint_Output_Cartesian(), ComputeEnergy(), ComputeEnstrophy(), ComputeMass(), Nektar::eTestIsolatedMountain, Nektar::eTestPlane, Nektar::eTestRossbyWave, Nektar::eTestSteadyZonal, Nektar::eTestUnstableJet, Nektar::eTestUnsteadyZonal, Nektar::SolverUtils::EquationSystem::GetTotPoints(), IsolatedMountainFlow(), m_Energy0, m_Enstrophy0, Nektar::SolverUtils::EquationSystem::m_fields, m_Mass0, Nektar::SolverUtils::EquationSystem::m_sessionName, m_TestType, m_Vorticity0, RossbyWave(), SteadyZonalFlow(), TestSWE2Dproblem(), UnstableJetFlow(), UnsteadyZonalFlow(), and Nektar::SolverUtils::EquationSystem::WriteFld().
◆ WallBoundary2D()
|
protected |
Definition at line 1432 of file MMFSWE.cpp.
1434{
1435
1436 int i;
1438 int nvariables = physarray.size();
1439
1440 // get physical values of the forward trace
1441 Array<OneD, Array<OneD, NekDouble>> Fwd0(nvariables);
1442 for (i = 0; i < nvariables; ++i)
1443 {
1444 Fwd0[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1445 m_fields[i]->ExtractTracePhys(physarray[i], Fwd0[i]);
1446 }
1447
1448 Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
1449 for (i = 0; i < nvariables; ++i)
1450 {
1451 Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1452 m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
1453 }
1454
1455 // Adjust the physical values of the trace to take
1456 // user defined boundaries into account
1457 int e, id1, id2, npts;
1458
1460 ++e)
1461 {
1462 npts = m_fields[0]
1463 ->GetBndCondExpansions()[bcRegion]
1464 ->GetExp(e)
1465 ->GetTotPoints();
1466 id1 = m_fields[0]->GetBndCondExpansions()[bcRegion]->GetPhys_Offset(e);
1467 id2 = m_fields[0]->GetTrace()->GetPhys_Offset(
1468 m_fields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt + e));
1469
1471 {
1472 case 1:
1473 {
1474 // negate the forward flux
1475 Vmath::Neg(npts, &Fwd[1][id2], 1);
1476 }
1477 break;
1478 case 2:
1479 {
1480 Array<OneD, NekDouble> tmp_n(npts);
1481 Array<OneD, NekDouble> tmp_t(npts);
1482
1484 &tmp_n[0], 1);
1486 &tmp_n[0], 1, &tmp_n[0], 1);
1487
1489 &tmp_t[0], 1);
1491 1, &tmp_t[0], 1, &tmp_t[0], 1);
1492
1493 // negate the normal flux
1494 Vmath::Neg(npts, tmp_n, 1);
1495
1496 Array<OneD, NekDouble> denom(npts);
1497 Array<OneD, NekDouble> tmp_u(npts);
1498 Array<OneD, NekDouble> tmp_v(npts);
1499
1500 // denom = (e^1 \cdot n ) (e^2 \cdot t) - (e^2 \cdot n ) (e^1
1501 // \cdot t)
1503 &m_nperpcdotMFFwd[0][id2], 1, &denom[0], 1);
1505 &m_nperpcdotMFFwd[1][id2], 1, &denom[0], 1,
1506 &denom[0], 1);
1507
1509 &tmp_u[0], 1);
1511 &tmp_u[0], 1, &tmp_u[0], 1);
1512 Vmath::Vdiv(npts, &tmp_u[0], 1, &denom[0], 1, &tmp_u[0], 1);
1513
1514 Vmath::Vcopy(npts, &tmp_u[0], 1, &Fwd[1][id2], 1);
1515
1517 &tmp_v[0], 1);
1519 &tmp_v[0], 1, &tmp_v[0], 1);
1520 Vmath::Vdiv(npts, &tmp_v[0], 1, &denom[0], 1, &tmp_v[0], 1);
1521
1522 Vmath::Vcopy(npts, &tmp_v[0], 1, &Fwd[2][id2], 1);
1523 }
1524 break;
1525
1526 default:
1528 }
1529
1530 // copy boundary adjusted values into the boundary expansion
1531 for (i = 0; i < nvariables; ++i)
1532 {
1533 Vmath::Vcopy(npts, &Fwd[i][id2], 1,
1534 &(m_fields[i]
1535 ->GetBndCondExpansions()[bcRegion]
1536 ->UpdatePhys())[id1],
1537 1);
1538 }
1539 }
1540}
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:744
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 472 of file MMFSWE.cpp.
475{
476 int i;
481
482 Array<OneD, Array<OneD, NekDouble>> fluxvector(m_shapedim);
483 Array<OneD, Array<OneD, NekDouble>> physfield(nvariables);
484
486 {
487 fluxvector[i] = Array<OneD, NekDouble>(nq);
488 }
489
490 // InField is Primitive
491 for (i = 0; i < nvariables; ++i)
492 {
493 physfield[i] = InField[i];
494 }
495
496 // Get the ith component of the flux vector in (physical space)
497 // fluxvector[0] = component for e^1 cdot \nabla \varphi
498 // fluxvector[1] = component for e^2 cdot \nabla \varphi
499 Array<OneD, NekDouble> fluxtmp(nq);
500 Array<OneD, NekDouble> tmp(ncoeffs);
501
502 // Compute Divergence Components
503 for (i = 0; i < nvariables; ++i)
504 {
505 GetSWEFluxVector(i, physfield, fluxvector);
506
507 OutField[i] = Array<OneD, NekDouble>(ncoeffs, 0.0);
509 {
510 // Directional derivation with respect to the j'th moving frame
511 // tmp_j = ( \nabla \phi, fluxvector[j] \mathbf{e}^j )
513 fluxvector[j], tmp);
514 Vmath::Vadd(ncoeffs, &tmp[0], 1, &OutField[i][0], 1,
515 &OutField[i][0], 1);
516 }
517 }
518
519 // Numerical Flux
520 Array<OneD, Array<OneD, NekDouble>> numfluxFwd(nvariables);
521 Array<OneD, Array<OneD, NekDouble>> numfluxBwd(nvariables);
522
523 for (i = 0; i < nvariables; ++i)
524 {
525 numfluxFwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
526 numfluxBwd[i] = Array<OneD, NekDouble>(nTracePointsTot);
527 }
528
529 NumericalSWEFlux(physfield, numfluxFwd, numfluxBwd);
530
531 // Evaulate <\phi, \hat{F}\cdot n> - OutField[i]
532 for (i = 0; i < nvariables; ++i)
533 {
534 Vmath::Neg(ncoeffs, OutField[i], 1);
535 m_fields[i]->AddFwdBwdTraceIntegral(numfluxFwd[i], numfluxBwd[i],
536 OutField[i]);
538 }
539}
void NumericalSWEFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numfluxFwd, Array< OneD, Array< OneD, NekDouble > > &numfluxBwd)
Definition: MMFSWE.cpp:667
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:739
References Nektar::SolverUtils::EquationSystem::GetNcoeffs(), Nektar::SolverUtils::EquationSystem::GetNpoints(), GetSWEFluxVector(), Nektar::SolverUtils::EquationSystem::GetTraceNpoints(), Nektar::SolverUtils::EquationSystem::m_fields, Nektar::SolverUtils::MMFSystem::m_movingframes, Nektar::SolverUtils::MMFSystem::m_shapedim, Vmath::Neg(), NumericalSWEFlux(), and Vmath::Vadd().
Referenced by DoOdeRhs().
Friends And Related Function Documentation
◆ MemoryManager< MMFSWE >
|
friend |
Member Data Documentation
◆ className
|
static |
Initial value:
=
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: BasicUtils/NekFactory.hpp:197
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: MMFSWE.h:66
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:99
Name of class.
◆ m_AddCoriolis
|
protected |
Definition at line 91 of file MMFSWE.h.
Referenced by DoOdeRhs(), and v_InitObject().
◆ m_AddRotation
|
protected |
Definition at line 91 of file MMFSWE.h.
Referenced by DoOdeRhs(), and v_InitObject().
◆ m_alpha
|
protected |
Definition at line 94 of file MMFSWE.h.
Referenced by EvaluateCoriolisForZonalFlow(), SteadyZonalFlow(), UnsteadyZonalFlow(), v_GenerateSummary(), and v_InitObject().
◆ m_angfreq
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by RossbyWave(), and v_InitObject().
◆ m_coriolis
Coriolis force.
Definition at line 89 of file MMFSWE.h.
Referenced by AddCoriolis(), ComputeEnstrophy(), and EvaluateCoriolis().
◆ m_depth
Still water depth.
Definition at line 85 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 86 of file MMFSWE.h.
Referenced by AddElevationEffect(), and EvaluateWaterDepth().
◆ m_en
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_Energy0
|
protected |
Definition at line 100 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_Enstrophy0
|
protected |
Definition at line 100 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_g
|
protected |
Definition at line 93 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 94 of file MMFSWE.h.
Referenced by Checkpoint_Output_Cartesian(), ComputeEnergy(), EvaluateWaterDepth(), RossbyWave(), SteadyZonalFlow(), and v_InitObject().
◆ m_hbar
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by UnstableJetFlow(), and v_InitObject().
◆ m_hs0
|
protected |
Definition at line 96 of file MMFSWE.h.
Referenced by EvaluateWaterDepth(), and v_InitObject().
◆ m_Hvar
|
protected |
Definition at line 94 of file MMFSWE.h.
Referenced by SteadyZonalFlow(), and v_InitObject().
◆ m_K
|
protected |
Definition at line 98 of file MMFSWE.h.
Referenced by RossbyWave(), and v_InitObject().
◆ m_k2
|
protected |
Definition at line 95 of file MMFSWE.h.
Referenced by EvaluateWaterDepth(), and v_InitObject().
◆ m_Mass0
|
protected |
Definition at line 100 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
◆ m_Omega
|
protected |
Definition at line 94 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 270 of file MMFSWE.h.
Referenced by UnstableJetFlow(), v_GenerateSummary(), and v_InitObject().
◆ m_RossbyDisturbance
|
private |
Definition at line 269 of file MMFSWE.h.
Referenced by RossbyWave(), v_GenerateSummary(), and v_InitObject().
◆ m_TestType
| TestType Nektar::MMFSWE::m_TestType |
Definition at line 78 of file MMFSWE.h.
Referenced by EvaluateCoriolis(), EvaluateWaterDepth(), v_EvaluateExactSolution(), v_GenerateSummary(), v_InitObject(), and v_SetInitialConditions().
◆ m_theta0
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_theta1
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_u0
|
protected |
Definition at line 94 of file MMFSWE.h.
Referenced by IsolatedMountainFlow(), SteadyZonalFlow(), UnsteadyZonalFlow(), and v_InitObject().
◆ m_uthetamax
|
protected |
Definition at line 97 of file MMFSWE.h.
Referenced by ComputeUnstableJetuphi(), and v_InitObject().
◆ m_veldotMF
◆ m_vellc
◆ m_velocity
Advection velocity.
Definition at line 106 of file MMFSWE.h.
Referenced by ComputeNablaCdotVelocity().
◆ m_Vorticity0
|
protected |
Definition at line 100 of file MMFSWE.h.
Referenced by v_DoSolve(), and v_SetInitialConditions().
