Skip to content

Commit

Permalink
Fix pretty egregious issue with storing the history object
Browse files Browse the repository at this point in the history
  • Loading branch information
@reverendbedford
reverendbedford committed Nov 16, 2022
1 parent 607703d commit 3e926c4
Show file tree
Hide file tree
Showing 5 changed files with 153 additions and 145 deletions.
11 changes: 7 additions & 4 deletions include/history.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,6 +281,8 @@ class NEML_EXPORT HistoryNEMLObject: public NEMLObject {
virtual void init_hist(History & h) const = 0;
/// This should be replaced at some point
virtual size_t nhist() const;
/// Dangerous nhist which assumes you cached
virtual size_t nh() const;

/// Setup a flat vector history
virtual void init_store(double * const h) const;
Expand All @@ -299,17 +301,18 @@ class NEML_EXPORT HistoryNEMLObject: public NEMLObject {
/// Cache the history to avoid recreating it every time
void cache_history_();

/// Quickly setup history
History gather_history_(double * data) const;
History gather_history_(const double * data) const;
History gather_blank_history_() const;
/// Quickly setup history
History gather_history_(double * data) const;
History gather_history_(const double * data) const;
History gather_blank_history_() const;

protected:
std::string prefix_;
History stored_hist_;

private:
bool cached_;
size_t ncache_;
};

} // namespace neml
Expand Down
38 changes: 18 additions & 20 deletions src/general_flow.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ TVPFlowRule::TVPFlowRule(ParameterSet & params) :
elastic_(params.get_object_parameter<LinearElasticModel>("elastic")),
flow_(params.get_object_parameter<ViscoPlasticFlowRule>("flow"))
{

cache_history_();
}

