diff --git a/modules/optimization/include/auxkernels/ReporterNearestPointAux.h b/modules/optimization/include/auxkernels/ReporterNearestPointAux.h new file mode 100644 index 000000000000..59a1af16f903 --- /dev/null +++ b/modules/optimization/include/auxkernels/ReporterNearestPointAux.h @@ -0,0 +1,46 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "AuxKernel.h" +#include "ReporterInterface.h" + +/** + &ant auxiliary value + */ +class ReporterNearestPointAux : public AuxKernel, public ReporterInterface +{ +public: + static InputParameters validParams(); + + ReporterNearestPointAux(const InputParameters & parameters); + + virtual void subdomainSetup() override; + +protected: + virtual Real computeValue() override; + + /// x-coordinates from reporter + const std::vector & _coordx; + /// y-coordinates from reporter + const std::vector & _coordy; + /// z-coordinates from reporter + const std::vector & _coordz; + /// time-coordinates from reporter + const std::vector & _coordt; + /// values from reporter + const std::vector & _values; + + /// Data for the current time + std::unordered_map _current_data; + +private: + const std::vector _empty_vec = {}; +}; diff --git a/modules/optimization/include/reporters/OptimizationData.h b/modules/optimization/include/reporters/OptimizationData.h index 40639e8c75a9..ff7836f3c5bb 100644 --- a/modules/optimization/include/reporters/OptimizationData.h +++ b/modules/optimization/include/reporters/OptimizationData.h @@ -26,6 +26,7 @@ class OptimizationData : public GeneralReporter std::vector & _measurement_xcoord; std::vector & _measurement_ycoord; std::vector & _measurement_zcoord; + std::vector & _measurement_time; std::vector & _measurement_values; std::vector & _simulation_values; std::vector & _misfit_values; @@ -38,4 +39,6 @@ class OptimizationData : public GeneralReporter void readMeasurementsFromFile(); void readMeasurementsFromInput(); void setSimuilationValuesForTesting(std::vector & data); + + const MooseVariableFieldBase * const _var; }; diff --git a/modules/optimization/src/auxkernels/ReporterNearestPointAux.C b/modules/optimization/src/auxkernels/ReporterNearestPointAux.C new file mode 100644 index 000000000000..4ef50a1e18b2 --- /dev/null +++ b/modules/optimization/src/auxkernels/ReporterNearestPointAux.C @@ -0,0 +1,154 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "ReporterNearestPointAux.h" + +registerMooseObject("isopodApp", ReporterNearestPointAux); + +InputParameters +ReporterNearestPointAux::validParams() +{ + InputParameters params = AuxKernel::validParams(); + params.addClassDescription( + "Assigns variable based on the nearest point to coordinates and values defined by a vector " + "reporter values, interpolates linearly in time with transient data."); + params.addParam( + "coord_x", + "Vector reporter value containing x-coordinate of points, default is assumed to be all 0s."); + params.addParam( + "coord_y", + "Vector reporter value containing y-coordinate of points, default is assumed to be all 0s."); + params.addParam( + "coord_z", + "Vector reporter value containing z-coordinate of points, default is assumed to be all 0s."); + params.addParam("time", + "Vector reporter value of time points if 'value' is a vector of " + "vectors with time-dependent data."); + params.addRequiredParam("value", "Reporter containing value data."); + return params; +} + +ReporterNearestPointAux::ReporterNearestPointAux(const InputParameters & parameters) + : AuxKernel(parameters), + ReporterInterface(this), + _coordx(isParamValid("coord_x") ? getReporterValue>("coord_x") : _empty_vec), + _coordy(isParamValid("coord_y") ? getReporterValue>("coord_y") : _empty_vec), + _coordz(isParamValid("coord_z") ? getReporterValue>("coord_z") : _empty_vec), + _coordt(isParamValid("time") ? getReporterValue>("time") : _empty_vec), + _values(getReporterValue>("value")) +{ +} + +void +ReporterNearestPointAux::subdomainSetup() +{ + // Do some size checks + const std::size_t nval = _values.size(); + if (nval == 0) + paramError("value", "Number of values must be greater than 0."); + else if (!_coordx.empty() && _coordx.size() != nval) + paramError("coord_x", + "Number of x coordinates (", + _coordx.size(), + ") does not match number of values (", + nval, + ")."); + else if (!_coordy.empty() && _coordy.size() != nval) + paramError("coord_y", + "Number of y coordinates (", + _coordy.size(), + ") does not match number of values (", + nval, + ")."); + else if (!_coordz.empty() && _coordz.size() != nval) + paramError("coord_z", + "Number of z coordinates (", + _coordz.size(), + ") does not match number of values (", + nval, + ")."); + else if (!_coordt.empty() && _coordt.size() != nval) + paramError("time", + "Number of times (", + _coordt.size(), + ") does not match number of values (", + nval, + ")."); + + // Find times for each unique coordinate + std::map>> coord_mapping; + for (const auto & i : make_range(nval)) + { + Point pt; + pt(0) = _coordx.empty() ? 0.0 : _coordx[i]; + pt(1) = _coordy.empty() ? 0.0 : _coordy[i]; + pt(2) = _coordz.empty() ? 0.0 : _coordz[i]; + const Real time = _coordt.empty() ? 0.0 : _coordt[i]; + + std::vector> * vec = nullptr; + for (auto & it : coord_mapping) + if (pt.absolute_fuzzy_equals(it.first)) + { + vec = &it.second; + break; + } + if (!vec) + vec = &coord_mapping[pt]; + + vec->emplace_back(time, _values[i]); + std::sort(vec->begin(), + vec->end(), + [](const std::pair & a, const std::pair & b) + { return a.first < b.first; }); + } + + // Linearly interpolate in time + _current_data.clear(); + for (const auto & it : coord_mapping) + { + auto & val = _current_data[it.first]; + const auto & tval = it.second; + if (tval.size() == 1) + val = tval[0].second; + else if (MooseUtils::absoluteFuzzyLessEqual(_t, tval[0].first)) + val = tval[0].second; + else if (MooseUtils::absoluteFuzzyGreaterEqual(_t, tval.back().first)) + val = tval.back().second; + else + for (std::size_t ti = 1; ti < tval.size(); ++ti) + { + if (MooseUtils::absoluteFuzzyEqual(_t, tval[ti].first)) + { + val = tval[ti].second; + break; + } + else if (tval[ti].first > _t) + { + const Real told = tval[ti - 1].first; + const Real tnew = tval[ti].first; + const Real vold = tval[ti - 1].second; + const Real vnew = tval[ti].second; + val = vold + (vnew - vold) * (_t - told) / (tnew - told); + break; + } + } + } +} + +Real +ReporterNearestPointAux::computeValue() +{ + const Point & pt = _q_point[_qp]; + const auto it = + std::min_element(_current_data.begin(), + _current_data.end(), + [&pt](const std::pair & p1, const std::pair & p2) + { return (pt - p1.first).norm_sq() < (pt - p2.first).norm_sq(); }); + return it->second; +} diff --git a/modules/optimization/src/reporters/OptimizationData.C b/modules/optimization/src/reporters/OptimizationData.C index 6c921c78ed95..80cd8a6c9006 100644 --- a/modules/optimization/src/reporters/OptimizationData.C +++ b/modules/optimization/src/reporters/OptimizationData.C @@ -9,6 +9,7 @@ #include "OptimizationData.h" #include "DelimitedFileReader.h" +#include "SystemBase.h" registerMooseObject("isopodApp", OptimizationData); @@ -16,6 +17,7 @@ InputParameters OptimizationData::validParams() { InputParameters params = GeneralReporter::validParams(); + params.addClassDescription( "Reporter to hold measurement and simulation data for optimization problems"); @@ -24,6 +26,8 @@ OptimizationData::validParams() "Measurement values collected from locations given by measurement_points"); params.addParam>("measurement_points", "Point locations corresponding to each measurement value"); + params.addParam>("measurement_times", + "Times corresponding to each measurement value"); params.addParam("measurement_file", "CSV file with measurement value and coordinates (value, x, y, z)."); @@ -33,8 +37,11 @@ OptimizationData::validParams() "file_ycoord", "y", "y coordinate column name from csv file being read in"); params.addParam( "file_zcoord", "z", "z coordinate column name from csv file being read in"); + params.addParam("file_time", "time", "time column name from csv file being read in"); params.addParam( "file_value", "value", "measurement value column name from csv file being read in"); + + params.addParam("variable", "Variable to sample at measurement points."); return params; } @@ -46,11 +53,20 @@ OptimizationData::OptimizationData(const InputParameters & parameters) declareValueByName>("measurement_ycoord", REPORTER_MODE_REPLICATED)), _measurement_zcoord( declareValueByName>("measurement_zcoord", REPORTER_MODE_REPLICATED)), + _measurement_time( + declareValueByName>("measurement_time", REPORTER_MODE_REPLICATED)), _measurement_values( declareValueByName>("measurement_values", REPORTER_MODE_REPLICATED)), _simulation_values( declareValueByName>("simulation_values", REPORTER_MODE_REPLICATED)), - _misfit_values(declareValueByName>("misfit_values", REPORTER_MODE_REPLICATED)) + _misfit_values( + declareValueByName>("misfit_values", REPORTER_MODE_REPLICATED)), + _var(isParamValid("variable") + ? &_fe_problem.getVariable(_tid, + getParam("variable"), + Moose::VarKindType::VAR_ANY, + Moose::VarFieldType::VAR_FIELD_STANDARD) + : nullptr) { if (isParamValid("measurement_file")) readMeasurementsFromFile(); @@ -61,6 +77,28 @@ OptimizationData::OptimizationData(const InputParameters & parameters) void OptimizationData::execute() { + if (!_var) + return; + + // FIXME: This is basically copied from PointValue. + // Implementation can be improved using the functionality in PointSamplerBase, + // but this will require changes in MOOSE to work for reporters. + + const std::size_t nvals = _measurement_xcoord.size(); + _simulation_values.resize(nvals); + _misfit_values.resize(nvals); + + const auto & sys = _var->sys().system(); + const auto vnum = _var->number(); + + for (const auto & i : make_range(nvals)) + if (MooseUtils::absoluteFuzzyEqual(_t, _measurement_time[i])) + { + const Point point(_measurement_xcoord[i], _measurement_ycoord[i], _measurement_zcoord[i]); + const Real val = sys.point_value(vnum, point, false); + _simulation_values[i] = val; + _misfit_values[i] = val - _measurement_values[i]; + } } void @@ -69,11 +107,13 @@ OptimizationData::readMeasurementsFromFile() std::string xName = getParam("file_xcoord"); std::string yName = getParam("file_ycoord"); std::string zName = getParam("file_zcoord"); + std::string tName = getParam("file_time"); std::string valueName = getParam("file_value"); bool found_x = false; bool found_y = false; bool found_z = false; + bool found_t = false; bool found_value = false; MooseUtils::DelimitedFileReader reader(getParam("measurement_file")); @@ -91,26 +131,27 @@ OptimizationData::readMeasurementsFromFile() if (names[i] == xName) { - for (auto && d : data[i]) - _measurement_xcoord.push_back(d); + _measurement_xcoord = data[i]; found_x = true; } else if (names[i] == yName) { - for (auto && d : data[i]) - _measurement_ycoord.push_back(d); + _measurement_ycoord = data[i]; found_y = true; } else if (names[i] == zName) { - for (auto && d : data[i]) - _measurement_zcoord.push_back(d); + _measurement_zcoord = data[i]; found_z = true; } + else if (names[i] == tName) + { + _measurement_time = data[i]; + found_t = true; + } else if (names[i] == valueName) { - for (auto && d : data[i]) - _measurement_values.push_back(d); + _measurement_values = data[i]; found_value = true; } } @@ -122,6 +163,8 @@ OptimizationData::readMeasurementsFromFile() paramError("measurement_file", "Column with name '", yName, "' missing from measurement file"); else if (!found_z) paramError("measurement_file", "Column with name '", zName, "' missing from measurement file"); + else if (!found_t) + _measurement_time.assign(rows, 0); else if (!found_value) paramError( "measurement_file", "Column with name '", valueName, "' missing from measurement file"); @@ -130,14 +173,18 @@ OptimizationData::readMeasurementsFromFile() void OptimizationData::readMeasurementsFromInput() { - std::vector measurement_points = getParam>("measurement_points"); - for (auto & p : measurement_points) + for (const auto & p : getParam>("measurement_points")) { _measurement_xcoord.push_back(p(0)); _measurement_ycoord.push_back(p(1)); _measurement_zcoord.push_back(p(2)); } + if (isParamValid("measurement_times")) + _measurement_time = getParam>("measurement_times"); + else + _measurement_time.assign(_measurement_xcoord.size(), 0.0); + if (isParamValid("measurement_values")) _measurement_values = getParam>("measurement_values"); else