Skip to content

Commit

Permalink
Adding objects for time dependent parameters and outputs idaholab#42
Browse files Browse the repository at this point in the history
  • Loading branch information
zachmprince authored and dschwen committed Aug 24, 2022
1 parent 3df597f commit bb3095e
Show file tree
Hide file tree
Showing 4 changed files with 261 additions and 11 deletions.
46 changes: 46 additions & 0 deletions modules/optimization/include/auxkernels/ReporterNearestPointAux.h
Original file line number Diff line number Diff line change
@@ -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<Real> & _coordx;
/// y-coordinates from reporter
const std::vector<Real> & _coordy;
/// z-coordinates from reporter
const std::vector<Real> & _coordz;
/// time-coordinates from reporter
const std::vector<Real> & _coordt;
/// values from reporter
const std::vector<Real> & _values;

/// Data for the current time
std::unordered_map<Point, Real> _current_data;

private:
const std::vector<Real> _empty_vec = {};
};
3 changes: 3 additions & 0 deletions modules/optimization/include/reporters/OptimizationData.h
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ class OptimizationData : public GeneralReporter
std::vector<Real> & _measurement_xcoord;
std::vector<Real> & _measurement_ycoord;
std::vector<Real> & _measurement_zcoord;
std::vector<Real> & _measurement_time;
std::vector<Real> & _measurement_values;
std::vector<Real> & _simulation_values;
std::vector<Real> & _misfit_values;
Expand All @@ -38,4 +39,6 @@ class OptimizationData : public GeneralReporter
void readMeasurementsFromFile();
void readMeasurementsFromInput();
void setSimuilationValuesForTesting(std::vector<Real> & data);

const MooseVariableFieldBase * const _var;
};
154 changes: 154 additions & 0 deletions modules/optimization/src/auxkernels/ReporterNearestPointAux.C
Original file line number Diff line number Diff line change
@@ -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<ReporterName>(
"coord_x",
"Vector reporter value containing x-coordinate of points, default is assumed to be all 0s.");
params.addParam<ReporterName>(
"coord_y",
"Vector reporter value containing y-coordinate of points, default is assumed to be all 0s.");
params.addParam<ReporterName>(
"coord_z",
"Vector reporter value containing z-coordinate of points, default is assumed to be all 0s.");
params.addParam<ReporterName>("time",
"Vector reporter value of time points if 'value' is a vector of "
"vectors with time-dependent data.");
params.addRequiredParam<ReporterName>("value", "Reporter containing value data.");
return params;
}

ReporterNearestPointAux::ReporterNearestPointAux(const InputParameters & parameters)
: AuxKernel(parameters),
ReporterInterface(this),
_coordx(isParamValid("coord_x") ? getReporterValue<std::vector<Real>>("coord_x") : _empty_vec),
_coordy(isParamValid("coord_y") ? getReporterValue<std::vector<Real>>("coord_y") : _empty_vec),
_coordz(isParamValid("coord_z") ? getReporterValue<std::vector<Real>>("coord_z") : _empty_vec),
_coordt(isParamValid("time") ? getReporterValue<std::vector<Real>>("time") : _empty_vec),
_values(getReporterValue<std::vector<Real>>("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<Point, std::vector<std::pair<Real, Real>>> 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<std::pair<Real, Real>> * 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<Real, Real> & a, const std::pair<Real, Real> & 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<Point, Real> & p1, const std::pair<Point, Real> & p2)
{ return (pt - p1.first).norm_sq() < (pt - p2.first).norm_sq(); });
return it->second;
}
69 changes: 58 additions & 11 deletions modules/optimization/src/reporters/OptimizationData.C
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,15 @@

#include "OptimizationData.h"
#include "DelimitedFileReader.h"
#include "SystemBase.h"

registerMooseObject("isopodApp", OptimizationData);

InputParameters
OptimizationData::validParams()
{
InputParameters params = GeneralReporter::validParams();

params.addClassDescription(
"Reporter to hold measurement and simulation data for optimization problems");

Expand All @@ -24,6 +26,8 @@ OptimizationData::validParams()
"Measurement values collected from locations given by measurement_points");
params.addParam<std::vector<Point>>("measurement_points",
"Point locations corresponding to each measurement value");
params.addParam<std::vector<Real>>("measurement_times",
"Times corresponding to each measurement value");

params.addParam<FileName>("measurement_file",
"CSV file with measurement value and coordinates (value, x, y, z).");
Expand All @@ -33,8 +37,11 @@ OptimizationData::validParams()
"file_ycoord", "y", "y coordinate column name from csv file being read in");
params.addParam<std::string>(
"file_zcoord", "z", "z coordinate column name from csv file being read in");
params.addParam<std::string>("file_time", "time", "time column name from csv file being read in");
params.addParam<std::string>(
"file_value", "value", "measurement value column name from csv file being read in");

params.addParam<VariableName>("variable", "Variable to sample at measurement points.");
return params;
}

Expand All @@ -46,11 +53,20 @@ OptimizationData::OptimizationData(const InputParameters & parameters)
declareValueByName<std::vector<Real>>("measurement_ycoord", REPORTER_MODE_REPLICATED)),
_measurement_zcoord(
declareValueByName<std::vector<Real>>("measurement_zcoord", REPORTER_MODE_REPLICATED)),
_measurement_time(
declareValueByName<std::vector<Real>>("measurement_time", REPORTER_MODE_REPLICATED)),
_measurement_values(
declareValueByName<std::vector<Real>>("measurement_values", REPORTER_MODE_REPLICATED)),
_simulation_values(
declareValueByName<std::vector<Real>>("simulation_values", REPORTER_MODE_REPLICATED)),
_misfit_values(declareValueByName<std::vector<Real>>("misfit_values", REPORTER_MODE_REPLICATED))
_misfit_values(
declareValueByName<std::vector<Real>>("misfit_values", REPORTER_MODE_REPLICATED)),
_var(isParamValid("variable")
? &_fe_problem.getVariable(_tid,
getParam<VariableName>("variable"),
Moose::VarKindType::VAR_ANY,
Moose::VarFieldType::VAR_FIELD_STANDARD)
: nullptr)
{
if (isParamValid("measurement_file"))
readMeasurementsFromFile();
Expand All @@ -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
Expand All @@ -69,11 +107,13 @@ OptimizationData::readMeasurementsFromFile()
std::string xName = getParam<std::string>("file_xcoord");
std::string yName = getParam<std::string>("file_ycoord");
std::string zName = getParam<std::string>("file_zcoord");
std::string tName = getParam<std::string>("file_time");
std::string valueName = getParam<std::string>("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<FileName>("measurement_file"));
Expand All @@ -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;
}
}
Expand All @@ -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");
Expand All @@ -130,14 +173,18 @@ OptimizationData::readMeasurementsFromFile()
void
OptimizationData::readMeasurementsFromInput()
{
std::vector<Point> measurement_points = getParam<std::vector<Point>>("measurement_points");
for (auto & p : measurement_points)
for (const auto & p : getParam<std::vector<Point>>("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<std::vector<Real>>("measurement_times");
else
_measurement_time.assign(_measurement_xcoord.size(), 0.0);

if (isParamValid("measurement_values"))
_measurement_values = getParam<std::vector<Real>>("measurement_values");
else
Expand Down

0 comments on commit bb3095e

Please sign in to comment.