diff --git a/framework/include/auxkernels/SolutionAux.h b/framework/include/auxkernels/SolutionAux.h index 1ab91effa858..147e6ed13bc5 100644 --- a/framework/include/auxkernels/SolutionAux.h +++ b/framework/include/auxkernels/SolutionAux.h @@ -11,7 +11,7 @@ #include "AuxKernel.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** * AuxKernel for reading a solution from file. @@ -41,7 +41,7 @@ class SolutionAux : public AuxKernel virtual Real computeValue() override; /// Reference to the SolutionUserObject storing the solution - const SolutionUserObject & _solution_object; + const SolutionUserObjectBase & _solution_object; /// The variable name of interest std::string _var_name; diff --git a/framework/include/auxkernels/SolutionScalarAux.h b/framework/include/auxkernels/SolutionScalarAux.h index 91b66380c3ed..80e4b33ec975 100644 --- a/framework/include/auxkernels/SolutionScalarAux.h +++ b/framework/include/auxkernels/SolutionScalarAux.h @@ -11,7 +11,7 @@ #include "AuxScalarKernel.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** * AuxScalarKernel for reading a solution from file. @@ -32,7 +32,7 @@ class SolutionScalarAux : public AuxScalarKernel virtual Real computeValue() override; /// Reference to the SolutionUserObject storing the solution - const SolutionUserObject & _solution_object; + const SolutionUserObjectBase & _solution_object; /// The variable name of interest std::string _var_name; diff --git a/framework/include/functions/Axisymmetric2D3DSolutionFunction.h b/framework/include/functions/Axisymmetric2D3DSolutionFunction.h index 1cec3d6984ab..3dece4349f07 100644 --- a/framework/include/functions/Axisymmetric2D3DSolutionFunction.h +++ b/framework/include/functions/Axisymmetric2D3DSolutionFunction.h @@ -11,7 +11,7 @@ #include "Function.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** * Function for reading a 2D axisymmetric solution from file and mapping it to a @@ -44,7 +44,7 @@ class Axisymmetric2D3DSolutionFunction : public Function protected: /// Pointer to SolutionUserObject containing the solution of interest - const SolutionUserObject * _solution_object_ptr; + const SolutionUserObjectBase * _solution_object_ptr; /// Factor to scale the solution by (default = 1) const Real _scale_factor; diff --git a/framework/include/functions/SolutionFunction.h b/framework/include/functions/SolutionFunction.h index 05ab6802777e..e648f25567b4 100644 --- a/framework/include/functions/SolutionFunction.h +++ b/framework/include/functions/SolutionFunction.h @@ -11,7 +11,7 @@ #include "Function.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** Function for reading a solution from file * Creates a function that extracts values from a solution read from a file, @@ -55,7 +55,7 @@ class SolutionFunction : public Function protected: /// Pointer to SolutionUserObject containing the solution of interest - const SolutionUserObject * _solution_object_ptr; + const SolutionUserObjectBase * _solution_object_ptr; /// The local SolutionUserObject index for the variable extracted from the file unsigned int _solution_object_var_index; diff --git a/framework/include/ics/ScalarSolutionIC.h b/framework/include/ics/ScalarSolutionIC.h index 927317959c1a..1f79ddb49485 100644 --- a/framework/include/ics/ScalarSolutionIC.h +++ b/framework/include/ics/ScalarSolutionIC.h @@ -11,7 +11,7 @@ #include "ScalarInitialCondition.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** * Class for reading an initial condition from a solution user object @@ -25,7 +25,7 @@ class ScalarSolutionIC : public ScalarInitialCondition protected: /// SolutionUserObject containing the solution of interest - const SolutionUserObject & _solution_object; + const SolutionUserObjectBase & _solution_object; /// The variable name extracted from the SolutionUserObject const VariableName & _solution_object_var_name; diff --git a/framework/include/ics/SolutionIC.h b/framework/include/ics/SolutionIC.h index 7bd64b9ae4c7..bf5f7372834c 100644 --- a/framework/include/ics/SolutionIC.h +++ b/framework/include/ics/SolutionIC.h @@ -11,7 +11,7 @@ #include "InitialCondition.h" -class SolutionUserObject; +class SolutionUserObjectBase; /** * Class for reading an initial condition from a solution user object @@ -27,7 +27,7 @@ class SolutionIC : public InitialCondition protected: /// SolutionUserObject containing the solution of interest - const SolutionUserObject & _solution_object; + const SolutionUserObjectBase & _solution_object; /// The variable name extracted from the SolutionUserObject const VariableName & _solution_object_var_name; diff --git a/framework/include/userobjects/SolutionUserObject.h b/framework/include/userobjects/SolutionUserObject.h index caa608425e36..799e68dc312e 100644 --- a/framework/include/userobjects/SolutionUserObject.h +++ b/framework/include/userobjects/SolutionUserObject.h @@ -10,503 +10,25 @@ #pragma once // MOOSE includes -#include "GeneralUserObject.h" - -// Forward declarations -namespace libMesh -{ -class ExodusII_IO; -class EquationSystems; -class System; -class MeshFunction; -template -class NumericVector; -} +#include "SolutionUserObjectBase.h" +#include "FunctionParserUtils.h" /** * User object that reads an existing solution from an input file and * uses it in the current simulation. */ -class SolutionUserObject : public GeneralUserObject +class SolutionUserObject : public SolutionUserObjectBase, public FunctionParserUtils { public: static InputParameters validParams(); SolutionUserObject(const InputParameters & parameters); - virtual ~SolutionUserObject(); // empty dtor required for unique_ptr with forward declarations - - /** - * When reading ExodusII files, this will update the interpolation times - */ - virtual void timestepSetup() override; - - /** - * Returns the local index for a given variable name - * @param var_name The name of the variable for which the index is located - * @return The local index of the variable - */ - unsigned int getLocalVarIndex(const std::string & var_name) const; - - /** - * Returns a value at a specific location and variable checking for multiple values and weighting - * these values to - * obtain a single unique value (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - Real pointValueWrapper(Real t, - const Point & p, - const std::string & var_name, - const MooseEnum & weighting_type = weightingType(), - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns a value at a specific location and variable (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param local_var_index The local index of the variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - Real pointValue(Real t, - const Point & p, - const unsigned int local_var_index, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns a value at a specific location and variable (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - Real pointValue(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns a value at a specific location and variable for cases where the solution is - * multivalued at element faces - * Use pointValue for continuous shape functions or if you are sure your point is within an - * element! - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param local_var_index The local index of the variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - std::map - discontinuousPointValue(Real t, - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns a value at a specific location and variable for cases where the solution is - * multivalued at element faces - * Use pointValue for continuous shape functions or if you are sure your point is within an - * element! - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - std::map - discontinuousPointValue(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns the gradient at a specific location and variable checking for multiple values and - * weighting these values to - * obtain a single unique value (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - libMesh::RealGradient - pointValueGradientWrapper(Real t, - const Point & p, - const std::string & var_name, - const MooseEnum & weighting_type = weightingType(), - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns the gradient at a specific location and variable (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - libMesh::RealGradient - pointValueGradient(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns the gradient at a specific location and variable (see SolutionFunction) - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param local_var_index The local index of the variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - libMesh::RealGradient - pointValueGradient(Real t, - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns the gradient at a specific location and variable for cases where the gradient is - * multivalued (e.g. at element faces) - * Use pointValueGradient for continuous gradients or if you are sure your point is within an - * element! - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param var_name The variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - std::map discontinuousPointValueGradient( - Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids = nullptr) const; - - /** - * Returns the gradient at a specific location and variable for cases where the gradient is - * multivalued (e.g. at element faces) - * Use pointValueGradient for continuous gradients or if you are sure your point is within an - * element! - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param p The location at which to return a value - * @param local_var_index The local index of the variable to be evaluated - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @return The desired value for the given variable at a location - */ - std::map discontinuousPointValueGradient( - Real t, - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids = nullptr) const; - - /** - * Return a value directly from a Node - * @param node A pointer to the node at which a value is desired - * @param var_name The variable from which to extract a value - * @return The desired value for the given node and variable name - */ - Real directValue(const Node * node, const std::string & var_name) const; - - /** - * Return a value from the centroid of an element - * @param elem A pointer to the element at which a value is desired - * @param var_name The variable from which to extract a value - * @return The desired value for the given element and variable name - */ - Real directValue(const Elem * elem, const std::string & var_name) const; - - /** - * Returns a value of a global variable - * @param t The time at which to extract (not used, it is handled automatically when reading the - * data) - * @param var_name The variable from which to extract a value - * @return The desired value for the given variable - */ - Real scalarValue(Real t, const std::string & var_name) const; - - // Required pure virtual function (not used) - virtual void initialize() override; - - // Required pure virtual function (not used) - virtual void finalize() override; - - // Required pure virtual function (not used) - virtual void execute() override; - - /// Initialize the System and Mesh objects for the solution being read - virtual void initialSetup() override; - - const std::vector & variableNames() const; - - bool isVariableNodal(const std::string & var_name) const; - - static MooseEnum weightingType() - { - return MooseEnum("found_first=1 average=2 smallest_element_id=4 largest_element_id=8", - "found_first"); - } - - /** - * Return the spatial dimension of the mesh file - * @return The spatial dimension of the mesh file - */ - unsigned int getMeshFileDimension() const { return _mesh->spatial_dimension(); } - - /** - * Return the name of the mesh file this object read the solution from - */ - const std::string getMeshFileName() const { return _mesh_file; } /** - * Get the map from block name to block ID. Only works for ExodusII files. - * - * @return Map from block name to block ID + * Get the time at which to sample the solution */ - const std::map & getBlockNamesToIds() const - { - return _block_name_to_id; - } - - /** - * Get the map from block id to block name. Only works for ExodusII files. - * - * @return Map from block ID to block name - */ - const std::map & getBlockIdsToNames() const - { - return _block_id_to_name; - } - - /** - * Get the type of file that was read - * @return Returns a MooseEnum that specifies the type of file read - */ - MooseEnum getSolutionFileType() const; - -protected: - /** - * Method for reading XDA mesh and equation systems file(s) - * This method is called by the constructor when 'file_type = xda' is set - * in the input file. - */ - void readXda(); - - /** - * Method for reading an ExodusII file, which is called when - * 'file_type = exodusII is set in the input file. - */ - void readExodusII(); - - /** - * Method for extracting value of solution based on the DOF, - * this is called by the public overloaded function that accept - * a node or element pointer. - * @param dof_index The DOF of the solution desired - * @return The solution at the given DOF - */ - virtual Real directValue(dof_id_type dof_index) const; - - /** - * Updates the times for interpolating ExodusII data - * @param time The new time - */ - void updateExodusTimeInterpolation(Real time); - - /** - * Updates the time indices to interpolate between for ExodusII data - * @param time The new time - */ - bool updateExodusBracketingTimeIndices(Real time); - - /** - * A wrapper method for calling the various MeshFunctions used for reading the data - * @param p The location at which data is desired - * @param local_var_index The local index of the variable to extract data from - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) - */ - Real evalMeshFunction(const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids = nullptr) const; - - /** - * A wrapper method for calling the various MeshFunctions that calls the mesh function - * functionality for evaluating discontinuous shape functions near a face (where it's multivalued) - * @param p The location at which data is desired - * @param local_var_index The local index of the variable to extract data from - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) - */ - std::map - evalMultiValuedMeshFunction(const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids = nullptr) const; - - /** - * A wrapper method interfacing with the libMesh mesh function for evaluating the gradient - * @param p The location at which data is desired - * @param local_var_index The local index of the variable to extract data from - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) - */ - libMesh::RealGradient - evalMeshFunctionGradient(const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids = nullptr) const; - - /** - * A wrapper method interfacing with the libMesh mesh function that calls the gradient - * functionality for evaluating potentially discontinuous gradients at element's faces (where it's - * multivalued) - * @param p The location at which data is desired - * @param local_var_index The local index of the variable to extract data from - * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere - * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) - */ - std::map evalMultiValuedMeshFunctionGradient( - const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids = nullptr) const; - - /** - * Read block ID map from the ExodusII file - */ - void readBlockIdMapFromExodusII(); - - /// File type to read (0 = xda; 1 = ExodusII) - MooseEnum _file_type; - - /// The XDA or ExodusII file that is being read - std::string _mesh_file; - - /// The XDA file that contians the EquationSystems data (xda only) - std::string _es_file; - - /// The system name to extract from the XDA file (xda only) - std::string _system_name; - - /// A list of variables to extract from the read system - std::vector _system_variables; - - /// Stores the local index need by MeshFunction - std::map _local_variable_index; - - /// Stores names of nodal variables - std::vector _nodal_variables; - - /// Stores names of elemental variables - std::vector _elemental_variables; - - /// Stores names of scalar variables - std::vector _scalar_variables; - - /// Current ExodusII time index - int _exodus_time_index; - - /// Flag for triggering interpolation of ExodusII data - bool _interpolate_times; - - /// Pointer the libMesh::mesh object - std::unique_ptr _mesh; - - /// Pointer to the libMesh::EquationSystems object - std::unique_ptr _es; - - /// Pointer libMesh::System class storing the read solution - libMesh::System * _system; - - /// Pointer the libMesh::MeshFunction object that the read data is stored - std::unique_ptr _mesh_function; - - /// Pointer to the libMesh::ExodusII used to read the files - std::unique_ptr _exodusII_io; - - /// Pointer to the serial solution vector - std::unique_ptr> _serialized_solution; - - /// Pointer to second libMesh::EquationSystems object, used for interpolation - std::unique_ptr _es2; - - /// Pointer to a second libMesh::System object, used for interpolation - libMesh::System * _system2; - - /// Pointer to second libMesh::MeshFuntion, used for interpolation - std::unique_ptr _mesh_function2; - - /// Pointer to second serial solution, used for interpolation - std::unique_ptr> _serialized_solution2; - - /// Interpolation time - Real _interpolation_time; - - /// Interpolation weight factor - Real _interpolation_factor; - - /// The times available in the ExodusII file - const std::vector * _exodus_times; - - /// Time index 1, used for interpolation - int _exodus_index1; - - /// Time index 2, used for interpolation - int _exodus_index2; - - /// Scale parameter - std::vector _scale; - - /// scale_multiplier parameter - std::vector _scale_multiplier; - - /// Translation - std::vector _translation; - - /// vector about which to rotate - RealVectorValue _rotation0_vector; - - /// angle (in degrees) which to rotate through about vector _rotation0_vector - Real _rotation0_angle; - - /// Rotation matrix that performs the "_rotation0_angle about rotation0_vector" - RealTensorValue _r0; - - /// vector about which to rotate - RealVectorValue _rotation1_vector; - - /// angle (in degrees) which to rotate through about vector _rotation1_vector - Real _rotation1_angle; - - /// Rotation matrix that performs the "_rotation1_angle about rotation1_vector" - RealTensorValue _r1; - - /// transformations (rotations, translation, scales) are performed in this order - MultiMooseEnum _transformation_order; - - /// True if initial_setup has executed - bool _initialized; - - /// Map from block ID to block names. Read from the ExodusII file - std::map _block_name_to_id; - - /// Map from block names to block IDs. Read from the ExodusII file - std::map _block_id_to_name; + virtual Real solutionSampleTime() override; -private: - static Threads::spin_mutex _solution_user_object_mutex; + /// function parser object for transforming the solution sample time + SymFunctionPtr _time_transformation; }; diff --git a/framework/include/userobjects/SolutionUserObjectBase.h b/framework/include/userobjects/SolutionUserObjectBase.h new file mode 100644 index 000000000000..797a918f5982 --- /dev/null +++ b/framework/include/userobjects/SolutionUserObjectBase.h @@ -0,0 +1,514 @@ +//* 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 "GeneralUserObject.h" + +// Forward declarations +namespace libMesh +{ +class ExodusII_IO; +class EquationSystems; +class System; +class MeshFunction; +template +class NumericVector; +} + +/** + * User object that reads an existing solution from an input file and + * uses it in the current simulation. + */ +class SolutionUserObjectBase : public GeneralUserObject +{ +public: + static InputParameters validParams(); + + SolutionUserObjectBase(const InputParameters & parameters); + + /** + * Get the time at which to sample the solution + */ + virtual Real solutionSampleTime() = 0; + + /** + * When reading ExodusII files, this will update the interpolation times + */ + virtual void timestepSetup() override; + + /** + * Returns the local index for a given variable name + * @param var_name The name of the variable for which the index is located + * @return The local index of the variable + */ + unsigned int getLocalVarIndex(const std::string & var_name) const; + + /** + * Returns a value at a specific location and variable checking for multiple values and weighting + * these values to + * obtain a single unique value (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + Real pointValueWrapper(Real t, + const Point & p, + const std::string & var_name, + const MooseEnum & weighting_type = weightingType(), + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns a value at a specific location and variable (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param local_var_index The local index of the variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + Real pointValue(Real t, + const Point & p, + const unsigned int local_var_index, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns a value at a specific location and variable (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + Real pointValue(Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns a value at a specific location and variable for cases where the solution is + * multivalued at element faces + * Use pointValue for continuous shape functions or if you are sure your point is within an + * element! + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param local_var_index The local index of the variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + std::map + discontinuousPointValue(Real t, + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns a value at a specific location and variable for cases where the solution is + * multivalued at element faces + * Use pointValue for continuous shape functions or if you are sure your point is within an + * element! + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + std::map + discontinuousPointValue(Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns the gradient at a specific location and variable checking for multiple values and + * weighting these values to + * obtain a single unique value (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + libMesh::RealGradient + pointValueGradientWrapper(Real t, + const Point & p, + const std::string & var_name, + const MooseEnum & weighting_type = weightingType(), + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns the gradient at a specific location and variable (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + libMesh::RealGradient + pointValueGradient(Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns the gradient at a specific location and variable (see SolutionFunction) + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param local_var_index The local index of the variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + libMesh::RealGradient + pointValueGradient(Real t, + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns the gradient at a specific location and variable for cases where the gradient is + * multivalued (e.g. at element faces) + * Use pointValueGradient for continuous gradients or if you are sure your point is within an + * element! + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param var_name The variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + std::map discontinuousPointValueGradient( + Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids = nullptr) const; + + /** + * Returns the gradient at a specific location and variable for cases where the gradient is + * multivalued (e.g. at element faces) + * Use pointValueGradient for continuous gradients or if you are sure your point is within an + * element! + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param p The location at which to return a value + * @param local_var_index The local index of the variable to be evaluated + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @return The desired value for the given variable at a location + */ + std::map discontinuousPointValueGradient( + Real t, + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids = nullptr) const; + + /** + * Return a value directly from a Node + * @param node A pointer to the node at which a value is desired + * @param var_name The variable from which to extract a value + * @return The desired value for the given node and variable name + */ + Real directValue(const Node * node, const std::string & var_name) const; + + /** + * Return a value from the centroid of an element + * @param elem A pointer to the element at which a value is desired + * @param var_name The variable from which to extract a value + * @return The desired value for the given element and variable name + */ + Real directValue(const Elem * elem, const std::string & var_name) const; + + /** + * Returns a value of a global variable + * @param t The time at which to extract (not used, it is handled automatically when reading the + * data) + * @param var_name The variable from which to extract a value + * @return The desired value for the given variable + */ + Real scalarValue(Real t, const std::string & var_name) const; + + // Required pure virtual function (not used) + virtual void initialize() override; + + // Required pure virtual function (not used) + virtual void finalize() override; + + // Required pure virtual function (not used) + virtual void execute() override; + + /// Initialize the System and Mesh objects for the solution being read + virtual void initialSetup() override; + + const std::vector & variableNames() const; + + bool isVariableNodal(const std::string & var_name) const; + + static MooseEnum weightingType() + { + return MooseEnum("found_first=1 average=2 smallest_element_id=4 largest_element_id=8", + "found_first"); + } + + /** + * Return the spatial dimension of the mesh file + * @return The spatial dimension of the mesh file + */ + unsigned int getMeshFileDimension() const { return _mesh->spatial_dimension(); } + + /** + * Return the name of the mesh file this object read the solution from + */ + const std::string getMeshFileName() const { return _mesh_file; } + + /** + * Get the map from block name to block ID. Only works for ExodusII files. + * + * @return Map from block name to block ID + */ + const std::map & getBlockNamesToIds() const + { + return _block_name_to_id; + } + + /** + * Get the map from block id to block name. Only works for ExodusII files. + * + * @return Map from block ID to block name + */ + const std::map & getBlockIdsToNames() const + { + return _block_id_to_name; + } + + /** + * Get the type of file that was read + * @return Returns a MooseEnum that specifies the type of file read + */ + MooseEnum getSolutionFileType() const; + +protected: + /** + * Method for reading XDA mesh and equation systems file(s) + * This method is called by the constructor when 'file_type = xda' is set + * in the input file. + */ + void readXda(); + + /** + * Method for reading an ExodusII file, which is called when + * 'file_type = exodusII is set in the input file. + */ + void readExodusII(); + + /** + * Method for extracting value of solution based on the DOF, + * this is called by the public overloaded function that accept + * a node or element pointer. + * @param dof_index The DOF of the solution desired + * @return The solution at the given DOF + */ + virtual Real directValue(dof_id_type dof_index) const; + + /** + * Updates the times for interpolating ExodusII data + */ + void updateExodusTimeInterpolation(); + + /** + * Updates the time indices to interpolate between for ExodusII data + */ + bool updateExodusBracketingTimeIndices(); + + /** + * A wrapper method for calling the various MeshFunctions used for reading the data + * @param p The location at which data is desired + * @param local_var_index The local index of the variable to extract data from + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) + */ + Real evalMeshFunction(const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids = nullptr) const; + + /** + * A wrapper method for calling the various MeshFunctions that calls the mesh function + * functionality for evaluating discontinuous shape functions near a face (where it's multivalued) + * @param p The location at which data is desired + * @param local_var_index The local index of the variable to extract data from + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) + */ + std::map + evalMultiValuedMeshFunction(const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids = nullptr) const; + + /** + * A wrapper method interfacing with the libMesh mesh function for evaluating the gradient + * @param p The location at which data is desired + * @param local_var_index The local index of the variable to extract data from + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) + */ + libMesh::RealGradient + evalMeshFunctionGradient(const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids = nullptr) const; + + /** + * A wrapper method interfacing with the libMesh mesh function that calls the gradient + * functionality for evaluating potentially discontinuous gradients at element's faces (where it's + * multivalued) + * @param p The location at which data is desired + * @param local_var_index The local index of the variable to extract data from + * @param subdomain_ids Subdomains IDs where to look for the value, if nullptr look everywhere + * @param func_num The MeshFunction index to use (1 = _mesh_function; 2 = _mesh_function2) + */ + std::map evalMultiValuedMeshFunctionGradient( + const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids = nullptr) const; + + /** + * Read block ID map from the ExodusII file + */ + void readBlockIdMapFromExodusII(); + + /// File type to read (0 = xda; 1 = ExodusII) + MooseEnum _file_type; + + /// The XDA or ExodusII file that is being read + std::string _mesh_file; + + /// The XDA file that contians the EquationSystems data (xda only) + std::string _es_file; + + /// The system name to extract from the XDA file (xda only) + std::string _system_name; + + /// A list of variables to extract from the read system + std::vector _system_variables; + + /// Stores the local index need by MeshFunction + std::map _local_variable_index; + + /// Stores names of nodal variables + std::vector _nodal_variables; + + /// Stores names of elemental variables + std::vector _elemental_variables; + + /// Stores names of scalar variables + std::vector _scalar_variables; + + /// Current ExodusII time index + int _exodus_time_index; + + /// Flag for triggering interpolation of ExodusII data + bool _interpolate_times; + + /// Pointer the libMesh::mesh object + std::unique_ptr _mesh; + + /// Pointer to the libMesh::EquationSystems object + std::unique_ptr _es; + + /// Pointer libMesh::System class storing the read solution + libMesh::System * _system; + + /// Pointer the libMesh::MeshFunction object that the read data is stored + std::unique_ptr _mesh_function; + + /// Pointer to the libMesh::ExodusII used to read the files + std::unique_ptr _exodusII_io; + + /// Pointer to the serial solution vector + std::unique_ptr> _serialized_solution; + + /// Pointer to second libMesh::EquationSystems object, used for interpolation + std::unique_ptr _es2; + + /// Pointer to a second libMesh::System object, used for interpolation + libMesh::System * _system2; + + /// Pointer to second libMesh::MeshFuntion, used for interpolation + std::unique_ptr _mesh_function2; + + /// Pointer to second serial solution, used for interpolation + std::unique_ptr> _serialized_solution2; + + /// Time in the current simulation at which the solution interpolation was last updated + Real _interpolation_time; + + /// Interpolation weight factor + Real _interpolation_factor; + + /// The times available in the ExodusII file + const std::vector * _exodus_times; + + /// Time index 1, used for interpolation + int _exodus_index1; + + /// Time index 2, used for interpolation + int _exodus_index2; + + /// Scale parameter + std::vector _scale; + + /// scale_multiplier parameter + std::vector _scale_multiplier; + + /// Translation + std::vector _translation; + + /// vector about which to rotate + RealVectorValue _rotation0_vector; + + /// angle (in degrees) which to rotate through about vector _rotation0_vector + Real _rotation0_angle; + + /// Rotation matrix that performs the "_rotation0_angle about rotation0_vector" + RealTensorValue _r0; + + /// vector about which to rotate + RealVectorValue _rotation1_vector; + + /// angle (in degrees) which to rotate through about vector _rotation1_vector + Real _rotation1_angle; + + /// Rotation matrix that performs the "_rotation1_angle about rotation1_vector" + RealTensorValue _r1; + + /// transformations (rotations, translation, scales) are performed in this order + MultiMooseEnum _transformation_order; + + /// True if initial_setup has executed + bool _initialized; + + /// Map from block ID to block names. Read from the ExodusII file + std::map _block_name_to_id; + + /// Map from block names to block IDs. Read from the ExodusII file + std::map _block_id_to_name; + +private: + static Threads::spin_mutex _solution_user_object_mutex; +}; diff --git a/framework/src/auxkernels/SolutionAux.C b/framework/src/auxkernels/SolutionAux.C index b3330958e354..f2ac0cc1fd24 100644 --- a/framework/src/auxkernels/SolutionAux.C +++ b/framework/src/auxkernels/SolutionAux.C @@ -9,7 +9,7 @@ #include "MooseError.h" #include "SolutionAux.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" registerMooseObject("MooseApp", SolutionAux); @@ -39,7 +39,7 @@ SolutionAux::validParams() SolutionAux::SolutionAux(const InputParameters & parameters) : AuxKernel(parameters), - _solution_object(getUserObject("solution")), + _solution_object(getUserObject("solution")), _direct(getParam("direct")), _scale_factor(getParam("scale_factor")), _add_factor(getParam("add_factor")) diff --git a/framework/src/auxkernels/SolutionScalarAux.C b/framework/src/auxkernels/SolutionScalarAux.C index 2bb3acc175c1..24fa66066828 100644 --- a/framework/src/auxkernels/SolutionScalarAux.C +++ b/framework/src/auxkernels/SolutionScalarAux.C @@ -9,7 +9,7 @@ #include "MooseError.h" #include "SolutionScalarAux.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" registerMooseObject("MooseApp", SolutionScalarAux); @@ -35,7 +35,7 @@ SolutionScalarAux::validParams() SolutionScalarAux::SolutionScalarAux(const InputParameters & parameters) : AuxScalarKernel(parameters), - _solution_object(getUserObject("solution")), + _solution_object(getUserObject("solution")), _scale_factor(getParam("scale_factor")), _add_factor(getParam("add_factor")) { diff --git a/framework/src/functions/Axisymmetric2D3DSolutionFunction.C b/framework/src/functions/Axisymmetric2D3DSolutionFunction.C index e88636b85cc0..e3a42e5f325c 100644 --- a/framework/src/functions/Axisymmetric2D3DSolutionFunction.C +++ b/framework/src/functions/Axisymmetric2D3DSolutionFunction.C @@ -9,7 +9,7 @@ #include "MooseError.h" #include "Axisymmetric2D3DSolutionFunction.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" registerMooseObject("MooseApp", Axisymmetric2D3DSolutionFunction); @@ -107,7 +107,7 @@ Axisymmetric2D3DSolutionFunction::initialSetup() // Get a pointer to the SolutionUserObject. A pointer is used because the UserObject is not // available during the // construction of the function - _solution_object_ptr = &getUserObject("solution"); + _solution_object_ptr = &getUserObject("solution"); // If 'from_variable' is not specified, get the value from the SolutionUserObject if (_var_names.size() == 0) diff --git a/framework/src/functions/SolutionFunction.C b/framework/src/functions/SolutionFunction.C index e5dd7937468a..f67b14333b8e 100644 --- a/framework/src/functions/SolutionFunction.C +++ b/framework/src/functions/SolutionFunction.C @@ -9,7 +9,7 @@ #include "MooseError.h" #include "SolutionFunction.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" #include "MooseMesh.h" registerMooseObject("MooseApp", SolutionFunction); @@ -57,7 +57,7 @@ SolutionFunction::initialSetup() // Get a pointer to the SolutionUserObject. A pointer is used because the UserObject is not // available during the // construction of the function - _solution_object_ptr = &getUserObject("solution"); + _solution_object_ptr = &getUserObject("solution"); std::string var_name; diff --git a/framework/src/ics/ScalarSolutionIC.C b/framework/src/ics/ScalarSolutionIC.C index 09323ca8a206..5e3243ce53f8 100644 --- a/framework/src/ics/ScalarSolutionIC.C +++ b/framework/src/ics/ScalarSolutionIC.C @@ -8,7 +8,7 @@ //* https://www.gnu.org/licenses/lgpl-2.1.html #include "ScalarSolutionIC.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" #include "MooseMesh.h" registerMooseObject("MooseApp", ScalarSolutionIC); @@ -34,7 +34,7 @@ ScalarSolutionIC::validParams() ScalarSolutionIC::ScalarSolutionIC(const InputParameters & parameters) : ScalarInitialCondition(parameters), - _solution_object(getUserObject("solution_uo")), + _solution_object(getUserObject("solution_uo")), _solution_object_var_name(getParam("from_variable")) { } diff --git a/framework/src/ics/SolutionIC.C b/framework/src/ics/SolutionIC.C index ab6d86352e44..90448860a672 100644 --- a/framework/src/ics/SolutionIC.C +++ b/framework/src/ics/SolutionIC.C @@ -8,7 +8,7 @@ //* https://www.gnu.org/licenses/lgpl-2.1.html #include "SolutionIC.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" #include "MooseMesh.h" #include "SystemBase.h" @@ -31,7 +31,7 @@ SolutionIC::validParams() SolutionIC::SolutionIC(const InputParameters & parameters) : InitialCondition(parameters), - _solution_object(getUserObject("solution_uo")), + _solution_object(getUserObject("solution_uo")), _solution_object_var_name(getParam("from_variable")) { } diff --git a/framework/src/postprocessors/ParsedPostprocessor.C b/framework/src/postprocessors/ParsedPostprocessor.C index 21b795e22d56..5a221eddd749 100644 --- a/framework/src/postprocessors/ParsedPostprocessor.C +++ b/framework/src/postprocessors/ParsedPostprocessor.C @@ -31,7 +31,7 @@ ParsedPostprocessor::validParams() {}, "Vector of values for the constants in constant_names (can be an FParser expression)"); params.addParam( - "use_t", false, "Make time (t) variables available in the function expression."); + "use_t", false, "Make time (t) variable available in the function expression."); params.addClassDescription("Computes a parsed expression with post-processors"); return params; diff --git a/framework/src/userobjects/SolutionUserObject.C b/framework/src/userobjects/SolutionUserObject.C index 2d8d9d47f203..a2bb3345fdd8 100644 --- a/framework/src/userobjects/SolutionUserObject.C +++ b/framework/src/userobjects/SolutionUserObject.C @@ -16,6 +16,7 @@ #include "MooseUtils.h" #include "MooseVariableFE.h" #include "RotationMatrix.h" +#include "Function.h" // libMesh includes #include "libmesh/equation_systems.h" @@ -35,1224 +36,38 @@ InputParameters SolutionUserObject::validParams() { // Get the input parameters from the parent class - InputParameters params = GeneralUserObject::validParams(); + InputParameters params = SolutionUserObjectBase::validParams(); + params += FunctionParserUtils::validParams(); - // Add required parameters - params.addRequiredParam( - "mesh", "The name of the mesh file (must be xda/xdr or exodusII file)."); - params.addParam>( - "system_variables", - std::vector(), - "The name of the nodal and elemental variables from the file you want to use for values"); + params.addCustomTypeParam( + "time_transformation", + "t", + "FunctionExpression", + "Expression to transform from current simulation time to time at " + "which to sample the solution."); - // When using XDA/XDR files the following must be defined - params.addParam( - "es", - "", - "The name of the file holding the equation system info in xda/xdr format (xda/xdr only)."); - params.addParam( - "system", - "nl0", - "The name of the system to pull values out of (xda/xdr only). The default name for the " - "nonlinear system is 'nl0', auxiliary system is 'aux0'"); - - // When using ExodusII a specific time is extracted - params.addParam("timestep", - "Index of the single timestep used or \"LATEST\" for " - "the last timestep (exodusII only). If not supplied, " - "time interpolation will occur."); - - // Add ability to perform coordinate transformation: scale, factor - params.addParam>( - "scale", std::vector(LIBMESH_DIM, 1), "Scale factor for points in the simulation"); - params.addParam>("scale_multiplier", - std::vector(LIBMESH_DIM, 1), - "Scale multiplying factor for points in the simulation"); - params.addParam>("translation", - std::vector(LIBMESH_DIM, 0), - "Translation factors for x,y,z coordinates of the simulation"); - params.addParam("rotation0_vector", - RealVectorValue(0, 0, 1), - "Vector about which to rotate points of the simulation."); - params.addParam( - "rotation0_angle", - 0.0, - "Anticlockwise rotation angle (in degrees) to use for rotation about rotation0_vector."); - params.addParam("rotation1_vector", - RealVectorValue(0, 0, 1), - "Vector about which to rotate points of the simulation."); - params.addParam( - "rotation1_angle", - 0.0, - "Anticlockwise rotation angle (in degrees) to use for rotation about rotation1_vector."); - - // following lines build the default_transformation_order - MultiMooseEnum default_transformation_order( - "rotation0 translation scale rotation1 scale_multiplier", "translation scale"); - params.addParam( - "transformation_order", - default_transformation_order, - "The order to perform the operations in. Define R0 to be the rotation matrix encoded by " - "rotation0_vector and rotation0_angle. Similarly for R1. Denote the scale by s, the " - "scale_multiplier by m, and the translation by t. Then, given a point x in the simulation, " - "if transformation_order = 'rotation0 scale_multiplier translation scale rotation1' then " - "form p = R1*(R0*x*m - t)/s. Then the values provided by the SolutionUserObject at point x " - "in the simulation are the variable values at point p in the mesh."); - params.addParamNamesToGroup("scale scale_multiplier translation rotation0_vector rotation0_angle " - "rotation1_angle transformation_order", - "Coordinate system transformation"); - params.addClassDescription("Reads a variable from a mesh in one simulation to another"); return params; } -// Static mutex definition -Threads::spin_mutex SolutionUserObject::_solution_user_object_mutex; - SolutionUserObject::SolutionUserObject(const InputParameters & parameters) - : GeneralUserObject(parameters), - _file_type(MooseEnum("xda=0 exodusII=1 xdr=2")), - _mesh_file(getParam("mesh")), - _es_file(getParam("es")), - _system_name(getParam("system")), - _system_variables(getParam>("system_variables")), - _exodus_time_index(-1), - _interpolate_times(false), - _system(nullptr), - _system2(nullptr), - _interpolation_time(0.0), - _interpolation_factor(0.0), - _exodus_times(nullptr), - _exodus_index1(-1), - _exodus_index2(-1), - _scale(getParam>("scale")), - _scale_multiplier(getParam>("scale_multiplier")), - _translation(getParam>("translation")), - _rotation0_vector(getParam("rotation0_vector")), - _rotation0_angle(getParam("rotation0_angle")), - _r0(RealTensorValue()), - _rotation1_vector(getParam("rotation1_vector")), - _rotation1_angle(getParam("rotation1_angle")), - _r1(RealTensorValue()), - _transformation_order(getParam("transformation_order")), - _initialized(false) -{ - // form rotation matrices with the specified angles - Real halfPi = std::acos(0.0); - Real a; - Real b; - - a = std::cos(halfPi * -_rotation0_angle / 90); - b = std::sin(halfPi * -_rotation0_angle / 90); - // the following is an anticlockwise rotation about z - RealTensorValue rot0_z(a, -b, 0, b, a, 0, 0, 0, 1); - // form the rotation matrix that will take rotation0_vector to the z axis - RealTensorValue vec0_to_z = RotationMatrix::rotVecToZ(_rotation0_vector); - // _r0 is then: rotate points so vec0 lies along z; then rotate about angle0; then rotate points - // back - _r0 = vec0_to_z.transpose() * (rot0_z * vec0_to_z); - - a = std::cos(halfPi * -_rotation1_angle / 90); - b = std::sin(halfPi * -_rotation1_angle / 90); - // the following is an anticlockwise rotation about z - RealTensorValue rot1_z(a, -b, 0, b, a, 0, 0, 0, 1); - // form the rotation matrix that will take rotation1_vector to the z axis - RealTensorValue vec1_to_z = RotationMatrix::rotVecToZ(_rotation1_vector); - // _r1 is then: rotate points so vec1 lies along z; then rotate about angle1; then rotate points - // back - _r1 = vec1_to_z.transpose() * (rot1_z * vec1_to_z); - - if (isParamValid("timestep") && getParam("timestep") == "-1") - mooseError("A \"timestep\" of -1 is no longer supported for interpolation. Instead simply " - "remove this parameter altogether for interpolation"); -} - -SolutionUserObject::~SolutionUserObject() {} - -void -SolutionUserObject::readXda() -{ - if (!isParamSetByUser("es")) - paramError("es", "Equation system file (.xda or .xdr) should have been specified"); - - // Check that the required files exist - MooseUtils::checkFileReadable(_es_file); - MooseUtils::checkFileReadable(_mesh_file); - - // Read the libmesh::mesh from the xda file - _mesh->read(_mesh_file); - - // Create the libmesh::EquationSystems - _es = std::make_unique(*_mesh); - - // Use new read syntax (binary) - if (_file_type == "xdr") - _es->read(_es_file, - libMesh::DECODE, - EquationSystems::READ_HEADER | EquationSystems::READ_DATA | - EquationSystems::READ_ADDITIONAL_DATA); - - // Use new read syntax - else if (_file_type == "xda") - _es->read(_es_file, - libMesh::READ, - EquationSystems::READ_HEADER | EquationSystems::READ_DATA | - EquationSystems::READ_ADDITIONAL_DATA); - - // This should never occur, just in case produce an error - else - mooseError("Failed to determine proper read method for XDA/XDR equation system file: ", - _es_file); - - // Update and store the EquationSystems name locally - _es->update(); - _system = &_es->get_system(_system_name); -} - -void -SolutionUserObject::readExodusII() -{ - // Define a default system name - if (_system_name == "") - _system_name = "SolutionUserObjectSystem"; - - // Read the Exodus file - _exodusII_io = std::make_unique(*_mesh); - _exodusII_io->read(_mesh_file); - readBlockIdMapFromExodusII(); - _exodus_times = &_exodusII_io->get_time_steps(); - - if (isParamValid("timestep")) - { - std::string s_timestep = getParam("timestep"); - int n_steps = _exodusII_io->get_num_time_steps(); - if (s_timestep == "LATEST") - _exodus_time_index = n_steps; - else - { - std::istringstream ss(s_timestep); - if (!((ss >> _exodus_time_index) && ss.eof()) || _exodus_time_index > n_steps) - mooseError("Invalid value passed as \"timestep\". Expected \"LATEST\" or a valid integer " - "less than ", - n_steps, - ", received ", - s_timestep); - } - } - else - // Interpolate between times rather than using values from a set timestep - _interpolate_times = true; - - // Check that the number of time steps is valid - int num_exo_times = _exodus_times->size(); - if (num_exo_times == 0) - mooseError("In SolutionUserObject, exodus file contains no timesteps."); - - // Account for parallel mesh - if (dynamic_cast(_mesh.get())) - { - _mesh->allow_renumbering(true); - _mesh->prepare_for_use(); - } - else - { - _mesh->allow_renumbering(false); - _mesh->prepare_for_use(); - } - - // Create EquationSystems object for solution - _es = std::make_unique(*_mesh); - _es->add_system(_system_name); - _system = &_es->get_system(_system_name); - - // Get the variable name lists as set; these need to be sets to perform set_intersection - const std::vector & all_nodal(_exodusII_io->get_nodal_var_names()); - const std::vector & all_elemental(_exodusII_io->get_elem_var_names()); - const std::vector & all_scalar(_exodusII_io->get_global_var_names()); - - // Build nodal/elemental variable lists, limit to variables listed in 'system_variables', if - // provided - // This function could be called more than once, so clear the member variables so we don't keep - // adding to the vectors - _nodal_variables.clear(); - _elemental_variables.clear(); - _scalar_variables.clear(); - if (!_system_variables.empty()) - { - for (const auto & var_name : _system_variables) - { - if (std::find(all_nodal.begin(), all_nodal.end(), var_name) != all_nodal.end()) - _nodal_variables.push_back(var_name); - if (std::find(all_elemental.begin(), all_elemental.end(), var_name) != all_elemental.end()) - _elemental_variables.push_back(var_name); - if (std::find(all_scalar.begin(), all_scalar.end(), var_name) != all_scalar.end()) - // Check if the scalar matches any field variables, and ignore the var if it does. This - // means its a Postprocessor. - if (std::find(begin(_nodal_variables), end(_nodal_variables), var_name) == - _nodal_variables.end() && - std::find(begin(_elemental_variables), end(_elemental_variables), var_name) == - _elemental_variables.end()) - _scalar_variables.push_back(var_name); - } - } - else - { - _nodal_variables = all_nodal; - _elemental_variables = all_elemental; - - for (auto var_name : all_scalar) - // Check if the scalar matches any field variables, and ignore the var if it does. This means - // its a Postprocessor. - if (std::find(begin(_nodal_variables), end(_nodal_variables), var_name) == - _nodal_variables.end() && - std::find(begin(_elemental_variables), end(_elemental_variables), var_name) == - _elemental_variables.end()) - _scalar_variables.push_back(var_name); - } - - // Add the variables to the system - for (const auto & var_name : _nodal_variables) - _system->add_variable(var_name, FIRST); - - for (const auto & var_name : _elemental_variables) - _system->add_variable(var_name, CONSTANT, MONOMIAL); - - for (const auto & var_name : _scalar_variables) - _system->add_variable(var_name, FIRST, SCALAR); - - // Initialize the equations systems - _es->init(); - - // Interpolate times - if (_interpolate_times) - { - // Create a second equation system - _es2 = std::make_unique(*_mesh); - _es2->add_system(_system_name); - _system2 = &_es2->get_system(_system_name); - - // Add the variables to the system - for (const auto & var_name : _nodal_variables) - _system2->add_variable(var_name, FIRST); - - for (const auto & var_name : _elemental_variables) - _system2->add_variable(var_name, CONSTANT, MONOMIAL); - - for (const auto & var_name : _scalar_variables) - _system2->add_variable(var_name, FIRST, SCALAR); - - // Initialize - _es2->init(); - - // Update the times for interpolation (initially start at 0) - updateExodusBracketingTimeIndices(0.0); - - // Copy the solutions from the first system - for (const auto & var_name : _nodal_variables) - { - _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1); - _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1); - } - - for (const auto & var_name : _elemental_variables) - { - _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1); - _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1); - } - - if (_scalar_variables.size() > 0) - { - _exodusII_io->copy_scalar_solution( - *_system, _scalar_variables, _scalar_variables, _exodus_index1 + 1); - _exodusII_io->copy_scalar_solution( - *_system2, _scalar_variables, _scalar_variables, _exodus_index2 + 1); - } - - // Update the systems - _system->update(); - _es->update(); - _system2->update(); - _es2->update(); - } - - // Non-interpolated times - else - { - if (_exodus_time_index > num_exo_times) - mooseError("In SolutionUserObject, timestep = ", - _exodus_time_index, - ", but there are only ", - num_exo_times, - " time steps."); - - // Copy the values from the ExodusII file - for (const auto & var_name : _nodal_variables) - _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_time_index); - - for (const auto & var_name : _elemental_variables) - _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_time_index); - - if (!_scalar_variables.empty()) - _exodusII_io->copy_scalar_solution( - *_system, _scalar_variables, _scalar_variables, _exodus_time_index); - - // Update the equations systems - _system->update(); - _es->update(); - } -} - -Real -SolutionUserObject::directValue(const Node * node, const std::string & var_name) const + : SolutionUserObjectBase(parameters), FunctionParserUtils(parameters) { - // Get the libmesh variable and system numbers - unsigned int var_num = _system->variable_number(var_name); - unsigned int sys_num = _system->number(); + // setup parsed expression for the time transformation + _time_transformation = std::make_shared(); + setParserFeatureFlags(_time_transformation); - // Get the node id and associated dof - dof_id_type node_id = node->id(); - const Node & sys_node = _system->get_mesh().node_ref(node_id); - mooseAssert(sys_node.n_dofs(sys_num, var_num) > 0, - "Variable " << var_name << " has no DoFs on node " << sys_node.id()); - dof_id_type dof_id = sys_node.dof_number(sys_num, var_num, 0); + // parse function + const auto & expression = getParam("time_transformation"); + if (_time_transformation->Parse(expression, "t") >= 0) + mooseError("Invalid parsed function\n", expression, "\n", _time_transformation->ErrorMsg()); - // Return the desired value for the dof - return directValue(dof_id); + // the only parameter is time + _func_params.resize(1); } Real -SolutionUserObject::directValue(const Elem * elem, const std::string & var_name) const -{ - // Get the libmesh variable and system numbers - unsigned int var_num = _system->variable_number(var_name); - unsigned int sys_num = _system->number(); - - // Get the element id and associated dof - dof_id_type elem_id = elem->id(); - const Elem & sys_elem = _system->get_mesh().elem_ref(elem_id); - mooseAssert(sys_elem.n_dofs(sys_num, var_num) > 0, - "Variable " << var_name << " has no DoFs on element " << sys_elem.id()); - dof_id_type dof_id = sys_elem.dof_number(sys_num, var_num, 0); - - // Return the desired value - return directValue(dof_id); -} - -void -SolutionUserObject::initialize() -{ -} - -void -SolutionUserObject::finalize() -{ -} - -void -SolutionUserObject::timestepSetup() -{ - // Update time interpolation for ExodusII solution - if (_file_type == 1 && _interpolate_times) - updateExodusTimeInterpolation(_t); -} - -void -SolutionUserObject::execute() -{ -} - -void -SolutionUserObject::initialSetup() -{ - - // Make sure this only happens once - if (_initialized) - return; - - // Create a libmesh::Mesh object for storing the loaded data. - // Several aspects of SolutionUserObject won't work with a DistributedMesh: - // .) ExodusII_IO::copy_nodal_solution() doesn't work in parallel. - // .) We don't know if directValue will be used, which may request - // a value on a Node we don't have. - // We force the Mesh used here to be a ReplicatedMesh. - _mesh = std::make_unique(_communicator); - - // ExodusII mesh file supplied - if (MooseUtils::hasExtension(_mesh_file, "e", /*strip_exodus_ext =*/true) || - MooseUtils::hasExtension(_mesh_file, "exo", true)) - { - _file_type = "exodusII"; - readExodusII(); - } - - // XDA mesh file supplied - else if (MooseUtils::hasExtension(_mesh_file, "xda")) - { - _file_type = "xda"; - readXda(); - } - - // XDR mesh file supplied - else if (MooseUtils::hasExtension(_mesh_file, "xdr")) - { - _file_type = "xdr"; - readXda(); - } - - // Produce an error for an unknown file type - else - mooseError("In SolutionUserObject, invalid file type: only .xda, .xdr, .exo and .e supported"); - - // Initialize the serial solution vector - _serialized_solution = NumericVector::build(_communicator); - _serialized_solution->init(_system->n_dofs(), false, libMesh::SERIAL); - - // Pull down a full copy of this vector on every processor so we can get values in parallel - _system->solution->localize(*_serialized_solution); - - // Vector of variable numbers to apply the MeshFunction to - std::vector var_nums; - - // If no variables were given, use all of them - if (_system_variables.empty()) - { - _system->get_all_variable_numbers(var_nums); - for (const auto & var_num : var_nums) - _system_variables.push_back(_system->variable_name(var_num)); - } - - // Otherwise, gather the numbers for the variables given - else - { - for (const auto & var_name : _system_variables) - var_nums.push_back(_system->variable_number(var_name)); - } - - // Create the MeshFunction for working with the solution data - _mesh_function = std::make_unique( - *_es, *_serialized_solution, _system->get_dof_map(), var_nums); - _mesh_function->init(); - - // Tell the MeshFunctions that we might be querying them outside the - // mesh, so we can handle any errors at the MOOSE rather than at the - // libMesh level. - DenseVector default_values; - _mesh_function->enable_out_of_mesh_mode(default_values); - - // Build second MeshFunction for interpolation - if (_interpolate_times) - { - // Need to pull down a full copy of this vector on every processor so we can get values in - // parallel - _serialized_solution2 = NumericVector::build(_communicator); - _serialized_solution2->init(_system2->n_dofs(), false, libMesh::SERIAL); - _system2->solution->localize(*_serialized_solution2); - - // Create the MeshFunction for the second copy of the data - _mesh_function2 = std::make_unique( - *_es2, *_serialized_solution2, _system2->get_dof_map(), var_nums); - _mesh_function2->init(); - _mesh_function2->enable_out_of_mesh_mode(default_values); - } - - // Populate the MeshFunction variable index - for (unsigned int i = 0; i < _system_variables.size(); ++i) - { - std::string name = _system_variables[i]; - _local_variable_index[name] = i; - } - - // Set initialization flag - _initialized = true; -} - -MooseEnum -SolutionUserObject::getSolutionFileType() const -{ - return _file_type; -} - -void -SolutionUserObject::updateExodusTimeInterpolation(Real time) -{ - if (time != _interpolation_time) - { - if (updateExodusBracketingTimeIndices(time)) - { - - for (const auto & var_name : _nodal_variables) - _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1); - for (const auto & var_name : _elemental_variables) - _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1); - if (_scalar_variables.size() > 0) - _exodusII_io->copy_scalar_solution( - *_system, _scalar_variables, _scalar_variables, _exodus_index1 + 1); - - _system->update(); - _es->update(); - _system->solution->localize(*_serialized_solution); - - for (const auto & var_name : _nodal_variables) - _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1); - for (const auto & var_name : _elemental_variables) - _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1); - if (_scalar_variables.size() > 0) - _exodusII_io->copy_scalar_solution( - *_system2, _scalar_variables, _scalar_variables, _exodus_index2 + 1); - - _system2->update(); - _es2->update(); - _system2->solution->localize(*_serialized_solution2); - } - _interpolation_time = time; - } -} - -bool -SolutionUserObject::updateExodusBracketingTimeIndices(Real time) -{ - if (_file_type != 1) - mooseError( - "In SolutionUserObject, getTimeInterpolationData only applicable for exodusII file type"); - - int old_index1 = _exodus_index1; - int old_index2 = _exodus_index2; - - int num_exo_times = _exodus_times->size(); - - if (time < (*_exodus_times)[0]) - { - _exodus_index1 = 0; - _exodus_index2 = 0; - _interpolation_factor = 0.0; - } - else - { - for (int i = 0; i < num_exo_times - 1; ++i) - { - if (time <= (*_exodus_times)[i + 1]) - { - _exodus_index1 = i; - _exodus_index2 = i + 1; - _interpolation_factor = - (time - (*_exodus_times)[i]) / ((*_exodus_times)[i + 1] - (*_exodus_times)[i]); - break; - } - else if (i == num_exo_times - 2) - { - _exodus_index1 = num_exo_times - 1; - _exodus_index2 = num_exo_times - 1; - _interpolation_factor = 1.0; - break; - } - } - } - - bool indices_modified(false); - - if (_exodus_index1 != old_index1 || _exodus_index2 != old_index2) - indices_modified = true; - - return indices_modified; -} - -unsigned int -SolutionUserObject::getLocalVarIndex(const std::string & var_name) const -{ - // Extract the variable index for the MeshFunction(s) - std::map::const_iterator it = _local_variable_index.find(var_name); - if (it == _local_variable_index.end()) - mooseError("Value requested for nonexistent variable '", - var_name, - "' in the '", - name(), - "' SolutionUserObject.\nSystem selected: ", - _system_name, - "\nAvailable variables:\n", - ConsoleUtils::formatString(Moose::stringify(_system_variables), "")); - return it->second; -} - -Real -SolutionUserObject::pointValueWrapper(Real t, - const Point & p, - const std::string & var_name, - const MooseEnum & weighting_type, - const std::set * subdomain_ids) const -{ - // first check if the FE type is continuous because in that case the value is - // unique and we can take a short cut, the default weighting_type found_first also - // shortcuts out - auto family = - _fe_problem - .getVariable( - _tid, var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD) - .feType() - .family; - - if (weighting_type == 1 || - (family != L2_LAGRANGE && family != MONOMIAL && family != L2_HIERARCHIC)) - return pointValue(t, p, var_name, subdomain_ids); - - // the shape function is discontinuous so we need to compute a suitable unique value - std::map values = discontinuousPointValue(t, p, var_name); - switch (weighting_type) - { - case 2: - { - Real average = 0.0; - for (auto & v : values) - average += v.second; - return average / Real(values.size()); - } - case 4: - { - Real smallest_elem_id_value = std::numeric_limits::max(); - dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId(); - for (auto & v : values) - if (v.first->id() < smallest_elem_id) - { - smallest_elem_id = v.first->id(); - smallest_elem_id_value = v.second; - } - return smallest_elem_id_value; - } - case 8: - { - Real largest_elem_id_value = std::numeric_limits::lowest(); - dof_id_type largest_elem_id = 0; - for (auto & v : values) - if (v.first->id() > largest_elem_id) - { - largest_elem_id = v.first->id(); - largest_elem_id_value = v.second; - } - return largest_elem_id_value; - } - } - - mooseError( - "SolutionUserObject::pointValueWrapper reaches line that it should not be able to reach."); - return 0.0; -} - -Real -SolutionUserObject::pointValue(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids) const -{ - const unsigned int local_var_index = getLocalVarIndex(var_name); - return pointValue(t, p, local_var_index, subdomain_ids); -} - -Real -SolutionUserObject::pointValue(Real libmesh_dbg_var(t), - const Point & p, - const unsigned int local_var_index, - const std::set * subdomain_ids) const -{ - // Create copy of point - Point pt(p); - - // do the transformations - for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) - { - if (_transformation_order[trans_num] == "rotation0") - pt = _r0 * pt; - else if (_transformation_order[trans_num] == "translation") - for (const auto i : make_range(Moose::dim)) - pt(i) -= _translation[i]; - else if (_transformation_order[trans_num] == "scale") - for (const auto i : make_range(Moose::dim)) - pt(i) /= _scale[i]; - else if (_transformation_order[trans_num] == "scale_multiplier") - for (const auto i : make_range(Moose::dim)) - pt(i) *= _scale_multiplier[i]; - else if (_transformation_order[trans_num] == "rotation1") - pt = _r1 * pt; - } - - // Extract the value at the current point - Real val = evalMeshFunction(pt, local_var_index, 1, subdomain_ids); - - // Interpolate - if (_file_type == 1 && _interpolate_times) - { - mooseAssert(t == _interpolation_time, - "Time passed into value() must match time at last call to timestepSetup()"); - Real val2 = evalMeshFunction(pt, local_var_index, 2, subdomain_ids); - val = val + (val2 - val) * _interpolation_factor; - } - - return val; -} - -std::map -SolutionUserObject::discontinuousPointValue(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids) const -{ - const unsigned int local_var_index = getLocalVarIndex(var_name); - return discontinuousPointValue(t, p, local_var_index, subdomain_ids); -} - -std::map -SolutionUserObject::discontinuousPointValue(Real libmesh_dbg_var(t), - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids) const -{ - // do the transformations - for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) - { - if (_transformation_order[trans_num] == "rotation0") - pt = _r0 * pt; - else if (_transformation_order[trans_num] == "translation") - for (const auto i : make_range(Moose::dim)) - pt(i) -= _translation[i]; - else if (_transformation_order[trans_num] == "scale") - for (const auto i : make_range(Moose::dim)) - pt(i) /= _scale[i]; - else if (_transformation_order[trans_num] == "scale_multiplier") - for (const auto i : make_range(Moose::dim)) - pt(i) *= _scale_multiplier[i]; - else if (_transformation_order[trans_num] == "rotation1") - pt = _r1 * pt; - } - - // Extract the value at the current point - std::map map = - evalMultiValuedMeshFunction(pt, local_var_index, 1, subdomain_ids); - - // Interpolate - if (_file_type == 1 && _interpolate_times) - { - mooseAssert(t == _interpolation_time, - "Time passed into value() must match time at last call to timestepSetup()"); - std::map map2 = evalMultiValuedMeshFunction(pt, local_var_index, 2); - - if (map.size() != map2.size()) - mooseError("In SolutionUserObject::discontinuousPointValue map and map2 have different size"); - - // construct the interpolated map - for (auto & k : map) - { - if (map2.find(k.first) == map2.end()) - mooseError( - "In SolutionUserObject::discontinuousPointValue map and map2 have differing keys"); - Real val = k.second; - Real val2 = map2[k.first]; - map[k.first] = val + (val2 - val) * _interpolation_factor; - } - } - - return map; -} - -RealGradient -SolutionUserObject::pointValueGradientWrapper( - Real t, - const Point & p, - const std::string & var_name, - const MooseEnum & weighting_type, - const std::set * subdomain_ids) const -{ - // the default weighting_type found_first shortcuts out - if (weighting_type == 1) - return pointValueGradient(t, p, var_name, subdomain_ids); - - // the shape function is discontinuous so we need to compute a suitable unique value - std::map values = - discontinuousPointValueGradient(t, p, var_name, subdomain_ids); - switch (weighting_type) - { - case 2: - { - RealGradient average = RealGradient(0.0, 0.0, 0.0); - for (auto & v : values) - average += v.second; - return average / Real(values.size()); - } - case 4: - { - RealGradient smallest_elem_id_value; - dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId(); - for (auto & v : values) - if (v.first->id() < smallest_elem_id) - { - smallest_elem_id = v.first->id(); - smallest_elem_id_value = v.second; - } - return smallest_elem_id_value; - } - case 8: - { - RealGradient largest_elem_id_value; - dof_id_type largest_elem_id = 0; - for (auto & v : values) - if (v.first->id() > largest_elem_id) - { - largest_elem_id = v.first->id(); - largest_elem_id_value = v.second; - } - return largest_elem_id_value; - } - } - - mooseError("SolutionUserObject::pointValueGradientWrapper reaches line that it should not be " - "able to reach."); - return RealGradient(0.0, 0.0, 0.0); -} - -RealGradient -SolutionUserObject::pointValueGradient(Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids) const -{ - const unsigned int local_var_index = getLocalVarIndex(var_name); - return pointValueGradient(t, p, local_var_index, subdomain_ids); -} - -RealGradient -SolutionUserObject::pointValueGradient(Real libmesh_dbg_var(t), - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids) const -{ - // do the transformations - for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) - { - if (_transformation_order[trans_num] == "rotation0") - pt = _r0 * pt; - else if (_transformation_order[trans_num] == "translation") - for (const auto i : make_range(Moose::dim)) - pt(i) -= _translation[i]; - else if (_transformation_order[trans_num] == "scale") - for (const auto i : make_range(Moose::dim)) - pt(i) /= _scale[i]; - else if (_transformation_order[trans_num] == "scale_multiplier") - for (const auto i : make_range(Moose::dim)) - pt(i) *= _scale_multiplier[i]; - else if (_transformation_order[trans_num] == "rotation1") - pt = _r1 * pt; - } - - // Extract the value at the current point - RealGradient val = evalMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); - - // Interpolate - if (_file_type == 1 && _interpolate_times) - { - mooseAssert(t == _interpolation_time, - "Time passed into value() must match time at last call to timestepSetup()"); - RealGradient val2 = evalMeshFunctionGradient(pt, local_var_index, 2, subdomain_ids); - val = val + (val2 - val) * _interpolation_factor; - } - - return val; -} - -std::map -SolutionUserObject::discontinuousPointValueGradient( - Real t, - const Point & p, - const std::string & var_name, - const std::set * subdomain_ids) const -{ - const unsigned int local_var_index = getLocalVarIndex(var_name); - return discontinuousPointValueGradient(t, p, local_var_index, subdomain_ids); -} - -std::map -SolutionUserObject::discontinuousPointValueGradient( - Real libmesh_dbg_var(t), - Point pt, - const unsigned int local_var_index, - const std::set * subdomain_ids) const -{ - // do the transformations - for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) - { - if (_transformation_order[trans_num] == "rotation0") - pt = _r0 * pt; - else if (_transformation_order[trans_num] == "translation") - for (const auto i : make_range(Moose::dim)) - pt(i) -= _translation[i]; - else if (_transformation_order[trans_num] == "scale") - for (const auto i : make_range(Moose::dim)) - pt(i) /= _scale[i]; - else if (_transformation_order[trans_num] == "scale_multiplier") - for (const auto i : make_range(Moose::dim)) - pt(i) *= _scale_multiplier[i]; - else if (_transformation_order[trans_num] == "rotation1") - pt = _r1 * pt; - } - - // Extract the value at the current point - std::map map = - evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); - - // Interpolate - if (_file_type == 1 && _interpolate_times) - { - mooseAssert(t == _interpolation_time, - "Time passed into value() must match time at last call to timestepSetup()"); - std::map map2 = - evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); - - if (map.size() != map2.size()) - mooseError("In SolutionUserObject::discontinuousPointValue map and map2 have different size"); - - // construct the interpolated map - for (auto & k : map) - { - if (map2.find(k.first) == map2.end()) - mooseError( - "In SolutionUserObject::discontinuousPointValue map and map2 have differing keys"); - RealGradient val = k.second; - RealGradient val2 = map2[k.first]; - map[k.first] = val + (val2 - val) * _interpolation_factor; - } - } - - return map; -} - -Real -SolutionUserObject::directValue(dof_id_type dof_index) const -{ - Real val = (*_serialized_solution)(dof_index); - if (_file_type == 1 && _interpolate_times) - { - Real val2 = (*_serialized_solution2)(dof_index); - val = val + (val2 - val) * _interpolation_factor; - } - return val; -} - -Real -SolutionUserObject::evalMeshFunction(const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids) const -{ - // Storage for mesh function output - DenseVector output; - - // Extract a value from the _mesh_function - { - Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); - if (func_num == 1) - (*_mesh_function)(p, 0.0, output, subdomain_ids); - - // Extract a value from _mesh_function2 - else if (func_num == 2) - (*_mesh_function2)(p, 0.0, output, subdomain_ids); - - else - mooseError("The func_num must be 1 or 2"); - } - - // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated - // outside the domain - if (output.size() == 0) - { - std::ostringstream oss; - p.print(oss); - mooseError("Failed to access the data for variable '", - _system_variables[local_var_index], - "' at point ", - oss.str(), - " in the '", - name(), - "' SolutionUserObject"); - } - return output(local_var_index); -} - -std::map -SolutionUserObject::evalMultiValuedMeshFunction( - const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids) const -{ - // Storage for mesh function output - std::map> temporary_output; - - // Extract a value from the _mesh_function - { - Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); - if (func_num == 1) - _mesh_function->discontinuous_value(p, 0.0, temporary_output, subdomain_ids); - - // Extract a value from _mesh_function2 - else if (func_num == 2) - _mesh_function2->discontinuous_value(p, 0.0, temporary_output, subdomain_ids); - - else - mooseError("The func_num must be 1 or 2"); - } - - // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated - // outside the domain - if (temporary_output.size() == 0) - { - std::ostringstream oss; - p.print(oss); - mooseError("Failed to access the data for variable '", - _system_variables[local_var_index], - "' at point ", - oss.str(), - " in the '", - name(), - "' SolutionUserObject"); - } - - // Fill the actual map that is returned - std::map output; - for (auto & k : temporary_output) - { - mooseAssert(k.second.size() > local_var_index, - "In SolutionUserObject::evalMultiValuedMeshFunction variable with local_var_index " - << local_var_index << " does not exist"); - output[k.first] = k.second(local_var_index); - } - - return output; -} - -RealGradient -SolutionUserObject::evalMeshFunctionGradient( - const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids) const -{ - // Storage for mesh function output - std::vector output; - - // Extract a value from the _mesh_function - { - Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); - if (func_num == 1) - _mesh_function->gradient(p, 0.0, output, subdomain_ids); - - // Extract a value from _mesh_function2 - else if (func_num == 2) - _mesh_function2->gradient(p, 0.0, output, subdomain_ids); - - else - mooseError("The func_num must be 1 or 2"); - } - - // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated - // outside the domain - if (output.size() == 0) - { - std::ostringstream oss; - p.print(oss); - mooseError("Failed to access the data for variable '", - _system_variables[local_var_index], - "' at point ", - oss.str(), - " in the '", - name(), - "' SolutionUserObject"); - } - return output[local_var_index]; -} - -std::map -SolutionUserObject::evalMultiValuedMeshFunctionGradient( - const Point & p, - const unsigned int local_var_index, - unsigned int func_num, - const std::set * subdomain_ids) const -{ - // Storage for mesh function output - std::map> temporary_output; - - // Extract a value from the _mesh_function - { - Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); - if (func_num == 1) - _mesh_function->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids); - - // Extract a value from _mesh_function2 - else if (func_num == 2) - _mesh_function2->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids); - - else - mooseError("The func_num must be 1 or 2"); - } - - // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated - // outside the domain - if (temporary_output.size() == 0) - { - std::ostringstream oss; - p.print(oss); - mooseError("Failed to access the data for variable '", - _system_variables[local_var_index], - "' at point ", - oss.str(), - " in the '", - name(), - "' SolutionUserObject"); - } - - // Fill the actual map that is returned - std::map output; - for (auto & k : temporary_output) - { - mooseAssert(k.second.size() > local_var_index, - "In SolutionUserObject::evalMultiValuedMeshFunction variable with local_var_index " - << local_var_index << " does not exist"); - output[k.first] = k.second[local_var_index]; - } - - return output; -} - -const std::vector & -SolutionUserObject::variableNames() const -{ - return _system_variables; -} - -bool -SolutionUserObject::isVariableNodal(const std::string & var_name) const -{ - return std::find(_nodal_variables.begin(), _nodal_variables.end(), var_name) != - _nodal_variables.end(); -} - -Real -SolutionUserObject::scalarValue(Real /*t*/, const std::string & var_name) const -{ - unsigned int var_num = _system->variable_number(var_name); - const DofMap & dof_map = _system->get_dof_map(); - std::vector dofs; - dof_map.SCALAR_dof_indices(dofs, var_num); - // We can handle only FIRST order scalar variables - return directValue(dofs[0]); -} - -void -SolutionUserObject::readBlockIdMapFromExodusII() +SolutionUserObject::solutionSampleTime() { -#ifdef LIBMESH_HAVE_EXODUS_API - libMesh::ExodusII_IO_Helper & exio_helper = _exodusII_io->get_exio_helper(); - const auto & id_to_block = exio_helper.id_to_block_names; - _block_name_to_id.clear(); - _block_id_to_name.clear(); - for (const auto & it : id_to_block) - { - _block_name_to_id[it.second] = it.first; - _block_id_to_name[it.first] = it.second; - } -#endif + _func_params[0] = _t; + return evaluate(_time_transformation); } diff --git a/framework/src/userobjects/SolutionUserObjectBase.C b/framework/src/userobjects/SolutionUserObjectBase.C new file mode 100644 index 000000000000..722afa09b599 --- /dev/null +++ b/framework/src/userobjects/SolutionUserObjectBase.C @@ -0,0 +1,1265 @@ +//* 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 "SolutionUserObjectBase.h" + +// MOOSE includes +#include "ConsoleUtils.h" +#include "MooseError.h" +#include "MooseMesh.h" +#include "MooseUtils.h" +#include "MooseVariableFE.h" +#include "RotationMatrix.h" +#include "Function.h" + +// libMesh includes +#include "libmesh/equation_systems.h" +#include "libmesh/mesh_function.h" +#include "libmesh/numeric_vector.h" +#include "libmesh/nonlinear_implicit_system.h" +#include "libmesh/transient_system.h" +#include "libmesh/parallel_mesh.h" +#include "libmesh/serial_mesh.h" +#include "libmesh/exodusII_io.h" +#include "libmesh/exodusII_io_helper.h" +#include "libmesh/enum_xdr_mode.h" + +InputParameters +SolutionUserObjectBase::validParams() +{ + // Get the input parameters from the parent class + InputParameters params = GeneralUserObject::validParams(); + + // Add required parameters + params.addRequiredParam( + "mesh", "The name of the mesh file (must be xda/xdr or exodusII file)."); + params.addParam>( + "system_variables", + std::vector(), + "The name of the nodal and elemental variables from the file you want to use for values"); + + // When using XDA/XDR files the following must be defined + params.addParam( + "es", + "", + "The name of the file holding the equation system info in xda/xdr format (xda/xdr only)."); + params.addParam( + "system", + "nl0", + "The name of the system to pull values out of (xda/xdr only). The default name for the " + "nonlinear system is 'nl0', auxiliary system is 'aux0'"); + + // When using ExodusII a specific time is extracted + params.addParam("timestep", + "Index of the single timestep used or \"LATEST\" for " + "the last timestep (exodusII only). If not supplied, " + "time interpolation will occur."); + + // Add ability to perform coordinate transformation: scale, factor + params.addParam>( + "scale", std::vector(LIBMESH_DIM, 1), "Scale factor for points in the simulation"); + params.addParam>("scale_multiplier", + std::vector(LIBMESH_DIM, 1), + "Scale multiplying factor for points in the simulation"); + params.addParam>("translation", + std::vector(LIBMESH_DIM, 0), + "Translation factors for x,y,z coordinates of the simulation"); + params.addParam("rotation0_vector", + RealVectorValue(0, 0, 1), + "Vector about which to rotate points of the simulation."); + params.addParam( + "rotation0_angle", + 0.0, + "Anticlockwise rotation angle (in degrees) to use for rotation about rotation0_vector."); + params.addParam("rotation1_vector", + RealVectorValue(0, 0, 1), + "Vector about which to rotate points of the simulation."); + params.addParam( + "rotation1_angle", + 0.0, + "Anticlockwise rotation angle (in degrees) to use for rotation about rotation1_vector."); + + // following lines build the default_transformation_order + MultiMooseEnum default_transformation_order( + "rotation0 translation scale rotation1 scale_multiplier", "translation scale"); + params.addParam( + "transformation_order", + default_transformation_order, + "The order to perform the operations in. Define R0 to be the rotation matrix encoded by " + "rotation0_vector and rotation0_angle. Similarly for R1. Denote the scale by s, the " + "scale_multiplier by m, and the translation by t. Then, given a point x in the simulation, " + "if transformation_order = 'rotation0 scale_multiplier translation scale rotation1' then " + "form p = R1*(R0*x*m - t)/s. Then the values provided by the SolutionUserObjectBase at " + "point x " + "in the simulation are the variable values at point p in the mesh."); + params.addParamNamesToGroup("scale scale_multiplier translation rotation0_vector rotation0_angle " + "rotation1_angle transformation_order", + "Coordinate system transformation"); + params.addClassDescription("Reads a variable from a mesh in one simulation to another"); + return params; +} + +// Static mutex definition +Threads::spin_mutex SolutionUserObjectBase::_solution_user_object_mutex; + +SolutionUserObjectBase::SolutionUserObjectBase(const InputParameters & parameters) + : GeneralUserObject(parameters), + _file_type(MooseEnum("xda=0 exodusII=1 xdr=2")), + _mesh_file(getParam("mesh")), + _es_file(getParam("es")), + _system_name(getParam("system")), + _system_variables(getParam>("system_variables")), + _exodus_time_index(-1), + _interpolate_times(false), + _system(nullptr), + _system2(nullptr), + _interpolation_time(0.0), + _interpolation_factor(0.0), + _exodus_times(nullptr), + _exodus_index1(-1), + _exodus_index2(-1), + _scale(getParam>("scale")), + _scale_multiplier(getParam>("scale_multiplier")), + _translation(getParam>("translation")), + _rotation0_vector(getParam("rotation0_vector")), + _rotation0_angle(getParam("rotation0_angle")), + _r0(RealTensorValue()), + _rotation1_vector(getParam("rotation1_vector")), + _rotation1_angle(getParam("rotation1_angle")), + _r1(RealTensorValue()), + _transformation_order(getParam("transformation_order")), + _initialized(false) +{ + // form rotation matrices with the specified angles + Real halfPi = libMesh::pi / 2.0; + Real a; + Real b; + + a = std::cos(halfPi * -_rotation0_angle / 90.0); + b = std::sin(halfPi * -_rotation0_angle / 90.0); + // the following is an anticlockwise rotation about z + RealTensorValue rot0_z(a, -b, 0, b, a, 0, 0, 0, 1); + // form the rotation matrix that will take rotation0_vector to the z axis + RealTensorValue vec0_to_z = RotationMatrix::rotVecToZ(_rotation0_vector); + // _r0 is then: rotate points so vec0 lies along z; then rotate about angle0; then rotate points + // back + _r0 = vec0_to_z.transpose() * (rot0_z * vec0_to_z); + + a = std::cos(halfPi * -_rotation1_angle / 90.0); + b = std::sin(halfPi * -_rotation1_angle / 90.0); + // the following is an anticlockwise rotation about z + RealTensorValue rot1_z(a, -b, 0, b, a, 0, 0, 0, 1); + // form the rotation matrix that will take rotation1_vector to the z axis + RealTensorValue vec1_to_z = RotationMatrix::rotVecToZ(_rotation1_vector); + // _r1 is then: rotate points so vec1 lies along z; then rotate about angle1; then rotate points + // back + _r1 = vec1_to_z.transpose() * (rot1_z * vec1_to_z); + + if (isParamValid("timestep") && getParam("timestep") == "-1") + mooseError("A \"timestep\" of -1 is no longer supported for interpolation. Instead simply " + "remove this parameter altogether for interpolation"); +} + +void +SolutionUserObjectBase::readXda() +{ + if (!isParamSetByUser("es")) + paramError("es", "Equation system file (.xda or .xdr) should have been specified"); + + // Check that the required files exist + MooseUtils::checkFileReadable(_es_file); + MooseUtils::checkFileReadable(_mesh_file); + + // Read the libmesh::mesh from the xda file + _mesh->read(_mesh_file); + + // Create the libmesh::EquationSystems + _es = std::make_unique(*_mesh); + + // Use new read syntax (binary) + if (_file_type == "xdr") + _es->read(_es_file, + libMesh::DECODE, + EquationSystems::READ_HEADER | EquationSystems::READ_DATA | + EquationSystems::READ_ADDITIONAL_DATA); + + // Use new read syntax + else if (_file_type == "xda") + _es->read(_es_file, + libMesh::READ, + EquationSystems::READ_HEADER | EquationSystems::READ_DATA | + EquationSystems::READ_ADDITIONAL_DATA); + + // This should never occur, just in case produce an error + else + mooseError("Failed to determine proper read method for XDA/XDR equation system file: ", + _es_file); + + // Update and store the EquationSystems name locally + _es->update(); + _system = &_es->get_system(_system_name); +} + +void +SolutionUserObjectBase::readExodusII() +{ + // Define a default system name + if (_system_name == "") + _system_name = "SolutionUserObjectSystem"; + + // Read the Exodus file + _exodusII_io = std::make_unique(*_mesh); + _exodusII_io->read(_mesh_file); + readBlockIdMapFromExodusII(); + _exodus_times = &_exodusII_io->get_time_steps(); + + if (isParamValid("timestep")) + { + std::string s_timestep = getParam("timestep"); + int n_steps = _exodusII_io->get_num_time_steps(); + if (s_timestep == "LATEST") + _exodus_time_index = n_steps; + else + { + std::istringstream ss(s_timestep); + if (!((ss >> _exodus_time_index) && ss.eof()) || _exodus_time_index > n_steps) + mooseError("Invalid value passed as \"timestep\". Expected \"LATEST\" or a valid integer " + "less than ", + n_steps, + ", received ", + s_timestep); + } + } + else + // Interpolate between times rather than using values from a set timestep + _interpolate_times = true; + + // Check that the number of time steps is valid + int num_exo_times = _exodus_times->size(); + if (num_exo_times == 0) + mooseError("In SolutionUserObjectBase, exodus file contains no timesteps."); + + // Account for parallel mesh + if (dynamic_cast(_mesh.get())) + { + _mesh->allow_renumbering(true); + _mesh->prepare_for_use(); + } + else + { + _mesh->allow_renumbering(false); + _mesh->prepare_for_use(); + } + + // Create EquationSystems object for solution + _es = std::make_unique(*_mesh); + _es->add_system(_system_name); + _system = &_es->get_system(_system_name); + + // Get the variable name lists as set; these need to be sets to perform set_intersection + const std::vector & all_nodal(_exodusII_io->get_nodal_var_names()); + const std::vector & all_elemental(_exodusII_io->get_elem_var_names()); + const std::vector & all_scalar(_exodusII_io->get_global_var_names()); + + // Build nodal/elemental variable lists, limit to variables listed in 'system_variables', if + // provided + // This function could be called more than once, so clear the member variables so we don't keep + // adding to the vectors + _nodal_variables.clear(); + _elemental_variables.clear(); + _scalar_variables.clear(); + if (!_system_variables.empty()) + { + for (const auto & var_name : _system_variables) + { + if (std::find(all_nodal.begin(), all_nodal.end(), var_name) != all_nodal.end()) + _nodal_variables.push_back(var_name); + if (std::find(all_elemental.begin(), all_elemental.end(), var_name) != all_elemental.end()) + _elemental_variables.push_back(var_name); + if (std::find(all_scalar.begin(), all_scalar.end(), var_name) != all_scalar.end()) + // Check if the scalar matches any field variables, and ignore the var if it does. This + // means its a Postprocessor. + if (std::find(begin(_nodal_variables), end(_nodal_variables), var_name) == + _nodal_variables.end() && + std::find(begin(_elemental_variables), end(_elemental_variables), var_name) == + _elemental_variables.end()) + _scalar_variables.push_back(var_name); + } + } + else + { + _nodal_variables = all_nodal; + _elemental_variables = all_elemental; + + for (auto var_name : all_scalar) + // Check if the scalar matches any field variables, and ignore the var if it does. This means + // its a Postprocessor. + if (std::find(begin(_nodal_variables), end(_nodal_variables), var_name) == + _nodal_variables.end() && + std::find(begin(_elemental_variables), end(_elemental_variables), var_name) == + _elemental_variables.end()) + _scalar_variables.push_back(var_name); + } + + // Add the variables to the system + for (const auto & var_name : _nodal_variables) + _system->add_variable(var_name, FIRST); + + for (const auto & var_name : _elemental_variables) + _system->add_variable(var_name, CONSTANT, MONOMIAL); + + for (const auto & var_name : _scalar_variables) + _system->add_variable(var_name, FIRST, SCALAR); + + // Initialize the equations systems + _es->init(); + + // Interpolate times + if (_interpolate_times) + { + // Create a second equation system + _es2 = std::make_unique(*_mesh); + _es2->add_system(_system_name); + _system2 = &_es2->get_system(_system_name); + + // Add the variables to the system + for (const auto & var_name : _nodal_variables) + _system2->add_variable(var_name, FIRST); + + for (const auto & var_name : _elemental_variables) + _system2->add_variable(var_name, CONSTANT, MONOMIAL); + + for (const auto & var_name : _scalar_variables) + _system2->add_variable(var_name, FIRST, SCALAR); + + // Initialize + _es2->init(); + + // Update the times for interpolation (initially start at 0) + updateExodusBracketingTimeIndices(); + + // Copy the solutions from the first system + for (const auto & var_name : _nodal_variables) + { + _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1); + _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1); + } + + for (const auto & var_name : _elemental_variables) + { + _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1); + _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1); + } + + if (_scalar_variables.size() > 0) + { + _exodusII_io->copy_scalar_solution( + *_system, _scalar_variables, _scalar_variables, _exodus_index1 + 1); + _exodusII_io->copy_scalar_solution( + *_system2, _scalar_variables, _scalar_variables, _exodus_index2 + 1); + } + + // Update the systems + _system->update(); + _es->update(); + _system2->update(); + _es2->update(); + } + + // Non-interpolated times + else + { + if (_exodus_time_index > num_exo_times) + mooseError("In SolutionUserObjectBase, timestep = ", + _exodus_time_index, + ", but there are only ", + num_exo_times, + " time steps."); + + // Copy the values from the ExodusII file + for (const auto & var_name : _nodal_variables) + _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_time_index); + + for (const auto & var_name : _elemental_variables) + _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_time_index); + + if (!_scalar_variables.empty()) + _exodusII_io->copy_scalar_solution( + *_system, _scalar_variables, _scalar_variables, _exodus_time_index); + + // Update the equations systems + _system->update(); + _es->update(); + } +} + +Real +SolutionUserObjectBase::directValue(const Node * node, const std::string & var_name) const +{ + // Get the libmesh variable and system numbers + unsigned int var_num = _system->variable_number(var_name); + unsigned int sys_num = _system->number(); + + // Get the node id and associated dof + dof_id_type node_id = node->id(); + const Node & sys_node = _system->get_mesh().node_ref(node_id); + mooseAssert(sys_node.n_dofs(sys_num, var_num) > 0, + "Variable " << var_name << " has no DoFs on node " << sys_node.id()); + dof_id_type dof_id = sys_node.dof_number(sys_num, var_num, 0); + + // Return the desired value for the dof + return directValue(dof_id); +} + +Real +SolutionUserObjectBase::directValue(const Elem * elem, const std::string & var_name) const +{ + // Get the libmesh variable and system numbers + unsigned int var_num = _system->variable_number(var_name); + unsigned int sys_num = _system->number(); + + // Get the element id and associated dof + dof_id_type elem_id = elem->id(); + const Elem & sys_elem = _system->get_mesh().elem_ref(elem_id); + mooseAssert(sys_elem.n_dofs(sys_num, var_num) > 0, + "Variable " << var_name << " has no DoFs on element " << sys_elem.id()); + dof_id_type dof_id = sys_elem.dof_number(sys_num, var_num, 0); + + // Return the desired value + return directValue(dof_id); +} + +void +SolutionUserObjectBase::initialize() +{ +} + +void +SolutionUserObjectBase::finalize() +{ +} + +void +SolutionUserObjectBase::timestepSetup() +{ + // Update time interpolation for ExodusII solution + if (_file_type == 1 && _interpolate_times) + updateExodusTimeInterpolation(); +} + +void +SolutionUserObjectBase::execute() +{ +} + +void +SolutionUserObjectBase::initialSetup() +{ + + // Make sure this only happens once + if (_initialized) + return; + + // Create a libmesh::Mesh object for storing the loaded data. + // Several aspects of SolutionUserObjectBase won't work with a DistributedMesh: + // .) ExodusII_IO::copy_nodal_solution() doesn't work in parallel. + // .) We don't know if directValue will be used, which may request + // a value on a Node we don't have. + // We force the Mesh used here to be a ReplicatedMesh. + _mesh = std::make_unique(_communicator); + + // ExodusII mesh file supplied + if (MooseUtils::hasExtension(_mesh_file, "e", /*strip_exodus_ext =*/true) || + MooseUtils::hasExtension(_mesh_file, "exo", true)) + { + _file_type = "exodusII"; + readExodusII(); + } + + // XDA mesh file supplied + else if (MooseUtils::hasExtension(_mesh_file, "xda")) + { + _file_type = "xda"; + readXda(); + } + + // XDR mesh file supplied + else if (MooseUtils::hasExtension(_mesh_file, "xdr")) + { + _file_type = "xdr"; + readXda(); + } + + // Produce an error for an unknown file type + else + mooseError( + "In SolutionUserObjectBase, invalid file type: only .xda, .xdr, .exo and .e supported"); + + // Initialize the serial solution vector + _serialized_solution = NumericVector::build(_communicator); + _serialized_solution->init(_system->n_dofs(), false, libMesh::SERIAL); + + // Pull down a full copy of this vector on every processor so we can get values in parallel + _system->solution->localize(*_serialized_solution); + + // Vector of variable numbers to apply the MeshFunction to + std::vector var_nums; + + // If no variables were given, use all of them + if (_system_variables.empty()) + { + _system->get_all_variable_numbers(var_nums); + for (const auto & var_num : var_nums) + _system_variables.push_back(_system->variable_name(var_num)); + } + + // Otherwise, gather the numbers for the variables given + else + { + for (const auto & var_name : _system_variables) + var_nums.push_back(_system->variable_number(var_name)); + } + + // Create the MeshFunction for working with the solution data + _mesh_function = std::make_unique( + *_es, *_serialized_solution, _system->get_dof_map(), var_nums); + _mesh_function->init(); + + // Tell the MeshFunctions that we might be querying them outside the + // mesh, so we can handle any errors at the MOOSE rather than at the + // libMesh level. + DenseVector default_values; + _mesh_function->enable_out_of_mesh_mode(default_values); + + // Build second MeshFunction for interpolation + if (_interpolate_times) + { + // Need to pull down a full copy of this vector on every processor so we can get values in + // parallel + _serialized_solution2 = NumericVector::build(_communicator); + _serialized_solution2->init(_system2->n_dofs(), false, libMesh::SERIAL); + _system2->solution->localize(*_serialized_solution2); + + // Create the MeshFunction for the second copy of the data + _mesh_function2 = std::make_unique( + *_es2, *_serialized_solution2, _system2->get_dof_map(), var_nums); + _mesh_function2->init(); + _mesh_function2->enable_out_of_mesh_mode(default_values); + } + + // Populate the MeshFunction variable index + for (unsigned int i = 0; i < _system_variables.size(); ++i) + { + std::string name = _system_variables[i]; + _local_variable_index[name] = i; + } + + // Set initialization flag + _initialized = true; +} + +MooseEnum +SolutionUserObjectBase::getSolutionFileType() const +{ + return _file_type; +} + +void +SolutionUserObjectBase::updateExodusTimeInterpolation() +{ + if (_t != _interpolation_time) + { + if (updateExodusBracketingTimeIndices()) + { + + for (const auto & var_name : _nodal_variables) + _exodusII_io->copy_nodal_solution(*_system, var_name, var_name, _exodus_index1 + 1); + for (const auto & var_name : _elemental_variables) + _exodusII_io->copy_elemental_solution(*_system, var_name, var_name, _exodus_index1 + 1); + if (_scalar_variables.size() > 0) + _exodusII_io->copy_scalar_solution( + *_system, _scalar_variables, _scalar_variables, _exodus_index1 + 1); + + _system->update(); + _es->update(); + _system->solution->localize(*_serialized_solution); + + for (const auto & var_name : _nodal_variables) + _exodusII_io->copy_nodal_solution(*_system2, var_name, var_name, _exodus_index2 + 1); + for (const auto & var_name : _elemental_variables) + _exodusII_io->copy_elemental_solution(*_system2, var_name, var_name, _exodus_index2 + 1); + if (_scalar_variables.size() > 0) + _exodusII_io->copy_scalar_solution( + *_system2, _scalar_variables, _scalar_variables, _exodus_index2 + 1); + + _system2->update(); + _es2->update(); + _system2->solution->localize(*_serialized_solution2); + } + _interpolation_time = _t; + } +} + +bool +SolutionUserObjectBase::updateExodusBracketingTimeIndices() +{ + if (_file_type != 1) + mooseError("In SolutionUserObjectBase, getTimeInterpolationData only applicable for exodusII " + "file type"); + + int old_index1 = _exodus_index1; + int old_index2 = _exodus_index2; + + int num_exo_times = _exodus_times->size(); + + const auto solution_time = solutionSampleTime(); + + if (solution_time < (*_exodus_times)[0]) + { + _exodus_index1 = 0; + _exodus_index2 = 0; + _interpolation_factor = 0.0; + } + else + { + for (int i = 0; i < num_exo_times - 1; ++i) + { + if (solution_time <= (*_exodus_times)[i + 1]) + { + _exodus_index1 = i; + _exodus_index2 = i + 1; + _interpolation_factor = + (solution_time - (*_exodus_times)[i]) / ((*_exodus_times)[i + 1] - (*_exodus_times)[i]); + break; + } + else if (i == num_exo_times - 2) + { + _exodus_index1 = num_exo_times - 1; + _exodus_index2 = num_exo_times - 1; + _interpolation_factor = 1.0; + break; + } + } + } + + bool indices_modified(false); + + if (_exodus_index1 != old_index1 || _exodus_index2 != old_index2) + indices_modified = true; + + return indices_modified; +} + +unsigned int +SolutionUserObjectBase::getLocalVarIndex(const std::string & var_name) const +{ + // Extract the variable index for the MeshFunction(s) + std::map::const_iterator it = _local_variable_index.find(var_name); + if (it == _local_variable_index.end()) + mooseError("Value requested for nonexistent variable '", + var_name, + "' in the '", + name(), + "' SolutionUserObjectBase.\nSystem selected: ", + _system_name, + "\nAvailable variables:\n", + ConsoleUtils::formatString(Moose::stringify(_system_variables), "")); + return it->second; +} + +Real +SolutionUserObjectBase::pointValueWrapper(Real t, + const Point & p, + const std::string & var_name, + const MooseEnum & weighting_type, + const std::set * subdomain_ids) const +{ + // first check if the FE type is continuous because in that case the value is + // unique and we can take a short cut, the default weighting_type found_first also + // shortcuts out + auto family = + _fe_problem + .getVariable( + _tid, var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD) + .feType() + .family; + + if (weighting_type == 1 || + (family != L2_LAGRANGE && family != MONOMIAL && family != L2_HIERARCHIC)) + return pointValue(t, p, var_name, subdomain_ids); + + // the shape function is discontinuous so we need to compute a suitable unique value + std::map values = discontinuousPointValue(t, p, var_name); + switch (weighting_type) + { + case 2: + { + Real average = 0.0; + for (auto & v : values) + average += v.second; + return average / Real(values.size()); + } + case 4: + { + Real smallest_elem_id_value = std::numeric_limits::max(); + dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId(); + for (auto & v : values) + if (v.first->id() < smallest_elem_id) + { + smallest_elem_id = v.first->id(); + smallest_elem_id_value = v.second; + } + return smallest_elem_id_value; + } + case 8: + { + Real largest_elem_id_value = std::numeric_limits::lowest(); + dof_id_type largest_elem_id = 0; + for (auto & v : values) + if (v.first->id() > largest_elem_id) + { + largest_elem_id = v.first->id(); + largest_elem_id_value = v.second; + } + return largest_elem_id_value; + } + } + + mooseError("SolutionUserObjectBase::pointValueWrapper reaches line that it should not be able to " + "reach."); + return 0.0; +} + +Real +SolutionUserObjectBase::pointValue(Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids) const +{ + const unsigned int local_var_index = getLocalVarIndex(var_name); + return pointValue(t, p, local_var_index, subdomain_ids); +} + +Real +SolutionUserObjectBase::pointValue(Real libmesh_dbg_var(t), + const Point & p, + const unsigned int local_var_index, + const std::set * subdomain_ids) const +{ + // Create copy of point + Point pt(p); + + // do the transformations + for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) + { + if (_transformation_order[trans_num] == "rotation0") + pt = _r0 * pt; + else if (_transformation_order[trans_num] == "translation") + for (const auto i : make_range(Moose::dim)) + pt(i) -= _translation[i]; + else if (_transformation_order[trans_num] == "scale") + for (const auto i : make_range(Moose::dim)) + pt(i) /= _scale[i]; + else if (_transformation_order[trans_num] == "scale_multiplier") + for (const auto i : make_range(Moose::dim)) + pt(i) *= _scale_multiplier[i]; + else if (_transformation_order[trans_num] == "rotation1") + pt = _r1 * pt; + } + + // Extract the value at the current point + Real val = evalMeshFunction(pt, local_var_index, 1, subdomain_ids); + + // Interpolate + if (_file_type == 1 && _interpolate_times) + { + mooseAssert(t == _interpolation_time, + "Time passed into value() must match time at last call to timestepSetup()"); + Real val2 = evalMeshFunction(pt, local_var_index, 2, subdomain_ids); + val = val + (val2 - val) * _interpolation_factor; + } + + return val; +} + +std::map +SolutionUserObjectBase::discontinuousPointValue( + Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids) const +{ + const unsigned int local_var_index = getLocalVarIndex(var_name); + return discontinuousPointValue(t, p, local_var_index, subdomain_ids); +} + +std::map +SolutionUserObjectBase::discontinuousPointValue( + Real libmesh_dbg_var(t), + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids) const +{ + // do the transformations + for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) + { + if (_transformation_order[trans_num] == "rotation0") + pt = _r0 * pt; + else if (_transformation_order[trans_num] == "translation") + for (const auto i : make_range(Moose::dim)) + pt(i) -= _translation[i]; + else if (_transformation_order[trans_num] == "scale") + for (const auto i : make_range(Moose::dim)) + pt(i) /= _scale[i]; + else if (_transformation_order[trans_num] == "scale_multiplier") + for (const auto i : make_range(Moose::dim)) + pt(i) *= _scale_multiplier[i]; + else if (_transformation_order[trans_num] == "rotation1") + pt = _r1 * pt; + } + + // Extract the value at the current point + std::map map = + evalMultiValuedMeshFunction(pt, local_var_index, 1, subdomain_ids); + + // Interpolate + if (_file_type == 1 && _interpolate_times) + { + mooseAssert(t == _interpolation_time, + "Time passed into value() must match time at last call to timestepSetup()"); + std::map map2 = evalMultiValuedMeshFunction(pt, local_var_index, 2); + + if (map.size() != map2.size()) + mooseError( + "In SolutionUserObjectBase::discontinuousPointValue map and map2 have different size"); + + // construct the interpolated map + for (auto & k : map) + { + if (map2.find(k.first) == map2.end()) + mooseError( + "In SolutionUserObjectBase::discontinuousPointValue map and map2 have differing keys"); + Real val = k.second; + Real val2 = map2[k.first]; + map[k.first] = val + (val2 - val) * _interpolation_factor; + } + } + + return map; +} + +RealGradient +SolutionUserObjectBase::pointValueGradientWrapper( + Real t, + const Point & p, + const std::string & var_name, + const MooseEnum & weighting_type, + const std::set * subdomain_ids) const +{ + // the default weighting_type found_first shortcuts out + if (weighting_type == 1) + return pointValueGradient(t, p, var_name, subdomain_ids); + + // the shape function is discontinuous so we need to compute a suitable unique value + std::map values = + discontinuousPointValueGradient(t, p, var_name, subdomain_ids); + switch (weighting_type) + { + case 2: + { + RealGradient average = RealGradient(0.0, 0.0, 0.0); + for (auto & v : values) + average += v.second; + return average / Real(values.size()); + } + case 4: + { + RealGradient smallest_elem_id_value; + dof_id_type smallest_elem_id = _fe_problem.mesh().maxElemId(); + for (auto & v : values) + if (v.first->id() < smallest_elem_id) + { + smallest_elem_id = v.first->id(); + smallest_elem_id_value = v.second; + } + return smallest_elem_id_value; + } + case 8: + { + RealGradient largest_elem_id_value; + dof_id_type largest_elem_id = 0; + for (auto & v : values) + if (v.first->id() > largest_elem_id) + { + largest_elem_id = v.first->id(); + largest_elem_id_value = v.second; + } + return largest_elem_id_value; + } + } + + mooseError("SolutionUserObjectBase::pointValueGradientWrapper reaches line that it should not be " + "able to reach."); + return RealGradient(0.0, 0.0, 0.0); +} + +RealGradient +SolutionUserObjectBase::pointValueGradient(Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids) const +{ + const unsigned int local_var_index = getLocalVarIndex(var_name); + return pointValueGradient(t, p, local_var_index, subdomain_ids); +} + +RealGradient +SolutionUserObjectBase::pointValueGradient(Real libmesh_dbg_var(t), + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids) const +{ + // do the transformations + for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) + { + if (_transformation_order[trans_num] == "rotation0") + pt = _r0 * pt; + else if (_transformation_order[trans_num] == "translation") + for (const auto i : make_range(Moose::dim)) + pt(i) -= _translation[i]; + else if (_transformation_order[trans_num] == "scale") + for (const auto i : make_range(Moose::dim)) + pt(i) /= _scale[i]; + else if (_transformation_order[trans_num] == "scale_multiplier") + for (const auto i : make_range(Moose::dim)) + pt(i) *= _scale_multiplier[i]; + else if (_transformation_order[trans_num] == "rotation1") + pt = _r1 * pt; + } + + // Extract the value at the current point + RealGradient val = evalMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); + + // Interpolate + if (_file_type == 1 && _interpolate_times) + { + mooseAssert(t == _interpolation_time, + "Time passed into value() must match time at last call to timestepSetup()"); + RealGradient val2 = evalMeshFunctionGradient(pt, local_var_index, 2, subdomain_ids); + val = val + (val2 - val) * _interpolation_factor; + } + + return val; +} + +std::map +SolutionUserObjectBase::discontinuousPointValueGradient( + Real t, + const Point & p, + const std::string & var_name, + const std::set * subdomain_ids) const +{ + const unsigned int local_var_index = getLocalVarIndex(var_name); + return discontinuousPointValueGradient(t, p, local_var_index, subdomain_ids); +} + +std::map +SolutionUserObjectBase::discontinuousPointValueGradient( + Real libmesh_dbg_var(t), + Point pt, + const unsigned int local_var_index, + const std::set * subdomain_ids) const +{ + // do the transformations + for (unsigned int trans_num = 0; trans_num < _transformation_order.size(); ++trans_num) + { + if (_transformation_order[trans_num] == "rotation0") + pt = _r0 * pt; + else if (_transformation_order[trans_num] == "translation") + for (const auto i : make_range(Moose::dim)) + pt(i) -= _translation[i]; + else if (_transformation_order[trans_num] == "scale") + for (const auto i : make_range(Moose::dim)) + pt(i) /= _scale[i]; + else if (_transformation_order[trans_num] == "scale_multiplier") + for (const auto i : make_range(Moose::dim)) + pt(i) *= _scale_multiplier[i]; + else if (_transformation_order[trans_num] == "rotation1") + pt = _r1 * pt; + } + + // Extract the value at the current point + std::map map = + evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); + + // Interpolate + if (_file_type == 1 && _interpolate_times) + { + mooseAssert(t == _interpolation_time, + "Time passed into value() must match time at last call to timestepSetup()"); + std::map map2 = + evalMultiValuedMeshFunctionGradient(pt, local_var_index, 1, subdomain_ids); + + if (map.size() != map2.size()) + mooseError( + "In SolutionUserObjectBase::discontinuousPointValue map and map2 have different size"); + + // construct the interpolated map + for (auto & k : map) + { + if (map2.find(k.first) == map2.end()) + mooseError( + "In SolutionUserObjectBase::discontinuousPointValue map and map2 have differing keys"); + RealGradient val = k.second; + RealGradient val2 = map2[k.first]; + map[k.first] = val + (val2 - val) * _interpolation_factor; + } + } + + return map; +} + +Real +SolutionUserObjectBase::directValue(dof_id_type dof_index) const +{ + Real val = (*_serialized_solution)(dof_index); + if (_file_type == 1 && _interpolate_times) + { + Real val2 = (*_serialized_solution2)(dof_index); + val = val + (val2 - val) * _interpolation_factor; + } + return val; +} + +Real +SolutionUserObjectBase::evalMeshFunction(const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids) const +{ + // Storage for mesh function output + DenseVector output; + + // Extract a value from the _mesh_function + { + Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); + if (func_num == 1) + (*_mesh_function)(p, 0.0, output, subdomain_ids); + + // Extract a value from _mesh_function2 + else if (func_num == 2) + (*_mesh_function2)(p, 0.0, output, subdomain_ids); + + else + mooseError("The func_num must be 1 or 2"); + } + + // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated + // outside the domain + if (output.size() == 0) + { + std::ostringstream oss; + p.print(oss); + mooseError("Failed to access the data for variable '", + _system_variables[local_var_index], + "' at point ", + oss.str(), + " in the '", + name(), + "' SolutionUserObjectBase"); + } + return output(local_var_index); +} + +std::map +SolutionUserObjectBase::evalMultiValuedMeshFunction( + const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids) const +{ + // Storage for mesh function output + std::map> temporary_output; + + // Extract a value from the _mesh_function + { + Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); + if (func_num == 1) + _mesh_function->discontinuous_value(p, 0.0, temporary_output, subdomain_ids); + + // Extract a value from _mesh_function2 + else if (func_num == 2) + _mesh_function2->discontinuous_value(p, 0.0, temporary_output, subdomain_ids); + + else + mooseError("The func_num must be 1 or 2"); + } + + // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated + // outside the domain + if (temporary_output.size() == 0) + { + std::ostringstream oss; + p.print(oss); + mooseError("Failed to access the data for variable '", + _system_variables[local_var_index], + "' at point ", + oss.str(), + " in the '", + name(), + "' SolutionUserObjectBase"); + } + + // Fill the actual map that is returned + std::map output; + for (auto & k : temporary_output) + { + mooseAssert( + k.second.size() > local_var_index, + "In SolutionUserObjectBase::evalMultiValuedMeshFunction variable with local_var_index " + << local_var_index << " does not exist"); + output[k.first] = k.second(local_var_index); + } + + return output; +} + +RealGradient +SolutionUserObjectBase::evalMeshFunctionGradient( + const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids) const +{ + // Storage for mesh function output + std::vector output; + + // Extract a value from the _mesh_function + { + Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); + if (func_num == 1) + _mesh_function->gradient(p, 0.0, output, subdomain_ids); + + // Extract a value from _mesh_function2 + else if (func_num == 2) + _mesh_function2->gradient(p, 0.0, output, subdomain_ids); + + else + mooseError("The func_num must be 1 or 2"); + } + + // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated + // outside the domain + if (output.size() == 0) + { + std::ostringstream oss; + p.print(oss); + mooseError("Failed to access the data for variable '", + _system_variables[local_var_index], + "' at point ", + oss.str(), + " in the '", + name(), + "' SolutionUserObjectBase"); + } + return output[local_var_index]; +} + +std::map +SolutionUserObjectBase::evalMultiValuedMeshFunctionGradient( + const Point & p, + const unsigned int local_var_index, + unsigned int func_num, + const std::set * subdomain_ids) const +{ + // Storage for mesh function output + std::map> temporary_output; + + // Extract a value from the _mesh_function + { + Threads::spin_mutex::scoped_lock lock(_solution_user_object_mutex); + if (func_num == 1) + _mesh_function->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids); + + // Extract a value from _mesh_function2 + else if (func_num == 2) + _mesh_function2->discontinuous_gradient(p, 0.0, temporary_output, subdomain_ids); + + else + mooseError("The func_num must be 1 or 2"); + } + + // Error if the data is out-of-range, which will be the case if the mesh functions are evaluated + // outside the domain + if (temporary_output.size() == 0) + { + std::ostringstream oss; + p.print(oss); + mooseError("Failed to access the data for variable '", + _system_variables[local_var_index], + "' at point ", + oss.str(), + " in the '", + name(), + "' SolutionUserObjectBase"); + } + + // Fill the actual map that is returned + std::map output; + for (auto & k : temporary_output) + { + mooseAssert( + k.second.size() > local_var_index, + "In SolutionUserObjectBase::evalMultiValuedMeshFunction variable with local_var_index " + << local_var_index << " does not exist"); + output[k.first] = k.second[local_var_index]; + } + + return output; +} + +const std::vector & +SolutionUserObjectBase::variableNames() const +{ + return _system_variables; +} + +bool +SolutionUserObjectBase::isVariableNodal(const std::string & var_name) const +{ + return std::find(_nodal_variables.begin(), _nodal_variables.end(), var_name) != + _nodal_variables.end(); +} + +Real +SolutionUserObjectBase::scalarValue(Real /*t*/, const std::string & var_name) const +{ + unsigned int var_num = _system->variable_number(var_name); + const DofMap & dof_map = _system->get_dof_map(); + std::vector dofs; + dof_map.SCALAR_dof_indices(dofs, var_num); + // We can handle only FIRST order scalar variables + return directValue(dofs[0]); +} + +void +SolutionUserObjectBase::readBlockIdMapFromExodusII() +{ +#ifdef LIBMESH_HAVE_EXODUS_API + libMesh::ExodusII_IO_Helper & exio_helper = _exodusII_io->get_exio_helper(); + const auto & id_to_block = exio_helper.id_to_block_names; + _block_name_to_id.clear(); + _block_id_to_name.clear(); + for (const auto & it : id_to_block) + { + _block_name_to_id[it.second] = it.first; + _block_id_to_name[it.first] = it.second; + } +#endif +} diff --git a/modules/optimization/include/userobjects/AdjointSolutionUserObject.h b/modules/optimization/include/userobjects/AdjointSolutionUserObject.h index 73eaa9e3ce07..e7bf89190842 100644 --- a/modules/optimization/include/userobjects/AdjointSolutionUserObject.h +++ b/modules/optimization/include/userobjects/AdjointSolutionUserObject.h @@ -9,19 +9,22 @@ #pragma once -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" -class AdjointSolutionUserObject : public SolutionUserObject +class AdjointSolutionUserObject : public SolutionUserObjectBase { public: static InputParameters validParams(); AdjointSolutionUserObject(const InputParameters & parameters); + virtual Real solutionSampleTime() override; + /** * Skipping parent class initialSetup since it will be called in timestepSetup */ virtual void initialSetup() override {} + /** * This will read a the files again if they have been re-written from optimization iteration */ @@ -30,8 +33,7 @@ class AdjointSolutionUserObject : public SolutionUserObject protected: /// Mapping between adjoint simulation time and adjoint simulation time const Real & _reverse_time_end; + /// The system time of the last instance the file was loaded std::time_t _file_mod_time; - /// The forward simulation time for last instance the solution was updated - Real _actual_interpolation_time; }; diff --git a/modules/optimization/src/userobjects/AdjointSolutionUserObject.C b/modules/optimization/src/userobjects/AdjointSolutionUserObject.C index 4998efa50286..63da95366324 100644 --- a/modules/optimization/src/userobjects/AdjointSolutionUserObject.C +++ b/modules/optimization/src/userobjects/AdjointSolutionUserObject.C @@ -9,6 +9,9 @@ #include "AdjointSolutionUserObject.h" +#include "libmesh/mesh_function.h" +#include "libmesh/exodusII_io.h" + #include registerMooseObject("OptimizationApp", AdjointSolutionUserObject); @@ -16,7 +19,7 @@ registerMooseObject("OptimizationApp", AdjointSolutionUserObject); InputParameters AdjointSolutionUserObject::validParams() { - InputParameters params = SolutionUserObject::validParams(); + InputParameters params = SolutionUserObjectBase::validParams(); params.addClassDescription( "Reads a variable from a mesh in one simulation to another specifically for loading forward " "solution in adjoint simulation during inverse optimization."); @@ -34,10 +37,9 @@ AdjointSolutionUserObject::validParams() } AdjointSolutionUserObject::AdjointSolutionUserObject(const InputParameters & parameters) - : SolutionUserObject(parameters), + : SolutionUserObjectBase(parameters), _reverse_time_end(getParam("reverse_time_end")), - _file_mod_time(std::numeric_limits::min()), - _actual_interpolation_time(0.0) + _file_mod_time(std::numeric_limits::min()) { if (!MooseUtils::hasExtension(_mesh_file, "e", /*strip_exodus_ext =*/true)) paramError("mesh", @@ -45,6 +47,12 @@ AdjointSolutionUserObject::AdjointSolutionUserObject(const InputParameters & par "solution is written and read from exodus format."); } +Real +AdjointSolutionUserObject::solutionSampleTime() +{ + return _reverse_time_end - _t + _dt; +} + void AdjointSolutionUserObject::timestepSetup() { @@ -63,24 +71,11 @@ AdjointSolutionUserObject::timestepSetup() _es2.reset(); // Read the exodus file - SolutionUserObject::initialSetup(); + SolutionUserObjectBase::initialSetup(); // Make sure to communicate what solution was actually loaded - _actual_interpolation_time = 0.0; _interpolation_time = -_dt; } - // Update time interpolation for ExodusII solution - const Real actual_time = _reverse_time_end - _t + _dt; - if (_actual_interpolation_time != actual_time) - { - // This will make sure the solutions are updated for the current forward time - _interpolation_time = _actual_interpolation_time; - // Update based on the forward simulation time - updateExodusTimeInterpolation(actual_time); - // This avoids errors in public functions when the time is given as the adjoint simulation time - _interpolation_time = _t; - // Ensures we only update the solution when necessary - _actual_interpolation_time = actual_time; - } + SolutionUserObjectBase::timestepSetup(); } diff --git a/modules/phase_field/src/userobjects/SolutionRasterizer.C b/modules/phase_field/src/userobjects/SolutionRasterizer.C index ee27603eeb43..5942485ac4eb 100644 --- a/modules/phase_field/src/userobjects/SolutionRasterizer.C +++ b/modules/phase_field/src/userobjects/SolutionRasterizer.C @@ -10,6 +10,8 @@ #include "SolutionRasterizer.h" #include +#include "libmesh/mesh_function.h" +#include "libmesh/exodusII_io.h" registerMooseObject("PhaseFieldApp", SolutionRasterizer); diff --git a/test/include/postprocessors/TestDiscontinuousValuePP.h b/test/include/postprocessors/TestDiscontinuousValuePP.h index b7b4aae37f14..b8fb198c9db2 100644 --- a/test/include/postprocessors/TestDiscontinuousValuePP.h +++ b/test/include/postprocessors/TestDiscontinuousValuePP.h @@ -12,7 +12,7 @@ #include "GeneralPostprocessor.h" // Forward Declarations -class SolutionUserObject; +class SolutionUserObjectBase; /** * Compute the value of a variable or the gradient at a specified location. @@ -52,5 +52,5 @@ class TestDiscontinuousValuePP : public GeneralPostprocessor unsigned int _gradient_component; /// Pointer to SolutionUserObject containing the solution of interest - const SolutionUserObject * _solution_object_ptr; + const SolutionUserObjectBase * _solution_object_ptr; }; diff --git a/test/src/postprocessors/TestDiscontinuousValuePP.C b/test/src/postprocessors/TestDiscontinuousValuePP.C index c67ca976e306..2effbf3b14a7 100644 --- a/test/src/postprocessors/TestDiscontinuousValuePP.C +++ b/test/src/postprocessors/TestDiscontinuousValuePP.C @@ -9,7 +9,7 @@ // MOOSE includes #include "TestDiscontinuousValuePP.h" -#include "SolutionUserObject.h" +#include "SolutionUserObjectBase.h" registerMooseObject("MooseTestApp", TestDiscontinuousValuePP); @@ -46,7 +46,7 @@ TestDiscontinuousValuePP::TestDiscontinuousValuePP(const InputParameters & param void TestDiscontinuousValuePP::initialSetup() { - _solution_object_ptr = &getUserObject("solution"); + _solution_object_ptr = &getUserObject("solution"); } Real diff --git a/test/tests/auxkernels/solution_aux/gold/solution_aux_exodus_interp_out2.e b/test/tests/auxkernels/solution_aux/gold/solution_aux_exodus_interp_out2.e new file mode 100644 index 000000000000..3967749bb955 Binary files /dev/null and b/test/tests/auxkernels/solution_aux/gold/solution_aux_exodus_interp_out2.e differ diff --git a/test/tests/auxkernels/solution_aux/tests b/test/tests/auxkernels/solution_aux/tests index 3e64f0ec14d1..27a873f1b54e 100644 --- a/test/tests/auxkernels/solution_aux/tests +++ b/test/tests/auxkernels/solution_aux/tests @@ -110,6 +110,15 @@ "temporal interpolation." [] + [exodus_time_transformation] + type = 'Exodiff' + input = 'solution_aux_exodus_interp.i' + exodiff = 'solution_aux_exodus_interp_out2.e' + cli_args = 'UserObjects/soln/time_transformation=t/2 Outputs/file_base=solution_aux_exodus_interp_out2' + requirement = "The SolutionAux object shall be capable of setting an auxiliary variable with " + "temporal interpolation and a time transformation function." + [] + [exodus_interp_restart] requirement = "The system shall be capable of initializing an auxiliary variable from an " "existing solution that"