Skip to content

Commit

Permalink
Update CustomTendencyTerms codes
Browse files Browse the repository at this point in the history 
- Address code suggestions
- Resolve an issue with pm-gpu
  • Loading branch information
@hyungyukang
hyungyukang committed Dec 7, 2024
1 parent 29582a8 commit 11846a0
Show file tree
Hide file tree
Showing 2 changed files with 112 additions and 104 deletions.
108 changes: 108 additions & 0 deletions components/omega/src/ocn/CustomTendencyTerms.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,114 @@ int ManufacturedSolution::init() {

} // end ManufacturedSolution init

//===--------------------------------------------------------------------===/
// Manufactured tendency term for the thickness equation
//===--------------------------------------------------------------------===/
void ManufacturedSolution::ManufacturedThicknessTendency::operator()(
Array2DReal ThicknessTend, const OceanState *State,
const AuxiliaryState *AuxState, int ThickTimeLevel, int VelTimeLevel,
TimeInstant Time) const {

// Get elapsed time since reference time
R8 ElapsedTimeSec;
TimeInterval ElapsedTimeInterval = Time - ReferenceTime;
ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds);

auto *Mesh = HorzMesh::getDefault();
auto NVertLevels = ThicknessTend.extent_int(1);

Array1DReal XCell = Mesh->XCell;
Array1DReal YCell = Mesh->YCell;

OMEGA_SCOPE(LocH0, H0);
OMEGA_SCOPE(LocEta0, Eta0);
OMEGA_SCOPE(LocKx, Kx);
OMEGA_SCOPE(LocKy, Ky);
OMEGA_SCOPE(LocAngFreq, AngFreq);

parallelFor(
{Mesh->NCellsAll, NVertLevels}, KOKKOS_LAMBDA(int ICell, int KLevel) {
R8 X = XCell(ICell);
R8 Y = YCell(ICell);
R8 Phase = LocKx * X + LocKy * Y - LocAngFreq * ElapsedTimeSec;
ThicknessTend(ICell, KLevel) +=
LocEta0 *
(-LocH0 * (LocKx + LocKy) * sin(Phase) - LocAngFreq * cos(Phase) +
LocEta0 * (LocKx + LocKy) * cos(2.0_Real * Phase));
});

} // end void ManufacturedThicknessTendency

//===--------------------------------------------------------------------===/
// Manufactured tendency term for the momentum equation
//===--------------------------------------------------------------------===/
void ManufacturedSolution::ManufacturedVelocityTendency::operator()(
Array2DReal NormalVelTend, const OceanState *State,
const AuxiliaryState *AuxState, int ThickTimeLevel, int VelTimeLevel,
TimeInstant Time) const {

// Get elapsed time since reference time
R8 ElapsedTimeSec;
TimeInterval ElapsedTimeInterval = Time - ReferenceTime;
ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds);

auto *Mesh = HorzMesh::getDefault();
auto NVertLevels = NormalVelTend.extent_int(1);

Array1DReal FEdge = Mesh->FEdge;
Array1DReal XEdge = Mesh->XEdge;
Array1DReal YEdge = Mesh->YEdge;
Array1DReal AngleEdge = Mesh->AngleEdge;

OMEGA_SCOPE(LocGrav, Grav);
OMEGA_SCOPE(LocEta0, Eta0);
OMEGA_SCOPE(LocKx, Kx);
OMEGA_SCOPE(LocKy, Ky);
OMEGA_SCOPE(LocAngFreq, AngFreq);
OMEGA_SCOPE(LocViscDel2, ViscDel2);
OMEGA_SCOPE(LocViscDel4, ViscDel4);
OMEGA_SCOPE(LocVelDiffTendencyEnable, VelDiffTendencyEnable);
OMEGA_SCOPE(LocVelHyperDiffTendencyEnable, VelHyperDiffTendencyEnable);

R8 LocKx2 = LocKx * LocKx;
R8 LocKy2 = LocKy * LocKy;
R8 LocKx4 = LocKx2 * LocKx2;
R8 LocKy4 = LocKy2 * LocKy2;

