DiracKernel for transient misfit idaholab#42

modules/optimization/doc/content/source/dirackernels/VectorPointSource.md

@@ -0,0 +1,22 @@
+# VectorPointSource
+
+!syntax parameters /DiracKernels/VectorPointSource
+
+## Overview
+
+A `VectorPointSource` reads in multiple point sources from a [Reporter](Reporters/index.md) or [VectorPostprocessor](VectorPostprocessors/index.md).  The point source values and coordinates are updated as the values are changed.
+
+!alert note
+It is important for the `VectorPointSource` to never use a [VectorPostprocessor](VectorPostprocessors/index.md) with [!param](/VectorPostprocessors/PointValueSampler/contains_complete_history)` = true`, as this can modify the ordering of the coordinates and points.
+
+## Example Input Syntax
+
+An example of a `VectorPointSource` using a [ConstantReporter](/ConstantReporter.md):
+
+!listing vector_point_source.i block=DiracKernels Reporters
+
+!syntax parameters /DiracKernels/VectorPointSource
+
+!syntax inputs /DiracKernels/VectorPointSource
+
+!syntax children /DiracKernels/VectorPointSource

modules/optimization/include/dirackernels/VectorPointSource.h

@@ -0,0 +1,43 @@
+//* 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
+
+// Moose Includes
+#include "DiracKernel.h"
+#include "ReporterInterface.h"
+
+class VectorPointSource : public DiracKernel, public ReporterInterface
+{
+public:
+  static InputParameters validParams();
+  VectorPointSource(const InputParameters & parameters);
+  virtual void addPoints() override;
+
+protected:
+  virtual Real computeQpResidual() override;
+
+private:
+  /// 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;
+  /// The final time when we want to reverse the time index in function evaluation
+  const Real & _reverse_time_end;
+  /// map to associate points with their index into the vpp value
+  std::map<Point, size_t> _point_to_index;
+
+  const std::vector<Real> _empty_vec = {};
+};

modules/optimization/src/dirackernels/VectorPointSource.C

@@ -0,0 +1,104 @@
+//* 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 "VectorPointSource.h"
+#include "MooseUtils.h"
+
+registerMooseObject("isopodApp", VectorPointSource);
+
+InputParameters
+VectorPointSource::validParams()
+{
+  InputParameters params = DiracKernel::validParams();
+  params.addClassDescription("Apply a point load defined by vectors.");
+  params.addParam<ReporterName>(
+      "coord_x",
+      "Vector value containing x-coordinate of points, default is assumed to be all 0s.");
+  params.addParam<ReporterName>(
+      "coord_y",
+      "Vector value containing y-coordinate of points, default is assumed to be all 0s.");
+  params.addParam<ReporterName>(
+      "coord_z",
+      "Vector value containing z-coordinate of points, default is assumed to be all 0s.");
+  params.addParam<ReporterName>("time",
+                                "Vector value containing time, default is assumed to be all 0s.");
+  params.addRequiredParam<ReporterName>("value", "Reporter containing value data.");
+  params.addParam<Real>("reverse_time_end", 0.0, "End time used for reversing the time values.");
+  return params;
+}
+
+VectorPointSource::VectorPointSource(const InputParameters & parameters)
+  : DiracKernel(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")),
+    _reverse_time_end(getParam<Real>("reverse_time_end"))
+{
+}
+
+void
+VectorPointSource::addPoints()
+{
+  // Do some size checks
+  const auto nval = _values.size();
+  if (nval == 0)
+    paramError("value", "Value vector must not be empty.");
+  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,
+               ").");
+
+  _point_to_index.clear();
+  const Real at =
+      MooseUtils::absoluteFuzzyEqual(_reverse_time_end, 0.0) ? _t : _reverse_time_end - _t + _dt;
+  unsigned int id = 0;
+  for (const auto & i : make_range(nval))
+    if (_coordt.empty() || MooseUtils::absoluteFuzzyEqual(at, _coordt[i]))
+    {
+      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];
+      _point_to_index[pt] = i;
+      addPoint(pt, id++);
+    }
+}
+
+Real
+VectorPointSource::computeQpResidual()
+{
+  return -_test[_i][_qp] * _values[_point_to_index[_current_point]];
+}

modules/optimization/test/tests/dirackernels/gold/vector_point_source_out_sample_0000.csv

