Line data Source code
1 : // Copyright 2024, UChicago Argonne, LLC
2 : // All Rights Reserved
3 : // Software Name: NEML2 -- the New Engineering material Model Library, version 2
4 : // By: Argonne National Laboratory
5 : // OPEN SOURCE LICENSE (MIT)
6 : //
7 : // Permission is hereby granted, free of charge, to any person obtaining a copy
8 : // of this software and associated documentation files (the "Software"), to deal
9 : // in the Software without restriction, including without limitation the rights
10 : // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
11 : // copies of the Software, and to permit persons to whom the Software is
12 : // furnished to do so, subject to the following conditions:
13 : //
14 : // The above copyright notice and this permission notice shall be included in
15 : // all copies or substantial portions of the Software.
16 : //
17 : // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
18 : // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
19 : // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
20 : // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
21 : // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
22 : // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
23 : // THE SOFTWARE.
24 :
25 : #include "neml2/models/SymmetricHermiteInterpolation.h"
26 : #include "neml2/tensors/functions/clamp.h"
27 : #include "neml2/tensors/functions/where.h"
28 : #include "neml2/tensors/assertions.h"
29 :
30 : namespace neml2
31 : {
32 : register_NEML2_object(SymmetricHermiteInterpolation);
33 :
34 : OptionSet
35 2 : SymmetricHermiteInterpolation::expected_options()
36 : {
37 2 : OptionSet options = Model::expected_options();
38 2 : options.doc() =
39 : "Define the symmetric Hermite interpolation function, taking the form of \\f$ (3c^2-4c^3) "
40 : "\\f$ for \\f$ 0 le c le 0.5 \\f$; \\f$ 3(1-c)^2 - 4(1-c)^3 \\f$ for \\f$ 0.5 le c le 1 "
41 : "\\f$, and 0.0 otherwise. Here, \\f$ c = \\frac{x-x_l}{x_h-x_l} \\f$ where \\f$x_l\\f$ and "
42 2 : "\\f$x_h\\f$ are the lower and upper bound for rescaling the input argument.";
43 :
44 4 : options.set_input("argument");
45 4 : options.set("argument").doc() = "Argument of the smooth step function";
46 :
47 4 : options.set<bool>("define_second_derivatives") = true;
48 :
49 4 : options.set_output("value");
50 4 : options.set("value").doc() = "Value of the smooth step function";
51 :
52 4 : options.set_buffer<TensorName<Scalar>>("lower_bound");
53 4 : options.set("lower_bound").doc() = "Lower bound of the argument";
54 :
55 4 : options.set_buffer<TensorName<Scalar>>("upper_bound");
56 2 : options.set("upper_bound").doc() = "Upper bound of the argument";
57 :
58 2 : return options;
59 0 : }
60 :
61 1 : SymmetricHermiteInterpolation::SymmetricHermiteInterpolation(const OptionSet & options)
62 : : Model(options),
63 1 : _x(declare_input_variable<Scalar>("argument")),
64 1 : _y(declare_output_variable<Scalar>("value")),
65 4 : _x0(declare_buffer<Scalar>("lb", "lower_bound")),
66 5 : _x1(declare_buffer<Scalar>("ub", "upper_bound"))
67 : {
68 1 : }
69 :
70 : void
71 3 : SymmetricHermiteInterpolation::set_value(bool out, bool dout_din, bool d2out_din2)
72 : {
73 3 : const auto eps = machine_precision(_x.scalar_type()).toDouble();
74 3 : const auto x = clamp((_x - _x0) / (_x1 - _x0), eps, 1.0 - eps);
75 :
76 3 : if (out)
77 : {
78 1 : auto f_xl = 3 * x * x - 4 * x * x * x;
79 1 : auto f_xh = 3 * (1 - x) * (1 - x) - 4 * (1 - x) * (1 - x) * (1 - x);
80 1 : _y = where(x < 0.5, f_xl, f_xh);
81 1 : }
82 :
83 3 : if (dout_din)
84 : {
85 1 : auto df_xl = 6 * x - 12 * x * x;
86 1 : auto df_xh = -6 * (1 - x) + 12 * (1 - x) * (1 - x);
87 :
88 1 : _y.d(_x) = where(x < 0.5, df_xl, df_xh);
89 1 : }
90 :
91 3 : if (d2out_din2)
92 : {
93 1 : auto df2_xl = 6 - 24 * x;
94 1 : auto df2_xh = 6 - 24 * (1 - x);
95 :
96 1 : const auto zeromask = Scalar(at::logical_and(at::lt(x, 1.0 - eps), at::gt(x, eps)));
97 1 : _y.d(_x, _x) = zeromask * where(x < 0.5, df2_xl, df2_xh);
98 1 : }
99 3 : }
100 : } // namespace neml2
|