Skip to content

Commit

Permalink
Merge pull request #94 from mmbell/magic7s
Browse files Browse the repository at this point in the history 
Reduce the number of 7s floating around the code.
  • Loading branch information
@cenamiller
cenamiller authored Aug 14, 2024
2 parents e342376 + 3fde384 commit b46153a
Show file tree
Hide file tree
Showing 5 changed files with 48 additions and 39 deletions.
38 changes: 20 additions & 18 deletions src/CostFunction3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@

#include "CostFunction3D.h"
#include "MetObs.h"
#include "VarDriver.h" // added
#include "VarDriver.h"
#include "timing/gptl.h"

#define INDEX(i, j, k, idim, jdim, vdim, var) ((vdim) * ((idim) * ((jdim) * (k) + j) + i) + var)
Expand Down Expand Up @@ -139,8 +139,9 @@ void CostFunction3D::initialize(HashMap* config,
real* bgU, real* obs, ReferenceState* ref)
{
// Initialize number of variables
varDim = 7;
derivDim = 4;
varDim = 7;
derivDim = 4;
obMetaSize = 7;
configHash = config;

/* Set the output path */
Expand Down Expand Up @@ -223,6 +224,7 @@ void CostFunction3D::initialize(HashMap* config,
DJrecip = 1. / DJ;
DKrecip = 1. / DK;


// Adjust the internal, variable domain to include boundaries
adjustInternalDomain(1);

Expand All @@ -246,7 +248,7 @@ void CostFunction3D::initialize(HashMap* config,
// These are local to this one
CTHTd = new real[nState];
stateU = new real[nState];
int64_t vector_size = mObs*(7+varDim*derivDim);
int64_t vector_size = mObs*(obMetaSize+varDim*derivDim);
obsVector = new real[vector_size];
obsData = new real[mObs];
HCq = new real[mObs+nodes];
Expand Down Expand Up @@ -347,14 +349,14 @@ void CostFunction3D::initState(const int iteration)
kFilter->setFilterLengthScale(kFilterScale);

// Set up the Fourier filter
for (int i = 0; i < 7; i++) {
for (int i = 0; i < varDim; i++) {
iMaxWavenumber[i] = -1.0;
jMaxWavenumber[i] = -1.0;
kMaxWavenumber[i] = -1.0;
}
if (configHash->exists("i_max_wavenumber")) {
// Set all the variables to the same filter
for (int i = 0; i < 7; i++) {
for (int i = 0; i < varDim; i++) {
iMaxWavenumber[i] = std::stof((*configHash)["i_max_wavenumber"]);
}
} else {
Expand All @@ -369,7 +371,7 @@ void CostFunction3D::initState(const int iteration)
}

if (configHash->exists("j_max_wavenumber")) {
for (int i = 0; i < 7; i++) {
for (int i = 0; i < varDim; i++) {
jMaxWavenumber[i] = std::stof((*configHash)["j_max_wavenumber"]);
}
} else {
Expand All @@ -383,7 +385,7 @@ void CostFunction3D::initState(const int iteration)
}

if (configHash->exists("k_max_wavenumber")) {
for (int i = 0; i < 7; i++) {
for (int i = 0; i < varDim; i++) {
kMaxWavenumber[i] = std::stof((*configHash)["k_max_wavenumber"]);
}
} else {
Expand Down Expand Up @@ -670,7 +672,7 @@ real CostFunction3D::funcValueAndGradient(real* state, real *gradient)
//#pragma omp parallel for reduction(+:obIP)
#pragma acc parallel loop reduction(+:obIP) private(m)
for (m = 0; m < mObs; m++) {
//int obIndex = m*(7+varDim*derivDim) + 1;
//int obIndex = m*(obMetaSize+varDim*derivDim) + 1;
//obIP += (HCq[m]-innovation[m])*(obsVector[obIndex])*(HCq[m]-innovation[m]);
obIP += (HCq[m]-innovation[m])*(obsData[m])*(HCq[m]-innovation[m]);
}
Expand Down Expand Up @@ -793,7 +795,7 @@ void CostFunction3D::calcInnovation()
#pragma omp parallel for reduction(+:innovationRMS) //[0]
#pragma acc parallel loop reduction(+:innovationRMS) //[0]
for (int m = 0; m < mObs; m++) {
innovation[m] = obsVector[m*(7+varDim*derivDim)] - HCq[m];
innovation[m] = obsVector[m*(obMetaSize+varDim*derivDim)] - HCq[m];
//innovation[m] = obsData[m] - HCq[m];
innovationRMS += (innovation[m]*innovation[m]);
HCq[m] = 0.0;
Expand Down Expand Up @@ -1647,8 +1649,8 @@ void CostFunction3D::obAdjustments() {
//JMD variable-interleave
// Load the obs locally and weight the nonlinear observation operators by interpolated bg fields
for (int m = 0; m < mObs; m++) {
int mi = m*(7+varDim*derivDim);
for (int ob = 0; ob < (7+varDim*derivDim); ob++) {
int mi = m*(obMetaSize+varDim*derivDim);
for (int ob = 0; ob < (obMetaSize+varDim*derivDim); ob++) {
obsVector[mi+ob] = rawObs[mi+ob];
}
real type = obsVector[mi+5];
Expand Down Expand Up @@ -1726,7 +1728,7 @@ void CostFunction3D::obAdjustments() {
}
}
for(int m=0;m<mObs;m++) {
obsData[m]=obsVector[m*(7+varDim*derivDim)+1];
obsData[m]=obsVector[m*(obMetaSize+varDim*derivDim)+1];
}
#pragma acc enter data copyin(obsData)
GPTLstop("CostFunction3D::obAdjustments");
Expand Down Expand Up @@ -2539,7 +2541,7 @@ void CostFunction3D::calcHmatrix()
//#pragma omp parallel for private(m,mi,hi,i,j,k,ii,iis,iie,jj,jjs,jje,kk,kks,kke,ibasis,jbasis,kbasis,iiNode,jjNode,kkNode,iNode,jNode,kNode,var,d,wgt_index,weight,cIndex) //[8.1]
#pragma acc parallel loop vector gang vector_length(32) private(m,mi,hi,iis,iie,jjs,jje,kks,kke,ibasis,jbasis,kbasis,iiNode,jjNode,kkNode,iNode,jNode,kNode,var,d,wgt_index,weight,cIndex) reduction(+:nnz)
for (m = 0; m < mObs; m++) {
mi = m*(7+varDim*derivDim);
mi = m*(obMetaSize+varDim*derivDim);
real i = obsVector[mi+2];
real j = obsVector[mi+3];
real k = obsVector[mi+4];
Expand All @@ -2552,7 +2554,7 @@ void CostFunction3D::calcHmatrix()
hi = 0;
for (var = 0; var < varDim; var++) {
for (d = 0; d < derivDim; d++) {
wgt_index = mi + (7*(d+1)) + var;
wgt_index = mi + (obMetaSize*(d+1)) + var;
if (!obsVector[wgt_index]) continue;
for (iiNode=iis;iiNode<=iie;++iiNode) {
iNode = iiNode;
Expand Down Expand Up @@ -2590,15 +2592,15 @@ void CostFunction3D::calcHmatrix()

for (m = 0; m < mObs; m++) {
IH[m]=hi;
mi = m*(7+varDim*derivDim);
mi = m*(obMetaSize+varDim*derivDim);
i = obsVector[mi+2]; j = obsVector[mi+3]; k = obsVector[mi+4];
ii = (int)((i - iMin)*DIrecip);iis=max(0,ii-1);iie=min(ii+2,iDim-1);
jj = (int)((j - jMin)*DJrecip);jjs=max(0,jj-1);jje=min(jj+2,jDim-1);
kk = (int)((k - kMin)*DKrecip);kks=max(0,kk-1);kke=min(kk+2,kDim-1);
ibasis = 0; jbasis = 0; kbasis = 0;
for (var = 0; var < varDim; var++) {
for (d = 0; d < derivDim; d++) {
wgt_index = mi + (7*(d+1)) + var;
wgt_index = mi + (obMetaSize*(d+1)) + var;
if (!obsVector[wgt_index]) continue;
for (iiNode=iis;iiNode<=iie;++iiNode) {
iNode = iiNode;
Expand Down Expand Up @@ -2654,7 +2656,7 @@ void CostFunction3D::calcHmatrix()
//
// copy H matrix stuff to the GPU Device
#pragma acc enter data copyin(H[:nnz])
std::cout << "Memory usage for [obsVector] (Mbytes): " << sizeof(real)*(mObs*(7+varDim*derivDim))/(1024.0*1024.0) << std::endl;
std::cout << "Memory usage for [obsVector] (Mbytes): " << sizeof(real)*(mObs*(obMetaSize+varDim*derivDim))/(1024.0*1024.0) << std::endl;
std::cout << "Memory usage for [obsData] (Mbytes): " << sizeof(real)*(mObs)/(1024.0*1024.0) << std::endl;
std::cout << "Memory usage for [HCq] (Mbytes): " << sizeof(real)*(mObs+(iDim*jDim*kDim))/(1024.0*1024.0) << std::endl;
std::cout << "Memory usage for [IH,JH] (Mbytes): " << (sizeof(uint64_t)*(mObs+1)+ sizeof(int32_t)*nnz)/(1024.*1024.) << std::endl;
Expand Down
3 changes: 3 additions & 0 deletions src/CostFunction3D.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,8 +138,11 @@ CostFunction3D(const Projection& proj, const int& numObs = 0, const int& stateSi
real* kGammaL;
real* kLL;
real* finalAnalysis;

uint64_t varDim; // NCAR: promoted to 64-bit, since it should auto-promote calculations with it to 64-bit
int derivDim;
int obMetaSize;

real bgError[7];
int iBCL[7], iBCR[7], jBCL[7], jBCR[7], kBCL[7], kBCR[7];
int derivative[4][3];
Expand Down
2 changes: 1 addition & 1 deletion src/CostFunctionRTZ.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -477,7 +477,7 @@ bool CostFunctionRTZ::outputAnalysis(const std::string& suffix, real* Astate)

ostream_iterator<real> od(qcstream, "\t ");
for (int m = 0; m < mObs; m++) {
int64_t mi = m*(7+varDim*derivDim);
int64_t mi = m*(obMetaSize+varDim*derivDim);
real i = obsVector[mi+2];
real j = obsVector[mi+3];
real k = obsVector[mi+4];
Expand Down
41 changes: 21 additions & 20 deletions src/VarDriver3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,9 @@
// Constructor
VarDriver3D::VarDriver3D() : VarDriver()
{
numVars = 7;
numVars = 7; // Number of variables on which to perform the anslysis
numDerivatives = 4;
obMetaSize = 7; // Size of the observation Meta data
bkgdAdapter = NULL;
bgU = NULL;
bgWeights = NULL;
Expand Down Expand Up @@ -2012,10 +2013,10 @@ bool VarDriver3D::preProcessMetObs()
#endif

// Load the observations into a vector
int64_t vector_size = (obVector.size() * (7 + numVars * numDerivatives));
int64_t vector_size = (obVector.size() * (obMetaSize + numVars * numDerivatives));
obs = new real[vector_size];
for (int64_t m=0; m < obVector.size(); m++) {
int64_t n = m * (7 + numVars * numDerivatives);
int64_t n = m * (obMetaSize + numVars * numDerivatives);
Observation ob = obVector.at(m);
obs[n] = ob.getOb();
real invError = ob.getInverseError();
Expand All @@ -2036,7 +2037,7 @@ bool VarDriver3D::preProcessMetObs()
obs[n+6] = ob.getTime();
for (unsigned int var = 0; var < numVars; var++) {
for (unsigned int d = 0; d < numDerivatives; ++d) {
int64_t wgt_index = n + ( 7 * (d + 1)) + var;
int64_t wgt_index = n + ( obMetaSize * (d + 1)) + var;
obs[wgt_index] = ob.getWeight(var, d);
}
}
Expand Down Expand Up @@ -2099,9 +2100,9 @@ bool VarDriver3D::loadPreProcessMetObs()
}

// Load the observations into the vector
obs = new real[obVector.size() * (7 + numVars * numDerivatives)];
obs = new real[obVector.size() * (obMetaSize + numVars * numDerivatives)];
for (int m=0; m < obVector.size(); m++) {
int n = m * (7 + numVars * numDerivatives);
int n = m * (obMetaSize + numVars * numDerivatives);
Observation ob = obVector.at(m);
obs[n] = ob.getOb();
real invError = ob.getInverseError();
Expand All @@ -2122,7 +2123,7 @@ bool VarDriver3D::loadPreProcessMetObs()
obs[n+6] = ob.getTime();
for (unsigned int var = 0; var < numVars; var++) {
for (unsigned int d = 0; d < numDerivatives; ++d) {
int wgt_index = n + (7 * (d + 1)) + var;
int wgt_index = n + (obMetaSize * (d + 1)) + var;
obs[wgt_index] = ob.getWeight(var, d);
}
}
Expand Down Expand Up @@ -2173,7 +2174,7 @@ int VarDriver3D::loadBackgroundObs()
return -1;
}

return bgIn.size() * 7 / 11;
return bgIn.size() * obMetaSize / 11;
}

bool VarDriver3D::adjustBackground()
Expand All @@ -2184,12 +2185,12 @@ bool VarDriver3D::adjustBackground()
START_TIMER(timeab);

// Load the observations into a vector
int numbgObs = bgIn.size() * 7 / 11;
int numbgObs = bgIn.size() * obMetaSize / 11;
if (std::stof(configHash["mc_weight"]) > 0) {
numbgObs += idim * jdim * kdim;
}
bgObs = new real[numbgObs * (7 + numVars * numDerivatives)];
for (unsigned int m = 0; m < numbgObs * (7 + numVars * numDerivatives); m++)
bgObs = new real[numbgObs * (obMetaSize + numVars * numDerivatives)];
for (unsigned int m = 0; m < numbgObs * (obMetaSize + numVars * numDerivatives); m++)
bgObs[m] = 0.;

int p = 0;
Expand Down Expand Up @@ -2247,7 +2248,7 @@ bool VarDriver3D::adjustBackground()
bgObs[p + 5] = -1;
bgObs[p + 6] = obTime;
bgObs[p + 7 + n] = 1.;
p += (7 + numVars * numDerivatives);
p += (obMetaSize + numVars * numDerivatives);
}
}

Expand All @@ -2271,19 +2272,19 @@ bool VarDriver3D::adjustBackground()
bgObs[p + 5] = -1;
bgObs[p + 6] = std::stoi(configHash["ref_time"]);
if (runMode == XYZ) {
bgObs[p + (7 * (1 + 1))] = 1.0;
bgObs[p + (7 * (2 + 1)) + 1] = 1.0;
bgObs[p + (7 * (3 + 1)) + 2] = 1.0;
bgObs[p + (obMetaSize * (1 + 1))] = 1.0;
bgObs[p + (obMetaSize * (2 + 1)) + 1] = 1.0;
bgObs[p + (obMetaSize * (3 + 1)) + 2] = 1.0;
} else if (runMode == RTZ) {
if (i > 0) {
real rInverse = 180.0/(i * Pi);
bgObs[p + 7] = 1.0/i;
bgObs[p + (7 * (1 + 1))] = 1.0;
bgObs[p + (7 * (2 + 1)) + 1] = rInverse;
bgObs[p + (7 * (3 + 1)) + 2] = 1.0;
bgObs[p + obMetaSize] = 1.0/i;
bgObs[p + (obMetaSize * (1 + 1))] = 1.0;
bgObs[p + (obMetaSize * (2 + 1)) + 1] = rInverse;
bgObs[p + (obMetaSize * (3 + 1)) + 2] = 1.0;
}
}
p += (7 + numVars * numDerivatives);
p += (obMetaSize + numVars * numDerivatives);
}
}
}
Expand Down
3 changes: 3 additions & 0 deletions src/VarDriver3D.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,8 @@ class VarDriver3D : public VarDriver
void setGridFlag(bool flag) { fixedGrid = flag; };
bool isFixedGrid() { return fixedGrid == true; };

int obMetaSize;

private:

typedef BSplineBase<real> SplineBase;
Expand Down Expand Up @@ -146,6 +148,7 @@ class VarDriver3D : public VarDriver
int64_t bStateSize;
int64_t uStateSize;


BkgdAdapter *bkgdAdapter;

// Some cost functions (COAMPS) might need access to the sigmas
Expand Down

0 comments on commit b46153a

Please sign in to comment.