@@ -0,0 +1,9 @@
+id,u,x,y,z
+0,0,0.25,0.25,0.25
+1,0,0.75,0.25,0.25
+2,0,0.25,0.75,0.25
+3,0,0.75,0.75,0.25
+4,0,0.25,0.25,0.75
+5,0,0.75,0.25,0.75
+6,0,0.25,0.75,0.75
+7,0,0.75,0.75,0.75

modules/optimization/test/tests/dirackernels/gold/vector_point_source_out_sample_0001.csv

@@ -0,0 +1,9 @@
+id,u,x,y,z
+0,0.1963745259433,0.25,0.25,0.25
+1,1.6178575975094,0.75,0.25,0.25
+2,3.0393406690787,0.25,0.75,0.25
+3,4.4608237406472,0.75,0.75,0.25
+4,5.8823068122165,0.25,0.25,0.75
+5,7.3037898837772,0.75,0.25,0.75
+6,8.725272955344,0.25,0.75,0.75
+7,10.146756026918,0.75,0.75,0.75

modules/optimization/test/tests/dirackernels/gold/vector_point_source_out_sample_0002.csv

@@ -0,0 +1,9 @@
+id,u,x,y,z
+0,12.101968747279,0.25,0.25,0.25
+1,13.616679344438,0.75,0.25,0.25
+2,15.1313899416,0.25,0.75,0.25
+3,16.646100538763,0.75,0.75,0.25
+4,18.160811135913,0.25,0.25,0.75
+5,19.675521733076,0.75,0.25,0.75
+6,21.190232330228,0.25,0.75,0.75
+7,22.704942927388,0.75,0.75,0.75

modules/optimization/test/tests/dirackernels/gold/vector_point_source_out_sample_0003.csv

@@ -0,0 +1,9 @@
+id,u,x,y,z
+0,24.897892321335,0.25,0.25,0.25
+1,26.420721184311,0.75,0.25,0.25
+2,27.943550047291,0.25,0.75,0.25
+3,29.466378910269,0.75,0.75,0.25
+4,30.989207773235,0.25,0.25,0.75
+5,32.512036636216,0.75,0.25,0.75
+6,34.034865499185,0.25,0.75,0.75
+7,35.557694362163,0.75,0.75,0.75

modules/optimization/test/tests/dirackernels/gold/vector_point_source_reverse_sample_0000.csv

@@ -0,0 +1,9 @@
+id,u,x,y,z
+0,0,0.25,0.25,0.25
+1,0,0.75,0.25,0.25
+2,0,0.25,0.75,0.25
+3,0,0.75,0.75,0.25
+4,0,0.25,0.25,0.75
+5,0,0.75,0.25,0.75
+6,0,0.25,0.75,0.75
+7,0,0.75,0.75,0.75

modules/optimization/test/tests/dirackernels/gold/vector_point_source_reverse_sample_0001.csv

@@ -0,0 +1,9 @@
+id,u,x,y,z
+0,23.837815789636,0.25,0.25,0.25
+1,25.259298861193,0.75,0.25,0.25
+2,26.680781932768,0.25,0.75,0.25
+3,28.102265004356,0.75,0.75,0.25
+4,29.523748075902,0.25,0.25,0.75
+5,30.945231147459,0.75,0.25,0.75
+6,32.366714219034,0.25,0.75,0.75
+7,33.78819729059,0.75,0.75,0.75

modules/optimization/test/tests/dirackernels/gold/vector_point_source_reverse_sample_0002.csv

@@ -0,0 +1,9 @@
+id,u,x,y,z
+0,13.98160270871,0.25,0.25,0.25
+1,15.496313305868,0.75,0.25,0.25
+2,17.011023903024,0.25,0.75,0.25
+3,18.525734500181,0.75,0.75,0.25
+4,20.040445097349,0.25,0.25,0.75
+5,21.555155694502,0.75,0.25,0.75
+6,23.069866291669,0.25,0.75,0.75
+7,24.584576888826,0.75,0.75,0.75

modules/optimization/test/tests/dirackernels/gold/vector_point_source_reverse_sample_0003.csv

@@ -0,0 +1,9 @@
+id,u,x,y,z
+0,1.5481077061452,0.25,0.25,0.25
+1,3.0709365691205,0.75,0.25,0.25
+2,4.5937654320952,0.25,0.75,0.25
+3,6.1165942950686,0.75,0.75,0.25
+4,7.6394231580537,0.25,0.25,0.75
+5,9.162252021024,0.75,0.25,0.75
+6,10.685080884007,0.25,0.75,0.75
+7,12.207909746981,0.75,0.75,0.75

