|
27 | 27 |
|
28 | 28 | #include "orca-jedi/geometry/Geometry.h" |
29 | 29 | #include "orca-jedi/state/State.h" |
| 30 | +#include "orca-jedi/increment/Increment.h" |
30 | 31 |
|
31 | 32 | #include "orca-jedi/interpolator/Interpolator.h" |
32 | 33 |
|
@@ -177,6 +178,128 @@ template void Interpolator::executeInterpolation<float>( |
177 | 178 | const std::vector<bool> & mask, |
178 | 179 | std::vector<double>::iterator& iter) const; |
179 | 180 |
|
| 181 | +/// \brief Interpolate from model space to observation space |
| 182 | +/// \param vars Oops variables |
| 183 | +/// \param inc Increment object (input) |
| 184 | +/// \param mask Mask vector in observation space |
| 185 | +/// \param result Result (output) vector in observation space |
| 186 | +void Interpolator::apply(const oops::Variables& vars, const Increment& inc, |
| 187 | + const std::vector<bool> & mask, |
| 188 | + std::vector<double>& result) const { |
| 189 | + // input is inc output is result |
| 190 | + const size_t nvars = vars.size(); |
| 191 | + |
| 192 | + for (size_t j=0; j < nvars; ++j) { |
| 193 | + if (!inc.variables().has(vars[j])) { |
| 194 | + std::stringstream err_stream; |
| 195 | + err_stream << "orcamodel::Interpolator::apply varname \" " |
| 196 | + << "\" " << vars[j] |
| 197 | + << " not found in the model increment." << std::endl; |
| 198 | + err_stream << " add the variable to the increment variables and " |
| 199 | + << "add a mapping from the geometry to that variable." |
| 200 | + << std::endl; |
| 201 | + throw eckit::BadParameter(err_stream.str(), Here()); |
| 202 | + } |
| 203 | + } |
| 204 | + |
| 205 | + const std::vector<size_t> varSizes = |
| 206 | + inc.geometry()->variableSizes(vars); |
| 207 | + size_t nvals = 0; |
| 208 | + for (size_t jvar=0; jvar < nvars; ++jvar) nvals += nlocs_ * varSizes[jvar]; |
| 209 | + result.resize(nvals); |
| 210 | + |
| 211 | + std::size_t out_idx = 0; |
| 212 | + for (size_t jvar=0; jvar < nvars; ++jvar) { |
| 213 | + auto gv_varname = vars[jvar].name(); |
| 214 | + atlas::Field tgt_field = atlasObsFuncSpace_.createField<double>( |
| 215 | + atlas::option::name(gv_varname) | |
| 216 | + atlas::option::levels(varSizes[jvar])); |
| 217 | + interpolator_.execute(inc.incrementFields()[gv_varname], tgt_field); |
| 218 | + auto field_view = atlas::array::make_view<double, 2>(tgt_field); |
| 219 | + atlas::field::MissingValue mv(inc.incrementFields()[gv_varname]); |
| 220 | + bool has_mv = static_cast<bool>(mv); |
| 221 | + for (std::size_t klev=0; klev < varSizes[jvar]; ++klev) { |
| 222 | + for (std::size_t iloc=0; iloc < nlocs_; iloc++) { |
| 223 | + if (has_mv && mv(field_view(iloc, klev))) { |
| 224 | + result[out_idx] = util::missingValue<double>(); |
| 225 | + } else { |
| 226 | + result[out_idx] = field_view(iloc, klev); |
| 227 | + } |
| 228 | + ++out_idx; |
| 229 | + } |
| 230 | + } |
| 231 | + } |
| 232 | +} |
| 233 | + |
| 234 | +/// \brief Interpolate from observation space to model space |
| 235 | +/// \param vars Oops variables |
| 236 | +/// \param inc Increment object (output) |
| 237 | +/// \param mask Mask (observation space) vector |
| 238 | +/// \param resultin Values (observation space) vector (input) |
| 239 | +void Interpolator::applyAD(const oops::Variables& vars, Increment& inc, |
| 240 | + const std::vector<bool> & mask, |
| 241 | + const std::vector<double> & resultin) const |
| 242 | +{ |
| 243 | + // input is resultin output is inc |
| 244 | + |
| 245 | + oops::Log::trace() << "orcamodel::Interpolator::applyAD start " |
| 246 | + << std::endl; |
| 247 | + |
| 248 | + const size_t nvars = vars.size(); |
| 249 | + |
| 250 | + for (size_t j=0; j < nvars; ++j) { |
| 251 | + if (!inc.variables().has(vars[j])) { |
| 252 | + std::stringstream err_stream; |
| 253 | + err_stream << "orcamodel::Interpolator::apply varname \" " |
| 254 | + << "\" " << vars[j] |
| 255 | + << " not found in the model increment." << std::endl; |
| 256 | + err_stream << " add the variable to the increment variables and " |
| 257 | + << "add a mapping from the geometry to that variable." |
| 258 | + << std::endl; |
| 259 | + throw eckit::BadParameter(err_stream.str(), Here()); |
| 260 | + } |
| 261 | + } |
| 262 | + |
| 263 | + const std::vector<size_t> varSizes = |
| 264 | + inc.geometry()->variableSizes(vars); |
| 265 | + size_t nvals = 0; |
| 266 | + |
| 267 | + for (size_t jvar=0; jvar < nvars; ++jvar) nvals += nlocs_ * varSizes[jvar]; |
| 268 | + |
| 269 | + std::size_t out_idx = 0; |
| 270 | + for (size_t jvar=0; jvar < nvars; ++jvar) { |
| 271 | + auto gv_varname = vars[jvar].name(); |
| 272 | + auto tgt_field = atlasObsFuncSpace_.createField<double>( |
| 273 | + atlas::option::name(gv_varname) | |
| 274 | + atlas::option::levels(varSizes[jvar])); |
| 275 | + |
| 276 | +// Copying observation array vector to an atlas observation field (tgt_field) |
| 277 | + auto field_view = atlas::array::make_view<double, 2>(tgt_field); |
| 278 | + atlas::field::MissingValue mv(inc.incrementFields()[gv_varname]); |
| 279 | + bool has_mv = static_cast<bool>(mv); |
| 280 | + |
| 281 | + for (std::size_t klev=0; klev < varSizes[jvar]; ++klev) { |
| 282 | + for (std::size_t iloc=0; iloc < nlocs_; iloc++) { |
| 283 | + if (has_mv && mv(field_view(iloc, klev))) { |
| 284 | + field_view(iloc, klev) = util::missingValue<double>(); |
| 285 | + } else { |
| 286 | + field_view(iloc, klev) = resultin[out_idx]; |
| 287 | + } |
| 288 | + ++out_idx; |
| 289 | + } |
| 290 | + } |
| 291 | + |
| 292 | +// halo exchange update ghost points |
| 293 | + std::shared_ptr<const Geometry> geom = inc.geometry(); |
| 294 | + geom->functionSpace().haloExchange(inc.incrementFields()[gv_varname]); |
| 295 | + |
| 296 | + interpolator_.execute_adjoint(inc.incrementFields()[gv_varname], tgt_field); |
| 297 | + } // jvar |
| 298 | + |
| 299 | + oops::Log::trace() << "orcamodel::Interpolator::applyAD done " |
| 300 | + << std::endl; |
| 301 | +} |
| 302 | + |
180 | 303 | void Interpolator::print(std::ostream & os) const { |
181 | 304 | os << "orcamodel::Interpolator: " << std::endl; |
182 | 305 | os << " Obs function space " << atlasObsFuncSpace_ << std::endl; |
|
0 commit comments