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/ScalarMultiplication.h"
26 : #include "neml2/tensors/Scalar.h"
27 : #include "neml2/misc/assertions.h"
28 : #include "neml2/tensors/functions/pow.h"
29 :
30 : namespace neml2
31 : {
32 : register_NEML2_object(ScalarMultiplication);
33 : OptionSet
34 2 : ScalarMultiplication::expected_options()
35 : {
36 :
37 2 : OptionSet options = Model::expected_options();
38 2 : options.doc() = "Calculate the multiplication (product) of multiple Scalar variable with a "
39 : "constant coefficient. Using reciprocal, one can have the reciprocity of "
40 2 : "variable 'a', aka. '1/a'";
41 :
42 4 : options.set<std::vector<VariableName>>("from_var");
43 4 : options.set("from_var").doc() = "Scalar variables to be multiplied";
44 :
45 4 : options.set_output("to_var");
46 4 : options.set("to_var").doc() = "The multiplicative product";
47 :
48 8 : options.set_parameter<TensorName<Scalar>>("coefficient") = {TensorName<Scalar>("1")};
49 2 : options.set("coefficient").doc() = "The coefficient multiply to the final product";
50 :
51 6 : options.set<std::vector<bool>>("reciprocal") = {false};
52 6 : options.set("reciprocal").doc() =
53 : "List of boolens, one for each variable, in which the reciprocity of the corresponding "
54 : "variable is taken. When the length of this list is 1, the same reciprocal condition applies "
55 2 : "to all variables.";
56 :
57 2 : options.set<bool>("define_second_derivatives") = true;
58 :
59 2 : return options;
60 0 : }
61 :
62 6 : ScalarMultiplication::ScalarMultiplication(const OptionSet & options)
63 : : Model(options),
64 6 : _to(declare_output_variable<Scalar>("to_var")),
65 36 : _A(declare_parameter<Scalar>("A", "coefficient"))
66 : {
67 36 : for (const auto & fv : options.get<std::vector<VariableName>>("from_var"))
68 30 : _from.push_back(&declare_input_variable<Scalar>(fv));
69 :
70 6 : _inv = options.get<std::vector<bool>>("reciprocal");
71 6 : neml_assert(_inv.size() == 1 || _inv.size() == _from.size(),
72 : "Expected 1 or ",
73 6 : _from.size(),
74 : " entries in reciprocal, got ",
75 6 : _inv.size(),
76 : ".");
77 :
78 : // Expand the list of booleans to match the number of coefficients
79 6 : if (_inv.size() == 1)
80 2 : _inv = std::vector<bool>(_from.size(), _inv[0]);
81 6 : }
82 :
83 : void
84 18 : ScalarMultiplication::set_value(bool out, bool dout_din, bool d2out_din2)
85 : {
86 18 : if (out)
87 : {
88 6 : auto r = _inv[0] ? _A / (*_from[0]) : _A * (*_from[0]);
89 24 : for (std::size_t i = 1; i < _from.size(); i++)
90 : {
91 18 : if (_inv[i])
92 9 : r = r / (*_from[i]);
93 : else
94 9 : r = r * (*_from[i]);
95 : }
96 6 : _to = r;
97 6 : }
98 :
99 18 : if (dout_din)
100 : {
101 30 : for (std::size_t i = 0; i < _from.size(); i++)
102 24 : if (_from[i]->is_dependent())
103 : {
104 24 : auto r = _inv[i] ? -_A / (*_from[i]) / (*_from[i]) : _A;
105 120 : for (std::size_t j = 0; j < _from.size(); j++)
106 96 : if (i != j)
107 : {
108 72 : if (_inv[j])
109 36 : r = r / (*_from[j]);
110 : else
111 36 : r = r * (*_from[j]);
112 : }
113 24 : _to.d(*_from[i]) = r;
114 24 : }
115 : }
116 :
117 18 : if (d2out_din2)
118 : {
119 30 : for (std::size_t i = 0; i < _from.size(); i++)
120 24 : if (_from[i]->is_dependent())
121 : {
122 24 : auto p = (_inv[i] ? -1.0 : 1.0);
123 120 : for (std::size_t j = 0; j < _from.size(); j++)
124 : {
125 96 : if (_from[j]->is_dependent())
126 : {
127 96 : auto q = (_inv[j] ? -1.0 : 1.0);
128 96 : if (i != j)
129 : {
130 72 : auto r = _A * p * pow(*_from[i], (p - 1)) * q * pow(*_from[j], (q - 1));
131 360 : for (std::size_t k = 0; k < _from.size(); k++)
132 288 : if (k != i && k != j)
133 144 : r = r * (_inv[k] ? 1. / (*_from[k]) : (*_from[k]));
134 72 : _to.d(*_from[i], *_from[j]) = r;
135 72 : }
136 24 : else if (_inv[i])
137 : {
138 12 : auto r = _A * p * (p - 1) * pow(*_from[i], (p - 2));
139 60 : for (std::size_t k = 0; k < _from.size(); k++)
140 48 : if (k != i)
141 36 : r = r * (_inv[k] ? 1. / (*_from[k]) : (*_from[k]));
142 12 : _to.d(*_from[i], *_from[j]) = r;
143 12 : }
144 : }
145 : }
146 : }
147 : }
148 18 : }
149 : } // namespace neml2
|