modules/optimization/test/tests/dirackernels/tests

@@ -0,0 +1,26 @@
+[Tests]
+  [vector_point_source]
+    design = 'VectorPointSource.md'
+    issues = '#42'
+    requirement = 'The system shall support point sources with locations, times, and values given by a vector values with'
+    [forward]
+      type = CSVDiff
+      input = vector_point_source.i
+      csvdiff = 'vector_point_source_out_sample_0000.csv
+                 vector_point_source_out_sample_0001.csv
+                 vector_point_source_out_sample_0002.csv
+                 vector_point_source_out_sample_0003.csv'
+      detail = 'forward time stepping;'
+    []
+    [reverse]
+      type = CSVDiff
+      input = vector_point_source.i
+      cli_args = 'DiracKernels/vpp_point_source/reverse_time_end=0.3 Outputs/file_base=vector_point_source_reverse'
+      csvdiff = 'vector_point_source_reverse_sample_0000.csv
+                 vector_point_source_reverse_sample_0001.csv
+                 vector_point_source_reverse_sample_0002.csv
+                 vector_point_source_reverse_sample_0003.csv'
+      detail = 'reverse time stepping;'
+    []
+  []
+[]

modules/optimization/test/tests/dirackernels/vector_point_source.i

@@ -0,0 +1,103 @@
+[Mesh]
+  [gmg]
+    type = GeneratedMeshGenerator
+    dim = 3
+    nx = 4
+    ny = 4
+    nz = 4
+  []
+[]
+
+[Variables]
+  [u]
+  []
+[]
+
+[Kernels]
+  [diff]
+    type = Diffusion
+    variable = u
+  []
+  [dot]
+    type = TimeDerivative
+    variable = u
+  []
+[]
+
+[DiracKernels]
+  [vpp_point_source]
+    type = VectorPointSource
+    variable = u
+    value = values4D/value
+    coord_x = values4D/coordx
+    coord_y = values4D/coordy
+    coord_z = values4D/coordz
+    time = values4D/time
+  []
+[]
+
+[Reporters]
+  [values4D]
+    type = ConstantReporter
+    real_vector_names = 'coordx coordy coordz time value'
+    real_vector_values = '0.25 0.75 0.25 0.75 0.25 0.75 0.25 0.75
+                          0.25 0.75 0.25 0.75 0.25 0.75 0.25 0.75
+                          0.25 0.75 0.25 0.75 0.25 0.75 0.25 0.75;
+
+                          0.25 0.25 0.75 0.75 0.25 0.25 0.75 0.75
+                          0.25 0.25 0.75 0.75 0.25 0.25 0.75 0.75
+                          0.25 0.25 0.75 0.75 0.25 0.25 0.75 0.75;
+
+                          0.25 0.25 0.25 0.25 0.75 0.75 0.75 0.75
+                          0.25 0.25 0.25 0.25 0.75 0.75 0.75 0.75
+                          0.25 0.25 0.25 0.25 0.75 0.75 0.75 0.75;
+
+                          0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10
+                          0.20 0.20 0.20 0.20 0.20 0.20 0.20 0.20
+                          0.30 0.30 0.30 0.30 0.30 0.30 0.30 0.30;
+
+                          0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00
+                          8.00 9.00 10.0 11.0 12.0 13.0 14.0 15.0
+                          16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0'
+    outputs = none
+  []
+[]
+
+[VectorPostprocessors]
+  [sample]
+    type = PointValueSampler
+    variable = u
+    points = '0.25 0.25 0.25
+              0.75 0.25 0.25
+              0.25 0.75 0.25
+              0.75 0.75 0.25
+              0.25 0.25 0.75
+              0.75 0.25 0.75
+              0.25 0.75 0.75
+              0.75 0.75 0.75'
+    sort_by = id
+    execute_on = 'initial timestep_end'
+  []
+[]
+
+[BCs]
+  [bc]
+    type = DirichletBC
+    variable = u
+    boundary = 'left right top bottom front back'
+    value = 0
+  []
+[]
+
+[Executioner]
+  type = Transient
+  dt = 0.1
+  num_steps = 3
+  solve_type = 'PJFNK'
+  nl_rel_tol = 1e-10
+[]
+
+[Outputs]
+  csv = true
+  execute_on = 'initial timestep_end'
+[]