parallelFor(
{Mesh->NEdgesAll, NVertLevels}, KOKKOS_LAMBDA(int IEdge, int KLevel) {
R8 X = XEdge(IEdge);
R8 Y = YEdge(IEdge);

R8 Phase = LocKx * X + LocKy * Y - LocAngFreq * ElapsedTimeSec;
R8 SourceTerm0 = LocAngFreq * sin(Phase) - 0.5_Real * LocEta0 *
(LocKx + LocKy) *
sin(2.0_Real * Phase);

R8 U = LocEta0 *
((-FEdge(IEdge) + LocGrav * LocKx) * cos(Phase) + SourceTerm0);
R8 V = LocEta0 *
((FEdge(IEdge) + LocGrav * LocKy) * cos(Phase) + SourceTerm0);

// Del2 and del4 source terms
if (LocVelDiffTendencyEnable) {
U += LocViscDel2 * LocEta0 * (LocKx2 + LocKy2) * cos(Phase);
V += LocViscDel2 * LocEta0 * (LocKx2 + LocKy2) * cos(Phase);
}
if (LocVelHyperDiffTendencyEnable) {
U -= LocViscDel4 * LocEta0 *
((LocKx4 + LocKy4 + LocKx2 * LocKy2) * cos(Phase));
V -= LocViscDel4 * LocEta0 *
((LocKx4 + LocKy4 + LocKx2 * LocKy2) * cos(Phase));
}

R8 NormalCompSourceTerm =
cos(AngleEdge(IEdge)) * U + sin(AngleEdge(IEdge)) * V;
NormalVelTend(IEdge, KLevel) += NormalCompSourceTerm;
});

} // end void ManufacturedVelocityTendency

} // end namespace OMEGA

//=-------------------------------------------------------------------------===/
108 changes: 4 additions & 104 deletions components/omega/src/ocn/CustomTendencyTerms.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ namespace OMEGA {
//===-----------------------------------------------------------------------===/
class ManufacturedSolution {

public:
//===--------------------------------------------------------------------===/
// Manufactured tendency term for the thickness equation
//===--------------------------------------------------------------------===/
Expand All @@ -39,40 +40,7 @@ class ManufacturedSolution {

void operator()(Array2DReal ThicknessTend, const OceanState *State,
const AuxiliaryState *AuxState, int ThickTimeLevel,
int VelTimeLevel, TimeInstant Time) {

// Get elapsed time since reference time
R8 ElapsedTimeSec;
TimeInterval ElapsedTimeInterval = Time - ReferenceTime;
ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds);

auto *Mesh = HorzMesh::getDefault();
auto NVertLevels = ThicknessTend.extent_int(1);

Array1DReal XCell = Mesh->XCell;
Array1DReal YCell = Mesh->YCell;

OMEGA_SCOPE(LocXCell, XCell);
OMEGA_SCOPE(LocYCell, YCell);
OMEGA_SCOPE(LocH0, H0);
OMEGA_SCOPE(LocEta0, Eta0);
OMEGA_SCOPE(LocKx, Kx);
OMEGA_SCOPE(LocKy, Ky);
OMEGA_SCOPE(LocAngFreq, AngFreq);

parallelFor(
{Mesh->NCellsAll, NVertLevels},
KOKKOS_LAMBDA(int ICell, int KLevel) {
R8 X = LocXCell(ICell);
R8 Y = LocYCell(ICell);
R8 Phase = LocKx * X + LocKy * Y - LocAngFreq * ElapsedTimeSec;
ThicknessTend(ICell, KLevel) +=
LocEta0 *
(-LocH0 * (LocKx + LocKy) * sin(Phase) -
LocAngFreq * cos(Phase) +
LocEta0 * (LocKx + LocKy) * cos(2.0_Real * Phase));
});
}
int VelTimeLevel, TimeInstant Time) const;
}; // end struct ManufacturedThicknessTendency