std::string TVPFlowRule::type()
Expand Down Expand Up @@ -156,7 +156,7 @@ void TVPFlowRule::ds_da(const double * const s, const double * const alpha,
double yv;
flow_->y(s, alpha, T, yv);

int sz = 6 * nhist();
int sz = 6 * nh();

std::vector<double> workv(sz);
double * work = &workv[0];
Expand All @@ -167,10 +167,10 @@ void TVPFlowRule::ds_da(const double * const s, const double * const alpha,

double t1[6];
flow_->g(s, alpha, T, t1);
std::vector<double> t2v(nhist());
std::vector<double> t2v(nh());
double * t2 = &t2v[0];
flow_->dy_da(s, alpha, T, t2);
outer_update_minus(t1, 6, t2, nhist(), work);
outer_update_minus(t1, 6, t2, nh(), work);

std::vector<double> t3v(sz);
double * t3 = &t3v[0];
Expand All @@ -187,9 +187,7 @@ void TVPFlowRule::ds_da(const double * const s, const double * const alpha,
double C[36];
elastic_->C(T, C);

mat_mat(6, nhist(), 6, C, work, d_sdot);


mat_mat(6, nh(), 6, C, work, d_sdot);
}

void TVPFlowRule::ds_de(const double * const s, const double * const alpha,
Expand All @@ -209,15 +207,15 @@ void TVPFlowRule::a(const double * const s, const double * const alpha,
flow_->y(s, alpha, T, dg);

flow_->h(s, alpha, T, adot);
for (size_t i=0; i<nhist(); i++) adot[i] *= dg;
for (size_t i=0; i<nh(); i++) adot[i] *= dg;

std::vector<double> tempv(nhist());
std::vector<double> tempv(nh());
double * temp = &tempv[0];
flow_->h_temp(s, alpha, T, temp);
for (size_t i=0; i<nhist(); i++) adot[i] += temp[i] * Tdot;
for (size_t i=0; i<nh(); i++) adot[i] += temp[i] * Tdot;

flow_->h_time(s, alpha, T, temp);
for (size_t i=0; i<nhist(); i++) adot[i] += temp[i];
for (size_t i=0; i<nh(); i++) adot[i] += temp[i];


}
Expand All @@ -230,19 +228,19 @@ void TVPFlowRule::da_ds(const double * const s, const double * const alpha,
double dg;
flow_->y(s, alpha, T, dg);

int sz = nhist() * 6;
int sz = nh() * 6;

flow_->dh_ds(s, alpha, T, d_adot);
for (int i=0; i<sz; i++) d_adot[i] *= dg;

std::vector<double> t1v(nhist());
std::vector<double> t1v(nh());
double * t1 = &t1v[0];
flow_->h(s, alpha, T, t1);

double t2[6];
flow_->dy_ds(s, alpha, T, t2);

outer_update(t1, nhist(), t2, 6, d_adot);
outer_update(t1, nh(), t2, 6, d_adot);

std::vector<double> t3v(sz);
double * t3 = &t3v[0];
Expand All @@ -263,21 +261,21 @@ void TVPFlowRule::da_da(const double * const s, const double * const alpha,
double dg;
flow_->y(s, alpha, T, dg);

int nh = nhist();
int sz = nh * nh;
int nhi = nh();
int sz = nhi * nhi;

flow_->dh_da(s, alpha, T, d_adot);
for (int i=0; i<sz; i++) d_adot[i] *= dg;

std::vector<double> t1v(nh);
std::vector<double> t1v(nhi);
double * t1 = &t1v[0];
flow_->h(s, alpha, T, t1);

std::vector<double> t2v(nh);
std::vector<double> t2v(nhi);
double * t2 = &t2v[0];
flow_->dy_da(s, alpha, T, t2);

outer_update(t1, nh, t2, nh, d_adot);
outer_update(t1, nhi, t2, nhi, d_adot);

std::vector<double> t3v(sz);
double * t3 = &t3v[0];
Expand All @@ -294,7 +292,7 @@ void TVPFlowRule::da_de(const double * const s, const double * const alpha,
double Tdot,
double * const d_adot)
{
std::fill(d_adot, d_adot+(nhist()*6), 0.0);
std::fill(d_adot, d_adot+(nh()*6), 0.0);

}

Expand Down
36 changes: 18 additions & 18 deletions src/hardening.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -973,7 +973,7 @@ ChabocheVoceRecovery::ChabocheVoceRecovery(ParameterSet & params) :
a_( params.get_object_parameter_vector<Interpolate>("a")),
noniso_(params.get_parameter<bool>("noniso"))
{

cache_history_();
}

std::string ChabocheVoceRecovery::type()
Expand Down Expand Up @@ -1045,7 +1045,7 @@ void ChabocheVoceRecovery::q(const double * const alpha, double T, double * cons

void ChabocheVoceRecovery::dq_da(const double * const alpha, double T, double * const qv) const
{
std::fill(qv, qv+(ninter()*nhist()), 0.0);
std::fill(qv, qv+(ninter()*nh()), 0.0);

// Isotropic part
qv[0] = -1.0;
Expand All @@ -1055,7 +1055,7 @@ void ChabocheVoceRecovery::dq_da(const double * const alpha, double T, double *

for (size_t i=0; i<n; i++) {
for (size_t j=0; j<6; j++) {
qv[CINDEX((j+1),(1+i*6+j), nhist())] = 1.0;
qv[CINDEX((j+1),(1+i*6+j), nh())] = 1.0;
}
}
}
Expand Down Expand Up @@ -1090,7 +1090,7 @@ void ChabocheVoceRecovery::h(const double * const s, const double * const alpha,
void ChabocheVoceRecovery::dh_ds(const double * const s, const double * const alpha, double T,
double * const dhv) const
{
std::fill(dhv, dhv + nhist()*6, 0.0);
std::fill(dhv, dhv + nh()*6, 0.0);

std::vector<double> c = eval_vector(c_, T);

Expand Down Expand Up @@ -1146,9 +1146,9 @@ void ChabocheVoceRecovery::dh_da(const double * const s, const double * const al
double * const dhv) const
{
// Again, there should be no earthly reason the compiler can't do this
size_t nh = nhist();
size_t nhi = nh();

std::fill(dhv, dhv + nh*nh, 0.0);
std::fill(dhv, dhv + nhi*nhi, 0.0);

// Isotropic contribution
dhv[0] = -theta0_->value(T) / Rmax_->value(T) * std::sqrt(2.0/3.0);
Expand Down Expand Up @@ -1180,7 +1180,7 @@ void ChabocheVoceRecovery::dh_da(const double * const s, const double * const al
// Fill in the gamma part
for (size_t i=0; i<n_; i++) {
for (size_t j=0; j<6; j++) {
dhv[CINDEX((1+i*6+j),(1+i*6+j),nh)] -= sqrt(2.0/3.0) * gmodels_[i]->gamma(alpha[0], T);
dhv[CINDEX((1+i*6+j),(1+i*6+j),nhi)] -= sqrt(2.0/3.0) * gmodels_[i]->gamma(alpha[0], T);
}
}

Expand All @@ -1189,7 +1189,7 @@ void ChabocheVoceRecovery::dh_da(const double * const s, const double * const al
for (size_t i=0; i<6; i++) {
for (size_t bj=0; bj<n_; bj++) {
for (size_t j=0; j<6; j++) {
dhv[CINDEX((1+bi*6+i),(1+bj*6+j),nh)] -= 2.0 / 3.0 * c[bi] *
dhv[CINDEX((1+bi*6+i),(1+bj*6+j),nhi)] -= 2.0 / 3.0 * c[bi] *
ss[CINDEX(i,j,6)];
}
}
Expand All @@ -1199,7 +1199,7 @@ void ChabocheVoceRecovery::dh_da(const double * const s, const double * const al
// Fill in the alpha part
for (size_t i=0; i<n_; i++) {
for (size_t j=0; j<6; j++) {
dhv[CINDEX((1+i*6+j),0,nhist())] = -sqrt(2.0/3.0) *
dhv[CINDEX((1+i*6+j),0,nhi)] = -sqrt(2.0/3.0) *
gmodels_[i]->dgamma(alpha[0], T) * alpha[1+i*6+j];
}
}
Expand All @@ -1209,7 +1209,7 @@ void ChabocheVoceRecovery::dh_da(const double * const s, const double * const al
void ChabocheVoceRecovery::h_time(const double * const s, const double * const alpha,
double T, double * const hv) const
{
std::fill(hv, hv+nhist(), 0.0);
std::fill(hv, hv+nh(), 0.0);

// Isotropic recovery term
hv[0] = r1_->value(T) * (Rmin_->value(T) - alpha[0]
Expand All @@ -1235,7 +1235,7 @@ void ChabocheVoceRecovery::h_time(const double * const s, const double * const a
void ChabocheVoceRecovery::dh_ds_time(const double * const s, const double * const alpha,
double T, double * const dhv) const
{
std::fill(dhv, dhv+nhist()*6, 0.0);
std::fill(dhv, dhv+nh()*6, 0.0);

// Also return if relax

Expand All @@ -1244,7 +1244,7 @@ void ChabocheVoceRecovery::dh_ds_time(const double * const s, const double * con
void ChabocheVoceRecovery::dh_da_time(const double * const s, const double * const alpha,
double T, double * const dhv) const
{
std::fill(dhv, dhv+nhist()*nhist(), 0.0);
std::fill(dhv, dhv+nh()*nh(), 0.0);

dhv[0] = std::copysign(r1_->value(T) * r2_->value(T) *
std::pow(std::fabs(Rmin_->value(T) - alpha[0]), r2_->value(T) -
Expand All @@ -1253,7 +1253,7 @@ void ChabocheVoceRecovery::dh_da_time(const double * const s, const double * con
std::vector<double> A = eval_vector(A_, T);
std::vector<double> a = eval_vector(a_, T);

size_t nh = nhist();
size_t nhi = nh();
size_t n = n_;

double XX[36];
Expand All @@ -1276,7 +1276,7 @@ void ChabocheVoceRecovery::dh_da_time(const double * const s, const double * con
else {
d = 0.0;
}
dhv[CINDEX(ia,ib,nh)] = -A[i] * pow( sqrt(3.0/2.0) * nXi, a[i]-1.0) * (
dhv[CINDEX(ia,ib,nhi)] = -A[i] * pow( sqrt(3.0/2.0) * nXi, a[i]-1.0) * (
d + (a[i] - 1.0) * XX[CINDEX(j,k,6)]);
}
}
Expand All @@ -1287,7 +1287,7 @@ void ChabocheVoceRecovery::dh_da_time(const double * const s, const double * con
void ChabocheVoceRecovery::h_temp(const double * const s, const double * const alpha, double T,
double * const hv) const
{
std::fill(hv, hv+nhist(), 0.0);
std::fill(hv, hv+nh(), 0.0);

std::vector<double> c = eval_vector(c_, T);
std::vector<double> dc = eval_deriv_vector(c_, T);
Expand All @@ -1303,13 +1303,13 @@ void ChabocheVoceRecovery::h_temp(const double * const s, const double * const a
void ChabocheVoceRecovery::dh_ds_temp(const double * const s, const double * const alpha, double T,
double * const dhv) const
{
std::fill(dhv, dhv+nhist()*6, 0.0);
std::fill(dhv, dhv+nh()*6, 0.0);
}

void ChabocheVoceRecovery::dh_da_temp(const double * const s, const double * const alpha, double T,
double * const dhv) const
{
std::fill(dhv, dhv+nhist()*nhist(), 0.0);
std::fill(dhv, dhv+nh()*nh(), 0.0);

std::vector<double> c = eval_vector(c_, T);
std::vector<double> dc = eval_deriv_vector(c_, T);
Expand All @@ -1318,7 +1318,7 @@ void ChabocheVoceRecovery::dh_da_temp(const double * const s, const double * con
if (c[i] == 0.0) continue;
for (size_t j=0; j<6; j++) {
size_t ci = 1 + i*6 + j;
dhv[CINDEX(ci,ci,nhist())] = - sqrt(2.0/3.0) * dc[i] / c[i];
dhv[CINDEX(ci,ci,nh())] = - sqrt(2.0/3.0) * dc[i] / c[i];
}
}
}
Expand Down
9 changes: 7 additions & 2 deletions src/history.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -432,19 +432,23 @@ std::vector<std::string> History::formatted_names() const
HistoryNEMLObject::HistoryNEMLObject(ParameterSet & params) :
NEMLObject(params), prefix_(""), cached_(false)
{

}

size_t HistoryNEMLObject::nhist() const
{
if (cached_)
return stored_hist_.size();
return ncache_;

History h;
populate_hist(h);
return h.size();
}

size_t HistoryNEMLObject::nh() const
{
return ncache_;
}

void HistoryNEMLObject::init_store(double * const h) const
{
History hist(h);
Expand Down Expand Up @@ -480,6 +484,7 @@ std::string HistoryNEMLObject::dprefix(std::string a, std::string b) const
void HistoryNEMLObject::cache_history_()
{
populate_hist(stored_hist_);
ncache_ = nhist();
cached_ = true;
}

Expand Down
Loading

0 comments on commit 3e926c4

Please sign in to comment.