//===--------------------------------------------------------------------===/
Expand All @@ -94,78 +62,10 @@ class ManufacturedSolution {

void operator()(Array2DReal NormalVelTend, const OceanState *State,
const AuxiliaryState *AuxState, int ThickTimeLevel,
int VelTimeLevel, TimeInstant Time) {

// Get elapsed time since reference time
R8 ElapsedTimeSec;
TimeInterval ElapsedTimeInterval = Time - ReferenceTime;
ElapsedTimeInterval.get(ElapsedTimeSec, TimeUnits::Seconds);

auto *Mesh = HorzMesh::getDefault();
auto NVertLevels = NormalVelTend.extent_int(1);

Array1DReal FEdge = Mesh->FEdge;
Array1DReal XEdge = Mesh->XEdge;
Array1DReal YEdge = Mesh->YEdge;
Array1DReal AngleEdge = Mesh->AngleEdge;

OMEGA_SCOPE(LocFEdge, FEdge);
OMEGA_SCOPE(LocXEdge, XEdge);
OMEGA_SCOPE(LocYEdge, YEdge);
OMEGA_SCOPE(LocAngleEdge, AngleEdge);
OMEGA_SCOPE(LocGrav, Grav);
OMEGA_SCOPE(LocEta0, Eta0);
OMEGA_SCOPE(LocKx, Kx);
OMEGA_SCOPE(LocKy, Ky);
OMEGA_SCOPE(LocAngFreq, AngFreq);
OMEGA_SCOPE(LocViscDel2, ViscDel2);
OMEGA_SCOPE(LocViscDel4, ViscDel4);
OMEGA_SCOPE(LocVelDiffTendencyEnable, VelDiffTendencyEnable);
OMEGA_SCOPE(LocVelHyperDiffTendencyEnable, VelHyperDiffTendencyEnable);

R8 LocKx2 = LocKx * LocKx;
R8 LocKy2 = LocKy * LocKy;
R8 LocKx4 = LocKx2 * LocKx2;
R8 LocKy4 = LocKy2 * LocKy2;

parallelFor(
{Mesh->NEdgesAll, NVertLevels},
KOKKOS_LAMBDA(int IEdge, int KLevel) {
R8 X = LocXEdge(IEdge);
R8 Y = LocYEdge(IEdge);

R8 Phase = LocKx * X + LocKy * Y - LocAngFreq * ElapsedTimeSec;
R8 SourceTerm0 = LocAngFreq * sin(Phase) -
0.5_Real * LocEta0 * (LocKx + LocKy) *
sin(2.0_Real * Phase);

R8 U = LocEta0 *
((-LocFEdge(IEdge) + LocGrav * LocKx) * cos(Phase) +
SourceTerm0);
R8 V = LocEta0 *
((LocFEdge(IEdge) + LocGrav * LocKy) * cos(Phase) +
SourceTerm0);

// Del2 and del4 source terms
if (LocVelDiffTendencyEnable) {
U += LocViscDel2 * LocEta0 * (LocKx2 + LocKy2) * cos(Phase);
V += LocViscDel2 * LocEta0 * (LocKx2 + LocKy2) * cos(Phase);
}
if (LocVelHyperDiffTendencyEnable) {
U -= LocViscDel4 * LocEta0 *
((LocKx4 + LocKy4 + LocKx2 * LocKy2) * cos(Phase));
V -= LocViscDel4 * LocEta0 *
((LocKx4 + LocKy4 + LocKx2 * LocKy2) * cos(Phase));
}

R8 NormalCompSourceTerm =
cos(LocAngleEdge(IEdge)) * U + sin(LocAngleEdge(IEdge)) * V;
NormalVelTend(IEdge, KLevel) += NormalCompSourceTerm;
});
}
int VelTimeLevel, TimeInstant Time) const;

}; // end struct ManufacturedVelocityTendency

public:
// Instances of manufactured tendencies
ManufacturedThicknessTendency ManufacturedThickTend;
ManufacturedVelocityTendency ManufacturedVelTend;
Expand Down

0 comments on commit 11846a0

Please sign in to comment.