From c49aa7b99a09aeb6088acb31725d7597ba08e2ee Mon Sep 17 00:00:00 2001 From: Sanath Date: Tue, 6 Oct 2026 17:35:51 +0200 Subject: [PATCH 1/8] first version --- CMakeLists.txt | 2 +- docs/NEML2.md | 2 + include/fans_plugin.h | 6 +++ include/matmodel.h | 27 ++++++++++++- include/solver.h | 56 +++++++++++++++++--------- include/solverCG.h | 30 ++++---------- include/solverFP.h | 6 +-- src/neml2_plugin.cpp | 72 ++++++++++++++++++++++++---------- src/plugin_material.cpp | 40 ++++++++++++++++--- src/tangent.cpp | 87 +++++++++++++++++++++++++++++++++++++++++ 10 files changed, 256 insertions(+), 72 deletions(-) create mode 100644 src/tangent.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 5b41a00..24bae8c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -156,7 +156,7 @@ target_include_directories(FANS_main PRIVATE "${PROJECT_BINARY_DIR}/include") # SOURCES # ############################################################################## -target_sources(FANS_FANS PRIVATE src/reader.cpp src/plugin_material.cpp) +target_sources(FANS_FANS PRIVATE src/reader.cpp src/plugin_material.cpp src/tangent.cpp) target_sources(FANS_main PRIVATE src/main.cpp) # ############################################################################## diff --git a/docs/NEML2.md b/docs/NEML2.md index 90e5974..390d648 100644 --- a/docs/NEML2.md +++ b/docs/NEML2.md @@ -32,6 +32,8 @@ pixi run -e dev-neml2 neml2-compile model.i "cpu cuda" This compiles the model block named like the file (`[model]`) for both devices into `compiled_models/model`, next to `model.i`. A third argument lists Python files the model needs, e.g. `"a.py b.py"`. +For `homogenized_tangent`, compile the model's flux-gradient derivative, e.g. `neml2-compile ... -d state/S:forces/E` for small strain. The current sensitivity solve uses PCG and therefore requires a symmetric, CG-compatible algorithmic tangent. + ## Input file ```json diff --git a/include/fans_plugin.h b/include/fans_plugin.h index 54a344e..5ea9fdd 100644 --- a/include/fans_plugin.h +++ b/include/fans_plugin.h @@ -58,6 +58,12 @@ int fans_plugin_evaluate(FANSPluginModel *model, size_t n_points, double t_old, const int *voxel, const int *phase, const double *history_old, double *flux, double *history_new, char *err, size_t errlen); +/* Consistent flux-gradient tangent [n_points][n_str][n_str]. */ +int fans_plugin_tangent(FANSPluginModel *model, size_t n_points, double t_old, double t, + const double *gradient, const int *voxel, const int *phase, + const double *history_old, const double *history_current, + double *tangent, char *err, size_t errlen); + void fans_plugin_free(FANSPluginModel *model); #ifdef __cplusplus diff --git a/include/matmodel.h b/include/matmodel.h index ddd7e5c..3222753 100644 --- a/include/matmodel.h +++ b/include/matmodel.h @@ -54,7 +54,32 @@ class Matmodel { } virtual void evaluate_batch(const vector &elems, const unsigned short *phase, const double *ue, double *sig_gp) {} - + virtual bool has_consistent_tangent() const { return false; } + virtual void compute_tangent_field(ptrdiff_t, const unsigned short *, vector &) + { + throw std::logic_error("Material tangent not implemented."); + } + Matrix &element_tangent_residual( + const Matrix &ue, const Matrix ¯o, const double *tangent) + { + res_e.setZero(); + for (int gp = 0; gp < n_gp; ++gp) { + const Matrix deps = B_int[gp] * ue + macro; + const Matrix dsig = Map>(tangent + gp * n_str * n_str) * deps; + res_e.noalias() += B_int[gp].transpose() * dsig * v_e / n_gp; + } + return res_e; + } + Matrix element_tangent_stress( + const Matrix &ue, const Matrix ¯o, const double *tangent) const + { + Matrix stress = Matrix::Zero(); + for (int gp = 0; gp < n_gp; ++gp) { + const Matrix deps = B_int[gp] * ue + macro; + stress.noalias() += Map>(tangent + gp * n_str * n_str) * deps / n_gp; + } + return stress; + } // A model without the element stiffnesses of a LinearModel can still say that its // flux is linear in the gradient: then the linear CG applies. virtual bool is_linear() const diff --git a/include/solver.h b/include/solver.h index b6f8789..dadfb8c 100644 --- a/include/solver.h +++ b/include/solver.h @@ -40,7 +40,6 @@ class Solver : private MixedBCController { RealArray v_r_real; // can't do the "classname()" intialization here, and Map doesn't have a default constructor RealArray v_u_real; - Map rhat; ArrayXd err_all; //!< Absolute error history Matrix fundamentalSolution; @@ -63,9 +62,9 @@ class Solver : private MixedBCController { void postprocess(Reader &reader, int load_idx, int time_idx); //!< Computes Strain and stress - void convolution(); + void apply_preconditioner(double *input, double *output); + double dotProduct(RealArray &a, RealArray &b); double compute_error(RealArray &r, const std::string &details = {}); - void CreateFFTWPlans(double *in, fftw_complex *transformed, double *out); VectorXd homogenized_strain; VectorXd homogenized_stress; @@ -101,14 +100,19 @@ class Solver : private MixedBCController { (svd.singularValues().array() > tolerance).select(svd.singularValues().array().inverse(), 0).matrix().asDiagonal() * svd.matrixU().adjoint(); } + double *fft_workspace; + fftw_plan planfft{nullptr}, planifft{nullptr}; + void createFFTWPlans(); void computeFundamentalSolution(); protected: - fftw_plan planfft, planifft; clock_t fft_time, buftime; size_t iter; }; +template +MatrixXd compute_consistent_homogenized_tangent(Solver &solver); + template Solver::Solver(Reader &reader, MaterialManager *matmgr) : reader(reader), @@ -135,8 +139,8 @@ Solver::Solver(Reader &reader, MaterialManager * v_u_real(v_u, n_z * howmany, local_n0 * n_y, OuterStride<>(n_z * howmany)), v_u_prev(fftw_alloc_real(local_n0 * n_y * n_z * howmany)), - rhat((std::complex *) v_r, local_n1 * n_x * (n_z / 2 + 1) * howmany), // actual initialization is below - buffer_padding(fftw_alloc_real(n_y * (n_z + 2) * howmany)) + buffer_padding(fftw_alloc_real(n_y * (n_z + 2) * howmany)), + fft_workspace(fftw_alloc_real(reader.alloc_local * 2)) { v_u_real.setZero(); for (ptrdiff_t i = local_n0 * n_y * n_z * howmany; i < (local_n0 + 1) * n_y * n_z * howmany; i++) { @@ -147,6 +151,7 @@ Solver::Solver(Reader &reader, MaterialManager * matmanager->initialize_internal_variables(local_n0 * n_y * n_z, matmanager->models[0]->n_gp); computeFundamentalSolution(); + createFFTWPlans(); } template @@ -208,7 +213,7 @@ void Solver::computeFundamentalSolution() } template -void Solver::CreateFFTWPlans(double *in, fftw_complex *transformed, double *out) +void Solver::createFFTWPlans() { int rank = 3; ptrdiff_t iblock = FFTW_MPI_DEFAULT_BLOCK; @@ -222,11 +227,8 @@ void Solver::CreateFFTWPlans(double *in, fftw_complex *transform // But, according to https://fftw.org/doc/MPI-Plan-Creation.html the BLOCK sizes must be the same: // "These must be the same block sizes as were passed to the corresponding ‘local_size’ function" const ptrdiff_t n[3] = {n_x, n_y, n_z}; - planfft = fftw_mpi_plan_many_dft_r2c(rank, n, howmany, iblock, oblock, in, transformed, communicator, FFTW_MEASURE | FFTW_MPI_TRANSPOSED_OUT); - planifft = fftw_mpi_plan_many_dft_c2r(rank, n, howmany, iblock, oblock, transformed, out, communicator, FFTW_MEASURE | FFTW_MPI_TRANSPOSED_IN); - - // see https://eigen.tuxfamily.org/dox/group__TutorialMapClass.html#title3 - new (&rhat) Map((std::complex *) transformed, local_n1 * n_x * (n_z / 2 + 1) * howmany); + planfft = fftw_mpi_plan_many_dft_r2c(rank, n, howmany, iblock, oblock, v_r, reinterpret_cast(fft_workspace), communicator, FFTW_MEASURE | FFTW_MPI_TRANSPOSED_OUT); + planifft = fftw_mpi_plan_many_dft_c2r(rank, n, howmany, iblock, oblock, reinterpret_cast(fft_workspace), v_r, communicator, FFTW_MEASURE | FFTW_MPI_TRANSPOSED_IN); } // TODO: possibly circumvent the padding problem by accessing r as a matrix? @@ -390,32 +392,42 @@ void Solver::iterateCubes(F f) } template -void Solver::convolution() +void Solver::apply_preconditioner(double *input, double *output) { // it is important that at least one of the dimensions n_x and n_z is divisible by two (or local_n1, but that can't be guaranteed from the outside) // discussion of real times complex: https://forum.kde.org/viewtopic.php?f=74&t=85678 clock_t dtime = clock(); - fftw_execute(planfft); + fftw_mpi_execute_dft_r2c(planfft, input, reinterpret_cast(fft_workspace)); fft_time += clock() - dtime; buftime = clock() - dtime; + Map spectrum(reinterpret_cast *>(fft_workspace), local_n1 * n_x * (n_z / 2 + 1) * howmany); Matrix, howmany, howmany> tmp; for (ptrdiff_t i = 0; i < (local_n1 * n_x * (n_z / 2 + 1)) / 2; i++) { tmp = fundamentalSolution.template middleCols(i * (howmany + 1)).template cast>(); - rhat.segment(2 * i * howmany) = tmp.template selfadjointView() * rhat.segment(2 * i * howmany); + spectrum.segment(2 * i * howmany) = tmp.template selfadjointView() * spectrum.segment(2 * i * howmany); tmp = fundamentalSolution.template middleCols(i * (howmany + 1) + 1).template cast>(); - rhat.segment((2 * i + 1) * howmany) = tmp.template selfadjointView() * rhat.segment((2 * i + 1) * howmany); + spectrum.segment((2 * i + 1) * howmany) = tmp.template selfadjointView() * spectrum.segment((2 * i + 1) * howmany); } dtime = clock(); - fftw_execute(planifft); + fftw_mpi_execute_dft_c2r(planifft, reinterpret_cast(fft_workspace), output); fft_time += clock() - dtime; buftime += clock() - dtime; } +template +double Solver::dotProduct(RealArray &a, RealArray &b) +{ + double local_value = (a * b).sum(); + double result; + MPI_Allreduce(&local_value, &result, 1, MPI_DOUBLE, MPI_SUM, communicator); + return result; +} + template double Solver::compute_error(RealArray &r, const std::string &details) { @@ -700,7 +712,8 @@ void Solver::postprocess(Reader &reader, int load_idx, int time_ << std::setprecision(12) << homogenized_tangent << '\n'; Log::logger().info("{}", output.str()); } - reader.writeData("homogenized_tangent", load_idx, time_idx, homogenized_tangent.data(), dims, 2); + const Matrix tangent_output = homogenized_tangent; + reader.writeData("homogenized_tangent", load_idx, time_idx, tangent_output.data(), dims, 2); } } @@ -741,6 +754,9 @@ VectorXd Solver::get_homogenized_stress() template MatrixXd Solver::get_homogenized_tangent(double pert_param) { + if (matmanager->models.size() == 1 && matmanager->models[0]->has_consistent_tangent()) + return compute_consistent_homogenized_tangent(*this); + homogenized_tangent = MatrixXd::Zero(n_str, n_str); VectorXd unperturbed_stress = get_homogenized_stress(); VectorXd perturbed_stress; @@ -795,6 +811,10 @@ Solver::~Solver() fftw_free(buffer_padding); buffer_padding = nullptr; } + if (fft_workspace) { + fftw_free(fft_workspace); + fft_workspace = nullptr; + } if (planfft) { fftw_destroy_plan(planfft); } diff --git a/include/solverCG.h b/include/solverCG.h index 684a778..13a5499 100644 --- a/include/solverCG.h +++ b/include/solverCG.h @@ -7,11 +7,9 @@ template class SolverCG : public Solver { public: - using Solver::n_x; using Solver::n_y; using Solver::n_z; using Solver::local_n0; - using Solver::local_n1; using Solver::v_u_real; using Solver::v_r_real; @@ -29,7 +27,6 @@ class SolverCG : public Solver { void internalSolve() override; std::string LineSearchSecant(); - double dotProduct(RealArray &a, RealArray &b); protected: using Solver::iter; @@ -47,18 +44,7 @@ SolverCG::SolverCG(Reader &reader, MaterialManager(n_z * howmany)) -{ - this->CreateFFTWPlans(this->v_r, (fftw_complex *) s, s); -} - -template -double SolverCG::dotProduct(RealArray &a, RealArray &b) -{ - double local_value = (a * b).sum(); - double result; - MPI_Allreduce(&local_value, &result, 1, MPI_DOUBLE, MPI_SUM, this->communicator); - return result; -} +{} template void SolverCG::internalSolve() @@ -84,13 +70,13 @@ void SolverCG::internalSolve() while ((iter < this->n_it) && (err_rel > this->TOL)) { - deltamid = dotProduct(v_r_real, s_real); + deltamid = this->dotProduct(v_r_real, s_real); - this->convolution(); + this->apply_preconditioner(this->v_r, s); s_real *= -1; delta0 = delta; - delta = dotProduct(v_r_real, s_real); + delta = this->dotProduct(v_r_real, s_real); std::string details; if (islinear && !this->isMixedBCActive()) { @@ -111,7 +97,7 @@ void SolverCG::internalSolve() rnew_real -= v_r_real; } - double alpha = delta / dotProduct(d_real, rnew_real); + double alpha = delta / this->dotProduct(d_real, rnew_real); v_r_real -= alpha * rnew_real; v_u_real -= alpha * d_real; } else { @@ -137,7 +123,7 @@ std::string SolverCG::LineSearchSecant() double alpha_prev = 0.0; double alpha_curr = alpha_warm; - double rpd = dotProduct(v_r_real, d_real); + double rpd = this->dotProduct(v_r_real, d_real); if (rpd >= 0.0) { ls_converged = false; alpha_warm = 0.1; @@ -146,7 +132,7 @@ std::string SolverCG::LineSearchSecant() v_u_real += d_real * alpha_curr; this->updateMixedBC(); this->template compute_residual<0>(rnew_real, v_u_real); - double r1pd = dotProduct(rnew_real, d_real); + double r1pd = this->dotProduct(rnew_real, d_real); const double r1pd0 = this->isMixedBCActive() ? r1pd : rpd; double denom, alpha_next; @@ -167,7 +153,7 @@ std::string SolverCG::LineSearchSecant() _iter++; this->template compute_residual<0>(rnew_real, v_u_real); - r1pd = dotProduct(rnew_real, d_real); + r1pd = this->dotProduct(rnew_real, d_real); } ls_converged = (fabs(r1pd) <= tol * fabs(r1pd0)); if (ls_converged || fabs(r1pd) < fabs(r1pd0)) { diff --git a/include/solverFP.h b/include/solverFP.h index ef3124c..545564b 100644 --- a/include/solverFP.h +++ b/include/solverFP.h @@ -25,9 +25,7 @@ class SolverFP : public Solver { template SolverFP::SolverFP(Reader &reader, MaterialManager *matmanager) : Solver(reader, matmanager) -{ - this->CreateFFTWPlans(this->v_r, (fftw_complex *) this->v_r, this->v_r); -} +{} template void SolverFP::internalSolve() @@ -41,7 +39,7 @@ void SolverFP::internalSolve() while ((iter < this->n_it) && (err_rel > this->TOL)) { - this->convolution(); + this->apply_preconditioner(this->v_r, this->v_r); v_u_real -= v_r_real; this->updateMixedBC(); this->template compute_residual<2>(v_r_real, v_u_real); diff --git a/src/neml2_plugin.cpp b/src/neml2_plugin.cpp index 5ef6482..6a6f27b 100644 --- a/src/neml2_plugin.cpp +++ b/src/neml2_plugin.cpp @@ -82,6 +82,31 @@ struct FANSPluginModel { { out.copy_(reorder(v.cpu(), var).reshape(out.sizes())); } + + at::Tensor view(const double *data, int64_t n, int64_t size) const + { + return at::from_blob(const_cast(data), {n, size}, kF64); + } + + std::map inputs(int64_t n, double t_old, double t, const double *gradient, + const int *voxel, const int *phase, const double *history_old, + const double *history_current) const + { + std::map in; + in[this->gradient.name] = to_neml2(view(gradient, n, this->gradient.size), this->gradient); + in["t"] = at::full({n}, t, kF64.device(device)); + in["t~1"] = at::full({n}, t_old, kF64.device(device)); + for (const Var &f : fields) { + const auto rows = at::from_blob(const_cast(f.per_phase ? phase : voxel), {n}, at::kInt).to(device, at::kLong); + in[f.name] = f.table.index_select(0, rows); + } + const auto old = view(history_old, n, n_history), current = view(history_current, n, n_history); + for (const Var &h : history) { + in[h.name + "~1"] = to_neml2(old.narrow(1, h.offset, h.size), h); + in[h.name] = to_neml2(current.narrow(1, h.offset, h.size), h); + } + return in; + } }; extern "C" { @@ -179,28 +204,12 @@ int fans_plugin_evaluate(FANSPluginModel *m, size_t n_points, double t_old, doub double *history_new, char *err, size_t errlen) { try { - // A FANS array as an [n][size] tensor, without copying - const int64_t n = int64_t(n_points); - auto view = [n](const double *p, int64_t size) { return at::from_blob(const_cast(p), {n, size}, kF64); }; - - std::map inputs; - inputs[m->gradient.name] = m->to_neml2(view(gradient, m->gradient.size), m->gradient); - inputs["t"] = at::full({n}, t, kF64.device(m->device)); - inputs["t~1"] = at::full({n}, t_old, kF64.device(m->device)); - for (const Var &f : m->fields) { // each point's row of the table - const auto rows = at::from_blob(const_cast(f.per_phase ? phase : voxel), {n}, at::kInt).to(m->device, at::kLong); - inputs[f.name] = f.table.index_select(0, rows); - } - - const at::Tensor old_history = view(history_old, m->n_history), new_history = view(history_new, m->n_history); - for (const Var &h : m->history) { - inputs[h.name + "~1"] = m->to_neml2(old_history.narrow(1, h.offset, h.size), h); - inputs[h.name] = m->to_neml2(new_history.narrow(1, h.offset, h.size), h); // initial guess - } + const int64_t n = int64_t(n_points); + const auto new_history = m->view(history_new, n, m->n_history); + const auto outputs = m->model->forward(m->inputs(n, t_old, t, gradient, voxel, phase, + history_old, history_new)); - const auto outputs = m->model->forward(inputs); - - m->to_fans(outputs.at(m->flux.name), m->flux, view(flux, m->flux.size)); + m->to_fans(outputs.at(m->flux.name), m->flux, m->view(flux, n, m->flux.size)); for (const Var &h : m->history) m->to_fans(outputs.at(h.name), h, new_history.narrow(1, h.offset, h.size)); return 0; @@ -210,6 +219,27 @@ int fans_plugin_evaluate(FANSPluginModel *m, size_t n_points, double t_old, doub } } +int fans_plugin_tangent(FANSPluginModel *m, size_t n_points, double t_old, double t, + const double *gradient, const int *voxel, const int *phase, + const double *history_old, const double *history_current, + double *tangent, char *err, size_t errlen) +{ + try { + const int64_t n = int64_t(n_points); + Var tangent_var = m->flux; + tangent_var.shape.insert(tangent_var.shape.end(), m->gradient.shape.begin() + 1, m->gradient.shape.end()); + tangent_var.size *= m->gradient.size; + const auto result = m->model->jacobian(m->inputs(n, t_old, t, gradient, voxel, phase, + history_old, history_current)); + m->to_fans(result.second.at(m->flux.name).at(m->gradient.name), tangent_var, + m->view(tangent, n, tangent_var.size)); + return 0; + } catch (const std::exception &e) { + std::snprintf(err, errlen, "%s", e.what()); + return 1; + } +} + void fans_plugin_free(FANSPluginModel *m) { delete m; diff --git a/src/plugin_material.cpp b/src/plugin_material.cpp index 18e49fc..6958e29 100644 --- a/src/plugin_material.cpp +++ b/src/plugin_material.cpp @@ -6,10 +6,8 @@ // Batching: before each sweep over the elements, evaluate_batch gets the flux // (the stress) at all the model's Gauss points from the plugin in large // batches (fast on a GPU); the sweep then takes these cached values. -// History: the plugin's internal variables at every Gauss point: `converged`, -// the state a step starts from, and `trial`, the state of the step itself. -// Once a step has converged, the solver evaluates it a last time (for the -// output) and `trial` becomes `converged`. +// History: `previous`, `converged`, and `trial` hold the states before, at, +// and during the current step. A converged step advances the first two. // Fields: further model inputs, e.g. an orientation, read from the // microstructure file per voxel or per phase (read_field). // Linear: "linear": true in the material properties says that the flux is @@ -34,6 +32,8 @@ struct Plugin { decltype(&fans_plugin_set_table) set_table; decltype(&fans_plugin_evaluate) evaluate; decltype(&fans_plugin_free) free_model; + + decltype(&fans_plugin_tangent) tangent; }; // "NEML2" -> libfans_neml2.so, found next to libFANS.so through its RPATH @@ -57,6 +57,7 @@ Plugin load_plugin(const string &name) sym(p.load, "fans_plugin_load"); sym(p.set_table, "fans_plugin_set_table"); sym(p.evaluate, "fans_plugin_evaluate"); + p.tangent = reinterpret_cast(dlsym(h, "fans_plugin_tangent")); sym(p.free_model, "fans_plugin_free"); return p; } @@ -156,7 +157,8 @@ class PluginModel : public Base { // The history of every Gauss point, and the arrays of one plugin call void initializeInternalVariables(ptrdiff_t num_elements, int num_gauss_points) override { - trial = converged = vector(num_elements * num_gauss_points * n_history, 0.0); + previous = trial = converged = vector(num_elements * num_gauss_points * n_history, 0.0); + gradient_field.resize(num_elements * num_gauss_points * n_str); gradient_buf.resize(chunk * n_gp * n_str); flux_buf.resize(chunk * n_gp * n_str); history_old_buf.resize(chunk * n_gp * n_history); @@ -168,6 +170,7 @@ class PluginModel : public Base { // A step has converged: its state is what the next step starts from void updateInternalVariables() override { + previous.swap(converged); converged = trial; } @@ -199,6 +202,8 @@ class PluginModel : public Base { return true; } + bool has_consistent_tangent() const override { return plugin.tangent != nullptr; } + // "linear": true in the material properties; FANS trusts it bool is_linear() const override { @@ -237,11 +242,35 @@ class PluginModel : public Base { for (size_t k = 0; k < ne; ++k) { const ptrdiff_t e = elems[b + k]; std::copy_n(&flux_buf[k * per_elem], per_elem, flux_all + e * per_elem); + std::copy_n(&gradient_buf[k * per_elem], per_elem, gradient_field.data() + e * per_elem); std::copy_n(history_new_buf.data() + k * hist_elem, hist_elem, trial.data() + e * hist_elem); } } } + void compute_tangent_field(ptrdiff_t n_elem, const unsigned short *phase, vector &tangent) override + { + if (!plugin.tangent) + throw std::runtime_error("The material plugin does not provide an algorithmic tangent."); + const size_t per_elem = size_t(n_gp) * n_str, tangent_elem = per_elem * n_str; + const size_t hist_elem = size_t(n_gp) * n_history; + tangent.resize(size_t(n_elem) * tangent_elem); + for (ptrdiff_t b = 0; b < n_elem; b += chunk) { + const size_t ne = std::min(chunk, n_elem - b); + for (size_t k = 0; k < ne; ++k) { + const ptrdiff_t e = b + k; + std::fill_n(&voxel_buf[k * n_gp], n_gp, int(e)); + std::fill_n(&phase_buf[k * n_gp], n_gp, phase[e]); + } + char err[FANS_PLUGIN_MSGLEN] = {0}; + if (plugin.tangent(model, ne * n_gp, this->time_old, this->time, gradient_field.data() + b * per_elem, + voxel_buf.data(), phase_buf.data(), previous.data() + b * hist_elem, + converged.data() + b * hist_elem, tangent.data() + b * tangent_elem, + err, sizeof(err)) != 0) + throw std::runtime_error(string("Plugin tangent evaluation failed: ") + err); + } + } + // Writes every history variable, e.g. state/internal/Ep as state_internal_Ep void postprocess(Solver &solver, Reader &reader, int load_idx, int time_idx) override { @@ -289,6 +318,7 @@ class PluginModel : public Base { size_t n_phase_rows{SIZE_MAX}; // phases below have a row in every per-phase field vector trial, converged; // history at every Gauss point vector gradient_buf, flux_buf, history_old_buf, history_new_buf; // one plugin call + vector previous, gradient_field; vector voxel_buf, phase_buf; const double *flux_cache{nullptr}; }; diff --git a/src/tangent.cpp b/src/tangent.cpp new file mode 100644 index 0000000..44d5e6e --- /dev/null +++ b/src/tangent.cpp @@ -0,0 +1,87 @@ +#include "solver.h" + +template +using MacroVector = Matrix; + +template +using ElementVector = Matrix; + +template +static void tangent_residual(Solver &solver, RealArray &residual, RealArray &delta_u, + const double *tangent, const MacroVector ¯o) +{ + auto *model = solver.matmanager->models[0]; + solver.template compute_residual_basic( + residual, delta_u, [&](ElementVector &ue, int, ptrdiff_t element) -> ElementVector & { + return model->element_tangent_residual(ue, macro, tangent + element * model->n_gp * n_str * n_str); + }); +} + +template +static MacroVector tangent_average(Solver &solver, RealArray &delta_u, RealArray &scratch, + const double *tangent, const MacroVector ¯o) +{ + auto *model = solver.matmanager->models[0]; + MacroVector average = MacroVector::Zero(); + ElementVector zero = ElementVector::Zero(); + solver.template compute_residual_basic<0>(scratch, delta_u, [&](ElementVector &ue, int, ptrdiff_t element) -> ElementVector & { + average += model->element_tangent_stress(ue, macro, tangent + element * model->n_gp * n_str * n_str); + return zero; + }); + MPI_Allreduce(MPI_IN_PLACE, average.data(), n_str, MPI_DOUBLE, MPI_SUM, solver.communicator); + return average / (solver.n_x * solver.n_y * solver.n_z); +} + +template +static void solve_tangent(Solver &solver, RealArray &delta_u, RealArray &r, RealArray &d, + RealArray &Kd, RealArray &z, const double *tangent, const MacroVector ¯o) +{ + tangent_residual(solver, r, delta_u, tangent, macro); + d.setZero(); + double rho = 0.0, error = std::sqrt(solver.dotProduct(r, r)); + const double limit = solver.reader.errorParameters["type"] == "relative" ? solver.TOL * error : solver.TOL; + for (int i = 0; i < solver.n_it && error > limit; ++i) { + solver.apply_preconditioner(r.data(), z.data()); + z *= -1.0; + const double old_rho = rho; + rho = solver.dotProduct(r, z); + d = z + (i ? rho / old_rho : 0.0) * d; + tangent_residual(solver, Kd, d, tangent, MacroVector::Zero()); + const double alpha = rho / solver.dotProduct(d, Kd); + delta_u -= alpha * d; + r -= alpha * Kd; + error = std::sqrt(solver.dotProduct(r, r)); + } +} + +template +MatrixXd compute_consistent_homogenized_tangent(Solver &solver) +{ + if (solver.isMixedBCActive()) + throw std::runtime_error("Homogenized tangent is not supported with mixed boundary conditions."); + const auto n_elem = solver.local_n0 * solver.n_y * solver.n_z; + const size_t n_work = (solver.local_n0 + 1) * solver.n_y * solver.n_z * howmany; + const size_t n_fft = std::max(solver.reader.alloc_local * 2, (solver.local_n0 + 1) * solver.n_y * (solver.n_z + 2) * howmany); + vector displacement(n_work), search_direction(n_work); + unique_ptr residual(fftw_alloc_real(n_fft), fftw_free), product(fftw_alloc_real(n_fft), fftw_free); + vector local_tangent; + auto *model = solver.matmanager->models[0]; + model->compute_tangent_field(n_elem, solver.ms, local_tangent); + RealArray delta_u(displacement.data(), solver.n_z * howmany, solver.local_n0 * solver.n_y, OuterStride<>(solver.n_z * howmany)); + RealArray r(residual.get(), solver.n_z * howmany, solver.local_n0 * solver.n_y, OuterStride<>((solver.n_z + 2) * howmany)); + RealArray d(search_direction.data(), solver.n_z * howmany, solver.local_n0 * solver.n_y, OuterStride<>(solver.n_z * howmany)); + RealArray Kd(product.get(), solver.n_z * howmany, solver.local_n0 * solver.n_y, OuterStride<>(solver.n_z * howmany)); + RealArray z(product.get(), solver.n_z * howmany, solver.local_n0 * solver.n_y, OuterStride<>((solver.n_z + 2) * howmany)); + MatrixXd homogenized_tangent(n_str, n_str); + for (int j = 0; j < n_str; ++j) { + delta_u.setZero(); + const MacroVector direction = MacroVector::Unit(j); + solve_tangent(solver, delta_u, r, d, Kd, z, local_tangent.data(), direction); + homogenized_tangent.col(j) = tangent_average(solver, delta_u, Kd, local_tangent.data(), direction); + } + return homogenized_tangent; +} + +template MatrixXd compute_consistent_homogenized_tangent(Solver<1, 3> &); +template MatrixXd compute_consistent_homogenized_tangent(Solver<3, 6> &); +template MatrixXd compute_consistent_homogenized_tangent(Solver<3, 9> &); From 371e8f04a7feb5ff35bc86c65e6384b215fec717 Mon Sep 17 00:00:00 2001 From: Sanath Date: Tue, 6 Oct 2026 18:59:59 +0200 Subject: [PATCH 2/8] simplified --- docs/NEML2.md | 6 +- include/MaterialManager.h | 7 ++ include/fans_plugin.h | 11 +-- include/matmodel.h | 30 ++------- include/solver.h | 59 ++++++---------- include/solverCG.h | 30 ++++++--- include/solverFP.h | 6 +- pixi.toml | 4 -- src/neml2_plugin.cpp | 83 +++++++++-------------- src/plugin_material.cpp | 58 +++++++++------- src/tangent.cpp | 138 ++++++++++++++++++-------------------- 11 files changed, 194 insertions(+), 238 deletions(-) diff --git a/docs/NEML2.md b/docs/NEML2.md index 390d648..c5ead69 100644 --- a/docs/NEML2.md +++ b/docs/NEML2.md @@ -27,12 +27,10 @@ This gives `test/FANS_neml2`. Without pixi, configure with `-DFANS_NEML2=ON` in ## Compile a model ```bash -pixi run -e dev-neml2 neml2-compile model.i "cpu cuda" +pixi run -e dev-neml2 neml2-compile model.i --model model --device cpu cuda --dtype float64 --output-dir compiled_models -d state/S:forces/E ``` -This compiles the model block named like the file (`[model]`) for both devices into `compiled_models/model`, next to `model.i`. A third argument lists Python files the model needs, e.g. `"a.py b.py"`. - -For `homogenized_tangent`, compile the model's flux-gradient derivative, e.g. `neml2-compile ... -d state/S:forces/E` for small strain. The current sensitivity solve uses PCG and therefore requires a symmetric, CG-compatible algorithmic tangent. +This compiles the block `[model]` of `model.i` for both devices into `compiled_models/model`. `--load a.py` (repeatable) imports a Python file the model needs. `-d flux:gradient` also compiles the model's tangent. ## Input file diff --git a/include/MaterialManager.h b/include/MaterialManager.h index 0c4ecfe..855de88 100644 --- a/include/MaterialManager.h +++ b/include/MaterialManager.h @@ -138,6 +138,13 @@ class MaterialManager { } } + // One linear model for all phases, e.g. the frozen tangent of src/tangent.cpp + MaterialManager(Matmodel *model, int n_phases) + : phase_to_info(new MaterialInfo[n_phases]), n_phases(n_phases), models{model}, all_stiffness(false) + { + std::fill_n(phase_to_info, n_phases, MaterialInfo{model, nullptr, 0, true}); + } + ~MaterialManager() { for (auto *m : models) diff --git a/include/fans_plugin.h b/include/fans_plugin.h index 5ea9fdd..5463f49 100644 --- a/include/fans_plugin.h +++ b/include/fans_plugin.h @@ -53,16 +53,11 @@ FANSPluginModel *fans_plugin_load(const char *config, char *msg, size_t msglen); int fans_plugin_set_table(FANSPluginModel *model, size_t field, int per_phase, size_t n_rows, const double *table, char *err, size_t errlen); -/* The step runs from time t_old to t. */ +/* The step runs from time t_old to t. Unless NULL, `tangent` gets the consistent + tangent d flux / d gradient, [n_points][n_str][n_str]. */ int fans_plugin_evaluate(FANSPluginModel *model, size_t n_points, double t_old, double t, const double *gradient, const int *voxel, const int *phase, const double *history_old, double *flux, - double *history_new, char *err, size_t errlen); - -/* Consistent flux-gradient tangent [n_points][n_str][n_str]. */ -int fans_plugin_tangent(FANSPluginModel *model, size_t n_points, double t_old, double t, - const double *gradient, const int *voxel, const int *phase, - const double *history_old, const double *history_current, - double *tangent, char *err, size_t errlen); + double *history_new, double *tangent, char *err, size_t errlen); void fans_plugin_free(FANSPluginModel *model); diff --git a/include/matmodel.h b/include/matmodel.h index 3222753..7d5e32d 100644 --- a/include/matmodel.h +++ b/include/matmodel.h @@ -54,32 +54,14 @@ class Matmodel { } virtual void evaluate_batch(const vector &elems, const unsigned short *phase, const double *ue, double *sig_gp) {} - virtual bool has_consistent_tangent() const { return false; } - virtual void compute_tangent_field(ptrdiff_t, const unsigned short *, vector &) + // The tangent d flux / d gradient of the converged step at the Gauss points of `elems`, + // [element][n_gp][n_str][n_str]; false if the model has none. + virtual bool evaluate_tangent(const vector &elems, const unsigned short *phase, + const double *ue, double *tangent_gp) { - throw std::logic_error("Material tangent not implemented."); - } - Matrix &element_tangent_residual( - const Matrix &ue, const Matrix ¯o, const double *tangent) - { - res_e.setZero(); - for (int gp = 0; gp < n_gp; ++gp) { - const Matrix deps = B_int[gp] * ue + macro; - const Matrix dsig = Map>(tangent + gp * n_str * n_str) * deps; - res_e.noalias() += B_int[gp].transpose() * dsig * v_e / n_gp; - } - return res_e; - } - Matrix element_tangent_stress( - const Matrix &ue, const Matrix ¯o, const double *tangent) const - { - Matrix stress = Matrix::Zero(); - for (int gp = 0; gp < n_gp; ++gp) { - const Matrix deps = B_int[gp] * ue + macro; - stress.noalias() += Map>(tangent + gp * n_str * n_str) * deps / n_gp; - } - return stress; + return false; } + // A model without the element stiffnesses of a LinearModel can still say that its // flux is linear in the gradient: then the linear CG applies. virtual bool is_linear() const diff --git a/include/solver.h b/include/solver.h index dadfb8c..77a601b 100644 --- a/include/solver.h +++ b/include/solver.h @@ -40,6 +40,7 @@ class Solver : private MixedBCController { RealArray v_r_real; // can't do the "classname()" intialization here, and Map doesn't have a default constructor RealArray v_u_real; + Map rhat; ArrayXd err_all; //!< Absolute error history Matrix fundamentalSolution; @@ -62,9 +63,9 @@ class Solver : private MixedBCController { void postprocess(Reader &reader, int load_idx, int time_idx); //!< Computes Strain and stress - void apply_preconditioner(double *input, double *output); - double dotProduct(RealArray &a, RealArray &b); + void convolution(); double compute_error(RealArray &r, const std::string &details = {}); + void CreateFFTWPlans(double *in, fftw_complex *transformed, double *out); VectorXd homogenized_strain; VectorXd homogenized_stress; @@ -73,6 +74,7 @@ class Solver : private MixedBCController { MatrixXd homogenized_tangent; MatrixXd get_homogenized_tangent(double pert_param); + bool consistent_homogenized_tangent(); //!< in src/tangent.cpp void enableMixedBC(const MixedBC &mbc, size_t step) { @@ -100,19 +102,14 @@ class Solver : private MixedBCController { (svd.singularValues().array() > tolerance).select(svd.singularValues().array().inverse(), 0).matrix().asDiagonal() * svd.matrixU().adjoint(); } - double *fft_workspace; - fftw_plan planfft{nullptr}, planifft{nullptr}; - void createFFTWPlans(); void computeFundamentalSolution(); protected: + fftw_plan planfft, planifft; clock_t fft_time, buftime; size_t iter; }; -template -MatrixXd compute_consistent_homogenized_tangent(Solver &solver); - template Solver::Solver(Reader &reader, MaterialManager *matmgr) : reader(reader), @@ -139,8 +136,8 @@ Solver::Solver(Reader &reader, MaterialManager * v_u_real(v_u, n_z * howmany, local_n0 * n_y, OuterStride<>(n_z * howmany)), v_u_prev(fftw_alloc_real(local_n0 * n_y * n_z * howmany)), - buffer_padding(fftw_alloc_real(n_y * (n_z + 2) * howmany)), - fft_workspace(fftw_alloc_real(reader.alloc_local * 2)) + rhat((std::complex *) v_r, local_n1 * n_x * (n_z / 2 + 1) * howmany), // actual initialization is below + buffer_padding(fftw_alloc_real(n_y * (n_z + 2) * howmany)) { v_u_real.setZero(); for (ptrdiff_t i = local_n0 * n_y * n_z * howmany; i < (local_n0 + 1) * n_y * n_z * howmany; i++) { @@ -151,7 +148,6 @@ Solver::Solver(Reader &reader, MaterialManager * matmanager->initialize_internal_variables(local_n0 * n_y * n_z, matmanager->models[0]->n_gp); computeFundamentalSolution(); - createFFTWPlans(); } template @@ -213,7 +209,7 @@ void Solver::computeFundamentalSolution() } template -void Solver::createFFTWPlans() +void Solver::CreateFFTWPlans(double *in, fftw_complex *transformed, double *out) { int rank = 3; ptrdiff_t iblock = FFTW_MPI_DEFAULT_BLOCK; @@ -227,8 +223,11 @@ void Solver::createFFTWPlans() // But, according to https://fftw.org/doc/MPI-Plan-Creation.html the BLOCK sizes must be the same: // "These must be the same block sizes as were passed to the corresponding ‘local_size’ function" const ptrdiff_t n[3] = {n_x, n_y, n_z}; - planfft = fftw_mpi_plan_many_dft_r2c(rank, n, howmany, iblock, oblock, v_r, reinterpret_cast(fft_workspace), communicator, FFTW_MEASURE | FFTW_MPI_TRANSPOSED_OUT); - planifft = fftw_mpi_plan_many_dft_c2r(rank, n, howmany, iblock, oblock, reinterpret_cast(fft_workspace), v_r, communicator, FFTW_MEASURE | FFTW_MPI_TRANSPOSED_IN); + planfft = fftw_mpi_plan_many_dft_r2c(rank, n, howmany, iblock, oblock, in, transformed, communicator, FFTW_MEASURE | FFTW_MPI_TRANSPOSED_OUT); + planifft = fftw_mpi_plan_many_dft_c2r(rank, n, howmany, iblock, oblock, transformed, out, communicator, FFTW_MEASURE | FFTW_MPI_TRANSPOSED_IN); + + // see https://eigen.tuxfamily.org/dox/group__TutorialMapClass.html#title3 + new (&rhat) Map((std::complex *) transformed, local_n1 * n_x * (n_z / 2 + 1) * howmany); } // TODO: possibly circumvent the padding problem by accessing r as a matrix? @@ -392,42 +391,32 @@ void Solver::iterateCubes(F f) } template -void Solver::apply_preconditioner(double *input, double *output) +void Solver::convolution() { // it is important that at least one of the dimensions n_x and n_z is divisible by two (or local_n1, but that can't be guaranteed from the outside) // discussion of real times complex: https://forum.kde.org/viewtopic.php?f=74&t=85678 clock_t dtime = clock(); - fftw_mpi_execute_dft_r2c(planfft, input, reinterpret_cast(fft_workspace)); + fftw_execute(planfft); fft_time += clock() - dtime; buftime = clock() - dtime; - Map spectrum(reinterpret_cast *>(fft_workspace), local_n1 * n_x * (n_z / 2 + 1) * howmany); Matrix, howmany, howmany> tmp; for (ptrdiff_t i = 0; i < (local_n1 * n_x * (n_z / 2 + 1)) / 2; i++) { tmp = fundamentalSolution.template middleCols(i * (howmany + 1)).template cast>(); - spectrum.segment(2 * i * howmany) = tmp.template selfadjointView() * spectrum.segment(2 * i * howmany); + rhat.segment(2 * i * howmany) = tmp.template selfadjointView() * rhat.segment(2 * i * howmany); tmp = fundamentalSolution.template middleCols(i * (howmany + 1) + 1).template cast>(); - spectrum.segment((2 * i + 1) * howmany) = tmp.template selfadjointView() * spectrum.segment((2 * i + 1) * howmany); + rhat.segment((2 * i + 1) * howmany) = tmp.template selfadjointView() * rhat.segment((2 * i + 1) * howmany); } dtime = clock(); - fftw_mpi_execute_dft_c2r(planifft, reinterpret_cast(fft_workspace), output); + fftw_execute(planifft); fft_time += clock() - dtime; buftime += clock() - dtime; } -template -double Solver::dotProduct(RealArray &a, RealArray &b) -{ - double local_value = (a * b).sum(); - double result; - MPI_Allreduce(&local_value, &result, 1, MPI_DOUBLE, MPI_SUM, communicator); - return result; -} - template double Solver::compute_error(RealArray &r, const std::string &details) { @@ -712,8 +701,8 @@ void Solver::postprocess(Reader &reader, int load_idx, int time_ << std::setprecision(12) << homogenized_tangent << '\n'; Log::logger().info("{}", output.str()); } - const Matrix tangent_output = homogenized_tangent; - reader.writeData("homogenized_tangent", load_idx, time_idx, tangent_output.data(), dims, 2); + const Matrix rows = homogenized_tangent; // as the file stores it + reader.writeData("homogenized_tangent", load_idx, time_idx, rows.data(), dims, 2); } } @@ -754,8 +743,8 @@ VectorXd Solver::get_homogenized_stress() template MatrixXd Solver::get_homogenized_tangent(double pert_param) { - if (matmanager->models.size() == 1 && matmanager->models[0]->has_consistent_tangent()) - return compute_consistent_homogenized_tangent(*this); + if (consistent_homogenized_tangent()) // from the materials' own tangents, if they all have one + return homogenized_tangent; homogenized_tangent = MatrixXd::Zero(n_str, n_str); VectorXd unperturbed_stress = get_homogenized_stress(); @@ -811,10 +800,6 @@ Solver::~Solver() fftw_free(buffer_padding); buffer_padding = nullptr; } - if (fft_workspace) { - fftw_free(fft_workspace); - fft_workspace = nullptr; - } if (planfft) { fftw_destroy_plan(planfft); } diff --git a/include/solverCG.h b/include/solverCG.h index 13a5499..684a778 100644 --- a/include/solverCG.h +++ b/include/solverCG.h @@ -7,9 +7,11 @@ template class SolverCG : public Solver { public: + using Solver::n_x; using Solver::n_y; using Solver::n_z; using Solver::local_n0; + using Solver::local_n1; using Solver::v_u_real; using Solver::v_r_real; @@ -27,6 +29,7 @@ class SolverCG : public Solver { void internalSolve() override; std::string LineSearchSecant(); + double dotProduct(RealArray &a, RealArray &b); protected: using Solver::iter; @@ -44,7 +47,18 @@ SolverCG::SolverCG(Reader &reader, MaterialManager(n_z * howmany)) -{} +{ + this->CreateFFTWPlans(this->v_r, (fftw_complex *) s, s); +} + +template +double SolverCG::dotProduct(RealArray &a, RealArray &b) +{ + double local_value = (a * b).sum(); + double result; + MPI_Allreduce(&local_value, &result, 1, MPI_DOUBLE, MPI_SUM, this->communicator); + return result; +} template void SolverCG::internalSolve() @@ -70,13 +84,13 @@ void SolverCG::internalSolve() while ((iter < this->n_it) && (err_rel > this->TOL)) { - deltamid = this->dotProduct(v_r_real, s_real); + deltamid = dotProduct(v_r_real, s_real); - this->apply_preconditioner(this->v_r, s); + this->convolution(); s_real *= -1; delta0 = delta; - delta = this->dotProduct(v_r_real, s_real); + delta = dotProduct(v_r_real, s_real); std::string details; if (islinear && !this->isMixedBCActive()) { @@ -97,7 +111,7 @@ void SolverCG::internalSolve() rnew_real -= v_r_real; } - double alpha = delta / this->dotProduct(d_real, rnew_real); + double alpha = delta / dotProduct(d_real, rnew_real); v_r_real -= alpha * rnew_real; v_u_real -= alpha * d_real; } else { @@ -123,7 +137,7 @@ std::string SolverCG::LineSearchSecant() double alpha_prev = 0.0; double alpha_curr = alpha_warm; - double rpd = this->dotProduct(v_r_real, d_real); + double rpd = dotProduct(v_r_real, d_real); if (rpd >= 0.0) { ls_converged = false; alpha_warm = 0.1; @@ -132,7 +146,7 @@ std::string SolverCG::LineSearchSecant() v_u_real += d_real * alpha_curr; this->updateMixedBC(); this->template compute_residual<0>(rnew_real, v_u_real); - double r1pd = this->dotProduct(rnew_real, d_real); + double r1pd = dotProduct(rnew_real, d_real); const double r1pd0 = this->isMixedBCActive() ? r1pd : rpd; double denom, alpha_next; @@ -153,7 +167,7 @@ std::string SolverCG::LineSearchSecant() _iter++; this->template compute_residual<0>(rnew_real, v_u_real); - r1pd = this->dotProduct(rnew_real, d_real); + r1pd = dotProduct(rnew_real, d_real); } ls_converged = (fabs(r1pd) <= tol * fabs(r1pd0)); if (ls_converged || fabs(r1pd) < fabs(r1pd0)) { diff --git a/include/solverFP.h b/include/solverFP.h index 545564b..ef3124c 100644 --- a/include/solverFP.h +++ b/include/solverFP.h @@ -25,7 +25,9 @@ class SolverFP : public Solver { template SolverFP::SolverFP(Reader &reader, MaterialManager *matmanager) : Solver(reader, matmanager) -{} +{ + this->CreateFFTWPlans(this->v_r, (fftw_complex *) this->v_r, this->v_r); +} template void SolverFP::internalSolve() @@ -39,7 +41,7 @@ void SolverFP::internalSolve() while ((iter < this->n_it) && (err_rel > this->TOL)) { - this->apply_preconditioner(this->v_r, this->v_r); + this->convolution(); v_u_real -= v_r_real; this->updateMixedBC(); this->template compute_residual<2>(v_r_real, v_u_real); diff --git a/pixi.toml b/pixi.toml index d5ca5d9..a41fd92 100644 --- a/pixi.toml +++ b/pixi.toml @@ -70,10 +70,6 @@ neml2 = ">=3.1.0, <4" [feature.neml2.tasks] build-fans_neml2 = "cmake -S . -B build-neml2 -DFANS_NEML2=ON -DFANS_ENABLE_NATIVE=ON && cmake --build build-neml2 -j" -[feature.neml2.tasks.neml2-compile] -args = ["file", { arg = "device", default = "cpu" }, { arg = "load", default = "" }] -cmd = "cd \"$INIT_CWD\" && neml2-compile {{ file }} --model $(basename {{ file }} .i) --device {{ device }} --dtype float64 --output-dir $(dirname {{ file }})/compiled_models{% for f in load | split %} --load {{ f }}{% endfor %}" - [tasks] test = { cmd = "pytest -v -s --from-pixi", cwd = "test/pytest", depends-on = ["test-fans"] } test-fans = { cmd = "./run_tests.sh -n {{n}}", args = [ diff --git a/src/neml2_plugin.cpp b/src/neml2_plugin.cpp index 6a6f27b..fbea787 100644 --- a/src/neml2_plugin.cpp +++ b/src/neml2_plugin.cpp @@ -82,31 +82,6 @@ struct FANSPluginModel { { out.copy_(reorder(v.cpu(), var).reshape(out.sizes())); } - - at::Tensor view(const double *data, int64_t n, int64_t size) const - { - return at::from_blob(const_cast(data), {n, size}, kF64); - } - - std::map inputs(int64_t n, double t_old, double t, const double *gradient, - const int *voxel, const int *phase, const double *history_old, - const double *history_current) const - { - std::map in; - in[this->gradient.name] = to_neml2(view(gradient, n, this->gradient.size), this->gradient); - in["t"] = at::full({n}, t, kF64.device(device)); - in["t~1"] = at::full({n}, t_old, kF64.device(device)); - for (const Var &f : fields) { - const auto rows = at::from_blob(const_cast(f.per_phase ? phase : voxel), {n}, at::kInt).to(device, at::kLong); - in[f.name] = f.table.index_select(0, rows); - } - const auto old = view(history_old, n, n_history), current = view(history_current, n, n_history); - for (const Var &h : history) { - in[h.name + "~1"] = to_neml2(old.narrow(1, h.offset, h.size), h); - in[h.name] = to_neml2(current.narrow(1, h.offset, h.size), h); - } - return in; - } }; extern "C" { @@ -201,15 +176,40 @@ int fans_plugin_set_table(FANSPluginModel *m, size_t field, int per_phase, size_ int fans_plugin_evaluate(FANSPluginModel *m, size_t n_points, double t_old, double t, const double *gradient, const int *voxel, const int *phase, const double *history_old, double *flux, - double *history_new, char *err, size_t errlen) + double *history_new, double *tangent, char *err, size_t errlen) { try { - const int64_t n = int64_t(n_points); - const auto new_history = m->view(history_new, n, m->n_history); - const auto outputs = m->model->forward(m->inputs(n, t_old, t, gradient, voxel, phase, - history_old, history_new)); + // A FANS array as an [n][size] tensor, without copying + const int64_t n = int64_t(n_points); + auto view = [n](const double *p, int64_t size) { return at::from_blob(const_cast(p), {n, size}, kF64); }; + + std::map inputs; + inputs[m->gradient.name] = m->to_neml2(view(gradient, m->gradient.size), m->gradient); + inputs["t"] = at::full({n}, t, kF64.device(m->device)); + inputs["t~1"] = at::full({n}, t_old, kF64.device(m->device)); + for (const Var &f : m->fields) { // each point's row of the table + const auto rows = at::from_blob(const_cast(f.per_phase ? phase : voxel), {n}, at::kInt).to(m->device, at::kLong); + inputs[f.name] = f.table.index_select(0, rows); + } + + const at::Tensor old_history = view(history_old, m->n_history), new_history = view(history_new, m->n_history); + for (const Var &h : m->history) { + inputs[h.name + "~1"] = m->to_neml2(old_history.narrow(1, h.offset, h.size), h); + inputs[h.name] = m->to_neml2(new_history.narrow(1, h.offset, h.size), h); // initial guess + } - m->to_fans(outputs.at(m->flux.name), m->flux, m->view(flux, n, m->flux.size)); + std::map outputs; + if (!tangent) { + outputs = m->model->forward(inputs); + } else { // the same evaluation with its Jacobian, whose flux-gradient block has the flux's shape, then the gradient's + auto [values, jacobian] = m->model->jacobian(inputs); + outputs = std::move(values); + Var d = m->flux; + d.shape.insert(d.shape.end(), m->gradient.shape.begin() + 1, m->gradient.shape.end()); + m->to_fans(jacobian.at(m->flux.name).at(m->gradient.name), d, view(tangent, m->flux.size * m->gradient.size)); + } + + m->to_fans(outputs.at(m->flux.name), m->flux, view(flux, m->flux.size)); for (const Var &h : m->history) m->to_fans(outputs.at(h.name), h, new_history.narrow(1, h.offset, h.size)); return 0; @@ -219,27 +219,6 @@ int fans_plugin_evaluate(FANSPluginModel *m, size_t n_points, double t_old, doub } } -int fans_plugin_tangent(FANSPluginModel *m, size_t n_points, double t_old, double t, - const double *gradient, const int *voxel, const int *phase, - const double *history_old, const double *history_current, - double *tangent, char *err, size_t errlen) -{ - try { - const int64_t n = int64_t(n_points); - Var tangent_var = m->flux; - tangent_var.shape.insert(tangent_var.shape.end(), m->gradient.shape.begin() + 1, m->gradient.shape.end()); - tangent_var.size *= m->gradient.size; - const auto result = m->model->jacobian(m->inputs(n, t_old, t, gradient, voxel, phase, - history_old, history_current)); - m->to_fans(result.second.at(m->flux.name).at(m->gradient.name), tangent_var, - m->view(tangent, n, tangent_var.size)); - return 0; - } catch (const std::exception &e) { - std::snprintf(err, errlen, "%s", e.what()); - return 1; - } -} - void fans_plugin_free(FANSPluginModel *m) { delete m; diff --git a/src/plugin_material.cpp b/src/plugin_material.cpp index 6958e29..b0e01b5 100644 --- a/src/plugin_material.cpp +++ b/src/plugin_material.cpp @@ -6,8 +6,11 @@ // Batching: before each sweep over the elements, evaluate_batch gets the flux // (the stress) at all the model's Gauss points from the plugin in large // batches (fast on a GPU); the sweep then takes these cached values. -// History: `previous`, `converged`, and `trial` hold the states before, at, -// and during the current step. A converged step advances the first two. +// History: the plugin's internal variables at every Gauss point: `converged`, +// the state a step starts from, and `trial`, the state of the step itself. +// Once a step has converged, the solver evaluates it a last time (for the +// output) and `trial` becomes `converged`; `previous` keeps the state that +// step started from, for its tangent (evaluate_tangent). // Fields: further model inputs, e.g. an orientation, read from the // microstructure file per voxel or per phase (read_field). // Linear: "linear": true in the material properties says that the flux is @@ -32,8 +35,6 @@ struct Plugin { decltype(&fans_plugin_set_table) set_table; decltype(&fans_plugin_evaluate) evaluate; decltype(&fans_plugin_free) free_model; - - decltype(&fans_plugin_tangent) tangent; }; // "NEML2" -> libfans_neml2.so, found next to libFANS.so through its RPATH @@ -57,7 +58,6 @@ Plugin load_plugin(const string &name) sym(p.load, "fans_plugin_load"); sym(p.set_table, "fans_plugin_set_table"); sym(p.evaluate, "fans_plugin_evaluate"); - p.tangent = reinterpret_cast(dlsym(h, "fans_plugin_tangent")); sym(p.free_model, "fans_plugin_free"); return p; } @@ -158,7 +158,6 @@ class PluginModel : public Base { void initializeInternalVariables(ptrdiff_t num_elements, int num_gauss_points) override { previous = trial = converged = vector(num_elements * num_gauss_points * n_history, 0.0); - gradient_field.resize(num_elements * num_gauss_points * n_str); gradient_buf.resize(chunk * n_gp * n_str); flux_buf.resize(chunk * n_gp * n_str); history_old_buf.resize(chunk * n_gp * n_history); @@ -202,8 +201,6 @@ class PluginModel : public Base { return true; } - bool has_consistent_tangent() const override { return plugin.tangent != nullptr; } - // "linear": true in the material properties; FANS trusts it bool is_linear() const override { @@ -236,39 +233,49 @@ class PluginModel : public Base { char err[FANS_PLUGIN_MSGLEN] = {0}; if (plugin.evaluate(model, ne * n_gp, this->time_old, this->time, gradient_buf.data(), voxel_buf.data(), phase_buf.data(), - history_old_buf.data(), flux_buf.data(), history_new_buf.data(), err, sizeof(err)) != 0) + history_old_buf.data(), flux_buf.data(), history_new_buf.data(), nullptr, err, sizeof(err)) != 0) throw std::runtime_error(string("Plugin material evaluation failed: ") + err); for (size_t k = 0; k < ne; ++k) { const ptrdiff_t e = elems[b + k]; std::copy_n(&flux_buf[k * per_elem], per_elem, flux_all + e * per_elem); - std::copy_n(&gradient_buf[k * per_elem], per_elem, gradient_field.data() + e * per_elem); std::copy_n(history_new_buf.data() + k * hist_elem, hist_elem, trial.data() + e * hist_elem); } } } - void compute_tangent_field(ptrdiff_t n_elem, const unsigned short *phase, vector &tangent) override + // The tangent of the converged step at all Gauss points of `elems`, into tangent_all: + // that step once more, from the state it started from (`previous`) to the one it reached + bool evaluate_tangent(const vector &elems, const unsigned short *phase, const double *ue_all, + double *tangent_all) override { - if (!plugin.tangent) - throw std::runtime_error("The material plugin does not provide an algorithmic tangent."); - const size_t per_elem = size_t(n_gp) * n_str, tangent_elem = per_elem * n_str; - const size_t hist_elem = size_t(n_gp) * n_history; - tangent.resize(size_t(n_elem) * tangent_elem); - for (ptrdiff_t b = 0; b < n_elem; b += chunk) { - const size_t ne = std::min(chunk, n_elem - b); + const size_t n_dof = size_t(howmany) * 8; + const size_t per_elem = size_t(n_gp) * n_str; + const size_t hist_elem = size_t(n_gp) * n_history; + vector tangent_buf(chunk * per_elem * n_str); + + for (size_t b = 0; b < elems.size(); b += chunk) { + const size_t ne = std::min(chunk, elems.size() - b); + for (size_t k = 0; k < ne; ++k) { - const ptrdiff_t e = b + k; + const ptrdiff_t e = elems[b + k]; + Eigen::Map(&gradient_buf[k * per_elem], per_elem).noalias() = + B * Eigen::Map(ue_all + e * n_dof, n_dof) + g0; std::fill_n(&voxel_buf[k * n_gp], n_gp, int(e)); std::fill_n(&phase_buf[k * n_gp], n_gp, phase[e]); + std::copy_n(previous.data() + e * hist_elem, hist_elem, history_old_buf.data() + k * hist_elem); + std::copy_n(converged.data() + e * hist_elem, hist_elem, history_new_buf.data() + k * hist_elem); // initial guess } + char err[FANS_PLUGIN_MSGLEN] = {0}; - if (plugin.tangent(model, ne * n_gp, this->time_old, this->time, gradient_field.data() + b * per_elem, - voxel_buf.data(), phase_buf.data(), previous.data() + b * hist_elem, - converged.data() + b * hist_elem, tangent.data() + b * tangent_elem, - err, sizeof(err)) != 0) - throw std::runtime_error(string("Plugin tangent evaluation failed: ") + err); + if (plugin.evaluate(model, ne * n_gp, this->time_old, this->time, gradient_buf.data(), voxel_buf.data(), phase_buf.data(), + history_old_buf.data(), flux_buf.data(), history_new_buf.data(), tangent_buf.data(), err, sizeof(err)) != 0) + throw std::runtime_error(string("Plugin material tangent failed: ") + err); + + for (size_t k = 0; k < ne; ++k) + std::copy_n(&tangent_buf[k * per_elem * n_str], per_elem * n_str, tangent_all + elems[b + k] * per_elem * n_str); } + return true; } // Writes every history variable, e.g. state/internal/Ep as state_internal_Ep @@ -316,9 +323,8 @@ class PluginModel : public Base { vector> history_vars; // name, doubles per point int n_history{0}; // history doubles per point size_t n_phase_rows{SIZE_MAX}; // phases below have a row in every per-phase field - vector trial, converged; // history at every Gauss point + vector trial, converged, previous; // history at every Gauss point vector gradient_buf, flux_buf, history_old_buf, history_new_buf; // one plugin call - vector previous, gradient_field; vector voxel_buf, phase_buf; const double *flux_cache{nullptr}; }; diff --git a/src/tangent.cpp b/src/tangent.cpp index 44d5e6e..1c4f9a9 100644 --- a/src/tangent.cpp +++ b/src/tangent.cpp @@ -1,87 +1,79 @@ +// The homogenized tangent of a converged step from the materials' own tangents. +// +// With the tangent C = d flux / d gradient of that step frozen at every Gauss +// point, the linearised microstructure is a linear material, flux = C gradient. +// FANS homogenizes it as any linear material: a solve per unit macroscopic +// gradient, whose mean flux is a column of the homogenized tangent. + #include "solver.h" +#include "LargeStrainMechModel.h" -template -using MacroVector = Matrix; +namespace { -template -using ElementVector = Matrix; +template +using Kinematics = std::conditional_t>; -template -static void tangent_residual(Solver &solver, RealArray &residual, RealArray &delta_u, - const double *tangent, const MacroVector ¯o) -{ - auto *model = solver.matmanager->models[0]; - solver.template compute_residual_basic( - residual, delta_u, [&](ElementVector &ue, int, ptrdiff_t element) -> ElementVector & { - return model->element_tangent_residual(ue, macro, tangent + element * model->n_gp * n_str * n_str); - }); -} +// flux = C gradient, with C [element][n_gp][n_str][n_str] +template +class FrozenTangent : public Kinematics { + public: + FrozenTangent(const Reader &reader, size_t n_elem) + : Kinematics(reader), C(n_elem * this->n_gp * n_str * n_str) {} -template -static MacroVector tangent_average(Solver &solver, RealArray &delta_u, RealArray &scratch, - const double *tangent, const MacroVector ¯o) -{ - auto *model = solver.matmanager->models[0]; - MacroVector average = MacroVector::Zero(); - ElementVector zero = ElementVector::Zero(); - solver.template compute_residual_basic<0>(scratch, delta_u, [&](ElementVector &ue, int, ptrdiff_t element) -> ElementVector & { - average += model->element_tangent_stress(ue, macro, tangent + element * model->n_gp * n_str * n_str); - return zero; - }); - MPI_Allreduce(MPI_IN_PLACE, average.data(), n_str, MPI_DOUBLE, MPI_SUM, solver.communicator); - return average / (solver.n_x * solver.n_y * solver.n_z); -} + vector C; -template -static void solve_tangent(Solver &solver, RealArray &delta_u, RealArray &r, RealArray &d, - RealArray &Kd, RealArray &z, const double *tangent, const MacroVector ¯o) -{ - tangent_residual(solver, r, delta_u, tangent, macro); - d.setZero(); - double rho = 0.0, error = std::sqrt(solver.dotProduct(r, r)); - const double limit = solver.reader.errorParameters["type"] == "relative" ? solver.TOL * error : solver.TOL; - for (int i = 0; i < solver.n_it && error > limit; ++i) { - solver.apply_preconditioner(r.data(), z.data()); - z *= -1.0; - const double old_rho = rho; - rho = solver.dotProduct(r, z); - d = z + (i ? rho / old_rho : 0.0) * d; - tangent_residual(solver, Kd, d, tangent, MacroVector::Zero()); - const double alpha = rho / solver.dotProduct(d, Kd); - delta_u -= alpha * d; - r -= alpha * Kd; - error = std::sqrt(solver.dotProduct(r, r)); + void get_sigma(int i, int, ptrdiff_t element_idx) override + { + const double *C_gp = &C[(element_idx * this->n_gp * n_str + i) * n_str]; + this->sigma.template segment(i).noalias() = + Map>(C_gp) * this->eps.template segment(i); } -} + Matrix get_reference_stiffness() override + { + throw std::logic_error("FrozenTangent: the solver already has its reference material"); + } +}; + +} // namespace + +// False if a material has no tangent of its own: then get_homogenized_tangent perturbs template -MatrixXd compute_consistent_homogenized_tangent(Solver &solver) +bool Solver::consistent_homogenized_tangent() { - if (solver.isMixedBCActive()) - throw std::runtime_error("Homogenized tangent is not supported with mixed boundary conditions."); - const auto n_elem = solver.local_n0 * solver.n_y * solver.n_z; - const size_t n_work = (solver.local_n0 + 1) * solver.n_y * solver.n_z * howmany; - const size_t n_fft = std::max(solver.reader.alloc_local * 2, (solver.local_n0 + 1) * solver.n_y * (solver.n_z + 2) * howmany); - vector displacement(n_work), search_direction(n_work); - unique_ptr residual(fftw_alloc_real(n_fft), fftw_free), product(fftw_alloc_real(n_fft), fftw_free); - vector local_tangent; - auto *model = solver.matmanager->models[0]; - model->compute_tangent_field(n_elem, solver.ms, local_tangent); - RealArray delta_u(displacement.data(), solver.n_z * howmany, solver.local_n0 * solver.n_y, OuterStride<>(solver.n_z * howmany)); - RealArray r(residual.get(), solver.n_z * howmany, solver.local_n0 * solver.n_y, OuterStride<>((solver.n_z + 2) * howmany)); - RealArray d(search_direction.data(), solver.n_z * howmany, solver.local_n0 * solver.n_y, OuterStride<>(solver.n_z * howmany)); - RealArray Kd(product.get(), solver.n_z * howmany, solver.local_n0 * solver.n_y, OuterStride<>(solver.n_z * howmany)); - RealArray z(product.get(), solver.n_z * howmany, solver.local_n0 * solver.n_y, OuterStride<>((solver.n_z + 2) * howmany)); - MatrixXd homogenized_tangent(n_str, n_str); + if (!matmanager->any_batched) + return false; + + // The materials' tangents at the converged step, whose nodal values batch_ue still holds + const size_t n_elem = local_n0 * n_y * n_z; + auto *frozen = new FrozenTangent(reader, n_elem); + MaterialManager linearised(frozen, reader.n_mat); // owns `frozen` + for (size_t m = 0; m < matmanager->models.size(); ++m) + if (!matmanager->models[m]->evaluate_tangent(batch_elems[m], ms, batch_ue.data(), frozen->C.data())) + return false; + + Log::logger().info("# Homogenized tangent from the materials' tangents: {} linear solves", n_str); + MaterialManager *materials = matmanager; + const vector u(v_u, v_u + n_elem * howmany); + matmanager = &linearised; + disableMixedBC(); + + homogenized_tangent.resize(n_str, n_str); for (int j = 0; j < n_str; ++j) { - delta_u.setZero(); - const MacroVector direction = MacroVector::Unit(j); - solve_tangent(solver, delta_u, r, d, Kd, z, local_tangent.data(), direction); - homogenized_tangent.col(j) = tangent_average(solver, delta_u, Kd, local_tangent.data(), direction); + vector unit(n_str, 0.0); + unit[j] = 1.0; + matmanager->set_gradient(unit); + std::fill_n(v_u, u.size(), 0.0); + solve(); + homogenized_tangent.col(j) = get_homogenized_stress(); } - return homogenized_tangent; + + // Back to the materials and the converged step + matmanager = materials; + std::copy(u.begin(), u.end(), v_u); + return true; } -template MatrixXd compute_consistent_homogenized_tangent(Solver<1, 3> &); -template MatrixXd compute_consistent_homogenized_tangent(Solver<3, 6> &); -template MatrixXd compute_consistent_homogenized_tangent(Solver<3, 9> &); +template bool Solver<1, 3>::consistent_homogenized_tangent(); +template bool Solver<3, 6>::consistent_homogenized_tangent(); +template bool Solver<3, 9>::consistent_homogenized_tangent(); From 0d73cc63ccbb1df4a3cba48e435b0655b5813fd6 Mon Sep 17 00:00:00 2001 From: Sanath Date: Tue, 6 Oct 2026 21:14:48 +0200 Subject: [PATCH 3/8] further simplifications to the design --- CMakeLists.txt | 3 +- include/matmodel.h | 36 ++++++++--- include/solver.h | 46 +------------- src/neml2_plugin.cpp | 20 +++--- src/plugin_material.cpp | 56 +++++------------ src/tangent.cpp | 133 +++++++++++++++++++++++++++++----------- 6 files changed, 156 insertions(+), 138 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 24bae8c..061fd3e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -156,7 +156,8 @@ target_include_directories(FANS_main PRIVATE "${PROJECT_BINARY_DIR}/include") # SOURCES # ############################################################################## -target_sources(FANS_FANS PRIVATE src/reader.cpp src/plugin_material.cpp src/tangent.cpp) +target_sources(FANS_FANS PRIVATE src/reader.cpp src/plugin_material.cpp + src/tangent.cpp) target_sources(FANS_main PRIVATE src/main.cpp) # ############################################################################## diff --git a/include/matmodel.h b/include/matmodel.h index 7d5e32d..1cc9ba4 100644 --- a/include/matmodel.h +++ b/include/matmodel.h @@ -27,6 +27,21 @@ class Matmodel { void getStrainStress(double *strain, double *stress, Matrix &ue, int mat_index, ptrdiff_t element_idx); void setGradient(vector _g0); + // As getStrainStress, for the tangents C = d sigma / d eps of the converged step at the + // element's Gauss points, [n_gp][n_str][n_str] + using Tangent = Map>; + void getTangent(double *tangent, Matrix &ue, int mat_index, ptrdiff_t element_idx) + { + eps.noalias() = B * ue + g0; + for (int i = 0; i < n_gp; ++i) + get_tangent(n_str * i, mat_index, element_idx, Tangent(tangent + i * n_str * n_str)); + } + // A nonlinear model that has a get_tangent says so here + virtual bool has_tangent() const + { + return false; + } + // Accessors for internal Gauss point data (populated after getStrainStress call) inline const double *get_eps_data() const { @@ -52,15 +67,10 @@ class Matmodel { { return false; } + // Unless null, tangent_gp gets the tangent d flux / d gradient of the converged + // step, [element][n_gp][n_str][n_str] (src/tangent.cpp). virtual void evaluate_batch(const vector &elems, const unsigned short *phase, - const double *ue, double *sig_gp) {} - // The tangent d flux / d gradient of the converged step at the Gauss points of `elems`, - // [element][n_gp][n_str][n_str]; false if the model has none. - virtual bool evaluate_tangent(const vector &elems, const unsigned short *phase, - const double *ue, double *tangent_gp) - { - return false; - } + const double *ue, double *sig_gp, double *tangent_gp = nullptr) {} // A model without the element stiffnesses of a LinearModel can still say that its // flux is linear in the gradient: then the linear CG applies. @@ -94,6 +104,16 @@ class Matmodel { void Construct_B(); virtual void get_sigma(int i, int mat_index, ptrdiff_t element_idx) = 0; + // As get_sigma, for the tangent C = d sigma / d eps at that Gauss point. This one is that + // of a linear model, its stresses at unit strains; a nonlinear model writes its own. + virtual void get_tangent(int i, int mat_index, ptrdiff_t element_idx, Tangent C) + { + for (int j = 0; j < n_str; ++j) { + eps.segment(i) = Matrix::Unit(j); + get_sigma(i, mat_index, element_idx); + C.col(j) = sigma.segment(i); + } + } }; template diff --git a/include/solver.h b/include/solver.h index 77a601b..a614111 100644 --- a/include/solver.h +++ b/include/solver.h @@ -4,8 +4,6 @@ #include "matmodel.h" #include "MaterialManager.h" -class J2Plasticity; - typedef Map, Unaligned, OuterStride<>> RealArray; template @@ -73,8 +71,7 @@ class Solver : private MixedBCController { void evaluate_batched_stress(double *u); //!< in src/plugin_material.cpp MatrixXd homogenized_tangent; - MatrixXd get_homogenized_tangent(double pert_param); - bool consistent_homogenized_tangent(); //!< in src/tangent.cpp + MatrixXd get_homogenized_tangent(double pert_param); //!< in src/tangent.cpp void enableMixedBC(const MixedBC &mbc, size_t step) { @@ -740,47 +737,6 @@ VectorXd Solver::get_homogenized_stress() return homogenized_stress; } -template -MatrixXd Solver::get_homogenized_tangent(double pert_param) -{ - if (consistent_homogenized_tangent()) // from the materials' own tangents, if they all have one - return homogenized_tangent; - - homogenized_tangent = MatrixXd::Zero(n_str, n_str); - VectorXd unperturbed_stress = get_homogenized_stress(); - VectorXd perturbed_stress; - vector pert_strain(n_str, 0.0); - vector g0 = matmanager->get_info(0).model->macroscale_loading; - bool islinear = matmanager->all_linear; - - for (auto *model : matmanager->models) { - if (dynamic_cast(model) != nullptr) { - throw std::runtime_error("Homogenized tangent computation not implemented for J2Plasticity models."); - } - } - // TODO: a deep copy of the solver object is needed here to avoid modifying the history of the solver object - - for (int i = 0; i < n_str; i++) { - if (islinear) { - pert_strain = vector(n_str, 0.0); - pert_strain[i] = 1.0; - } else { - pert_strain = g0; - pert_strain[i] += pert_param; - } - - matmanager->set_gradient(pert_strain); - disableMixedBC(); - solve(); - perturbed_stress = get_homogenized_stress(); - - homogenized_tangent.col(i) = islinear ? perturbed_stress : (perturbed_stress - unperturbed_stress) / pert_param; - } - - homogenized_tangent = 0.5 * (homogenized_tangent + homogenized_tangent.transpose()).eval(); - return homogenized_tangent; -} - template Solver::~Solver() { diff --git a/src/neml2_plugin.cpp b/src/neml2_plugin.cpp index fbea787..8e9b24b 100644 --- a/src/neml2_plugin.cpp +++ b/src/neml2_plugin.cpp @@ -58,6 +58,7 @@ struct FANSPluginModel { std::unique_ptr model; at::Device device{at::kCPU}; Var gradient, flux; + Var tangent; // d flux / d gradient: the flux's shape, then the gradient's std::vector history, fields; int64_t n_history{0}; // history doubles per point at::Tensor perm = at::tensor({0, 1, 2, 5, 4, 3}, at::kLong); // FANS <-> NEML2 Mandel order @@ -117,8 +118,11 @@ FANSPluginModel *fans_plugin_load(const char *config, char *msg, size_t msglen) throw std::runtime_error("the artifact has no input '" + gradient + "' (the \"gradient\")"); if (fi < 0) throw std::runtime_error("the artifact has no output '" + flux + "' (the \"flux\")"); - h->gradient = make_var(gradient, in_shape[gi]); - h->flux = make_var(flux, out_shape[fi]); + h->gradient = make_var(gradient, in_shape[gi]); + h->flux = make_var(flux, out_shape[fi]); + std::vector tangent_shape = out_shape[fi]; + tangent_shape.insert(tangent_shape.end(), in_shape[gi].begin(), in_shape[gi].end()); + h->tangent = make_var("tangent", tangent_shape); if (h->flux.size != h->gradient.size) throw std::runtime_error("the gradient '" + gradient + "' and the flux '" + flux + "' differ in size"); @@ -199,14 +203,12 @@ int fans_plugin_evaluate(FANSPluginModel *m, size_t n_points, double t_old, doub } std::map outputs; - if (!tangent) { + if (tangent) { // the same evaluation, with its Jacobian + neml2::aoti::VariablePairJacobian jacobian; + std::tie(outputs, jacobian) = m->model->jacobian(inputs); + m->to_fans(jacobian.at(m->flux.name).at(m->gradient.name), m->tangent, view(tangent, m->tangent.size)); + } else { outputs = m->model->forward(inputs); - } else { // the same evaluation with its Jacobian, whose flux-gradient block has the flux's shape, then the gradient's - auto [values, jacobian] = m->model->jacobian(inputs); - outputs = std::move(values); - Var d = m->flux; - d.shape.insert(d.shape.end(), m->gradient.shape.begin() + 1, m->gradient.shape.end()); - m->to_fans(jacobian.at(m->flux.name).at(m->gradient.name), d, view(tangent, m->flux.size * m->gradient.size)); } m->to_fans(outputs.at(m->flux.name), m->flux, view(flux, m->flux.size)); diff --git a/src/plugin_material.cpp b/src/plugin_material.cpp index b0e01b5..2cd7518 100644 --- a/src/plugin_material.cpp +++ b/src/plugin_material.cpp @@ -10,7 +10,7 @@ // the state a step starts from, and `trial`, the state of the step itself. // Once a step has converged, the solver evaluates it a last time (for the // output) and `trial` becomes `converged`; `previous` keeps the state that -// step started from, for its tangent (evaluate_tangent). +// step started from, for its tangent. // Fields: further model inputs, e.g. an orientation, read from the // microstructure file per voxel or per phase (read_field). // Linear: "linear": true in the material properties says that the flux is @@ -201,21 +201,30 @@ class PluginModel : public Base { return true; } + bool has_tangent() const override + { + return true; + } + // "linear": true in the material properties; FANS trusts it bool is_linear() const override { return linear; } - // The flux at all Gauss points of `elems`, into flux_all, and their state, into `trial` + // The flux at all Gauss points of `elems`, into flux_all, and their state, into `trial`. + // With tangent_all, the converged step once more, from the state it started from, for its tangent void evaluate_batch(const vector &elems, const unsigned short *phase, const double *ue_all, - double *flux_all) override + double *flux_all, double *tangent_all) override { const size_t n_dof = size_t(howmany) * 8; // nodal values per element const size_t per_elem = size_t(n_gp) * n_str; const size_t hist_elem = size_t(n_gp) * n_history; flux_cache = flux_all; + const vector &old = tangent_all ? previous : converged; + vector tangent_buf(tangent_all ? chunk * per_elem * n_str : 0); + for (size_t b = 0; b < elems.size(); b += chunk) { const size_t ne = std::min(chunk, elems.size() - b); @@ -227,57 +236,26 @@ class PluginModel : public Base { std::fill_n(&phase_buf[k * n_gp], n_gp, phase[e]); if (phase[e] >= n_phase_rows) throw std::runtime_error("Phase " + std::to_string(phase[e]) + " has no row in the per-phase fields of its plugin material."); - std::copy_n(converged.data() + e * hist_elem, hist_elem, history_old_buf.data() + k * hist_elem); + std::copy_n(old.data() + e * hist_elem, hist_elem, history_old_buf.data() + k * hist_elem); std::copy_n(trial.data() + e * hist_elem, hist_elem, history_new_buf.data() + k * hist_elem); // initial guess } char err[FANS_PLUGIN_MSGLEN] = {0}; if (plugin.evaluate(model, ne * n_gp, this->time_old, this->time, gradient_buf.data(), voxel_buf.data(), phase_buf.data(), - history_old_buf.data(), flux_buf.data(), history_new_buf.data(), nullptr, err, sizeof(err)) != 0) + history_old_buf.data(), flux_buf.data(), history_new_buf.data(), tangent_all ? tangent_buf.data() : nullptr, + err, sizeof(err)) != 0) throw std::runtime_error(string("Plugin material evaluation failed: ") + err); for (size_t k = 0; k < ne; ++k) { const ptrdiff_t e = elems[b + k]; std::copy_n(&flux_buf[k * per_elem], per_elem, flux_all + e * per_elem); std::copy_n(history_new_buf.data() + k * hist_elem, hist_elem, trial.data() + e * hist_elem); + if (tangent_all) + std::copy_n(&tangent_buf[k * per_elem * n_str], per_elem * n_str, tangent_all + e * per_elem * n_str); } } } - // The tangent of the converged step at all Gauss points of `elems`, into tangent_all: - // that step once more, from the state it started from (`previous`) to the one it reached - bool evaluate_tangent(const vector &elems, const unsigned short *phase, const double *ue_all, - double *tangent_all) override - { - const size_t n_dof = size_t(howmany) * 8; - const size_t per_elem = size_t(n_gp) * n_str; - const size_t hist_elem = size_t(n_gp) * n_history; - vector tangent_buf(chunk * per_elem * n_str); - - for (size_t b = 0; b < elems.size(); b += chunk) { - const size_t ne = std::min(chunk, elems.size() - b); - - for (size_t k = 0; k < ne; ++k) { - const ptrdiff_t e = elems[b + k]; - Eigen::Map(&gradient_buf[k * per_elem], per_elem).noalias() = - B * Eigen::Map(ue_all + e * n_dof, n_dof) + g0; - std::fill_n(&voxel_buf[k * n_gp], n_gp, int(e)); - std::fill_n(&phase_buf[k * n_gp], n_gp, phase[e]); - std::copy_n(previous.data() + e * hist_elem, hist_elem, history_old_buf.data() + k * hist_elem); - std::copy_n(converged.data() + e * hist_elem, hist_elem, history_new_buf.data() + k * hist_elem); // initial guess - } - - char err[FANS_PLUGIN_MSGLEN] = {0}; - if (plugin.evaluate(model, ne * n_gp, this->time_old, this->time, gradient_buf.data(), voxel_buf.data(), phase_buf.data(), - history_old_buf.data(), flux_buf.data(), history_new_buf.data(), tangent_buf.data(), err, sizeof(err)) != 0) - throw std::runtime_error(string("Plugin material tangent failed: ") + err); - - for (size_t k = 0; k < ne; ++k) - std::copy_n(&tangent_buf[k * per_elem * n_str], per_elem * n_str, tangent_all + elems[b + k] * per_elem * n_str); - } - return true; - } - // Writes every history variable, e.g. state/internal/Ep as state_internal_Ep void postprocess(Solver &solver, Reader &reader, int load_idx, int time_idx) override { diff --git a/src/tangent.cpp b/src/tangent.cpp index 1c4f9a9..266ef65 100644 --- a/src/tangent.cpp +++ b/src/tangent.cpp @@ -1,19 +1,23 @@ -// The homogenized tangent of a converged step from the materials' own tangents. +// The homogenized tangent of a converged step. // -// With the tangent C = d flux / d gradient of that step frozen at every Gauss -// point, the linearised microstructure is a linear material, flux = C gradient. -// FANS homogenizes it as any linear material: a solve per unit macroscopic -// gradient, whose mean flux is a column of the homogenized tangent. +// If every material has a tangent C = d flux / d gradient, the tangents of the +// step, frozen at every Gauss point, are a linear material: flux = C gradient. +// Its homogenized stiffness is the consistent homogenized tangent, and FANS +// finds it as for any linear material: a solve per unit macroscopic gradient, +// whose mean flux is a column. +// Otherwise the materials themselves are perturbed, by a finite difference. #include "solver.h" #include "LargeStrainMechModel.h" +#include "material_models/small_strain/J2Plasticity.h" namespace { +// The B matrices of the problem template using Kinematics = std::conditional_t>; -// flux = C gradient, with C [element][n_gp][n_str][n_str] +// The linear material flux = C gradient, C being [element][n_gp][n_str][n_str] template class FrozenTangent : public Kinematics { public: @@ -24,56 +28,113 @@ class FrozenTangent : public Kinematics { void get_sigma(int i, int, ptrdiff_t element_idx) override { - const double *C_gp = &C[(element_idx * this->n_gp * n_str + i) * n_str]; - this->sigma.template segment(i).noalias() = - Map>(C_gp) * this->eps.template segment(i); + const size_t gp = element_idx * this->n_gp + i / n_str; + Map> C_gp(&C[gp * n_str * n_str]); + this->sigma.template segment(i).noalias() = C_gp * this->eps.template segment(i); } + // Not needed: the solver has its reference material already Matrix get_reference_stiffness() override { - throw std::logic_error("FrozenTangent: the solver already has its reference material"); + return Matrix::Zero(); } }; +// The tangent of every material at the converged step, into C +template +void collect_tangents(Solver &solver, vector &C) +{ + const auto &materials = *solver.matmanager; + const ptrdiff_t n_slab = solver.n_y * solver.n_z * howmany; // nodal values of a slab of the grid + const size_t per_elem = C.size() / (solver.local_n0 * solver.n_y * solver.n_z); + double *u = solver.v_u; + + // The native materials element by element, as Solver::get_homogenized_stress goes through them + MPI_Sendrecv(u, n_slab, MPI_DOUBLE, (solver.world_rank + solver.world_size - 1) % solver.world_size, 0, + u + solver.local_n0 * n_slab, n_slab, MPI_DOUBLE, (solver.world_rank + 1) % solver.world_size, 0, solver.communicator, MPI_STATUS_IGNORE); + Matrix ue; + solver.template iterateCubes<0>([&](ptrdiff_t *idx, ptrdiff_t *) { + const MaterialInfo &info = materials.get_info(solver.ms[idx[0]]); + if (info.model->wants_batch()) + return; + for (int i = 0; i < 8; ++i) + for (int j = 0; j < howmany; ++j) + ue(howmany * i + j) = u[howmany * idx[i] + j]; + info.model->getTangent(&C[idx[0] * per_elem], ue, info.local_mat_id, idx[0]); + }); + + // The batched materials in batches, from the nodal values of the step, which batch_ue still holds + for (size_t m = 0; m < materials.models.size(); ++m) + if (materials.models[m]->wants_batch()) + materials.models[m]->evaluate_batch(solver.batch_elems[m], solver.ms, solver.batch_ue.data(), solver.batch_gp_stress.data(), C.data()); +} + } // namespace -// False if a material has no tangent of its own: then get_homogenized_tangent perturbs template -bool Solver::consistent_homogenized_tangent() +MatrixXd Solver::get_homogenized_tangent(double pert_param) { - if (!matmanager->any_batched) - return false; - - // The materials' tangents at the converged step, whose nodal values batch_ue still holds - const size_t n_elem = local_n0 * n_y * n_z; - auto *frozen = new FrozenTangent(reader, n_elem); - MaterialManager linearised(frozen, reader.n_mat); // owns `frozen` - for (size_t m = 0; m < matmanager->models.size(); ++m) - if (!matmanager->models[m]->evaluate_tangent(batch_elems[m], ms, batch_ue.data(), frozen->C.data())) - return false; - - Log::logger().info("# Homogenized tangent from the materials' tangents: {} linear solves", n_str); + // The step, to come back to MaterialManager *materials = matmanager; - const vector u(v_u, v_u + n_elem * howmany); - matmanager = &linearised; + const vector g0 = matmanager->models[0]->macroscale_loading; + const vector u(v_u, v_u + local_n0 * n_y * n_z * howmany); + + // Does every material have a tangent? Linear models with element stiffnesses are + // left as they are: perturbing them is exact, and faster. + bool consistent = !matmanager->all_stiffness; + for (int phase = 0; phase < reader.n_mat; ++phase) { + const MaterialInfo &info = matmanager->get_info(phase); + consistent = consistent && (info.is_linear || info.model->has_tangent()); + } + + // If so, the frozen tangents take the place of the materials + std::unique_ptr> linearised; + if (consistent) { + auto *frozen = new FrozenTangent(reader, local_n0 * n_y * n_z); + collect_tangents(*this, frozen->C); + linearised = std::make_unique>(frozen, reader.n_mat); // owns `frozen` + matmanager = linearised.get(); + } else { + for (auto *model : matmanager->models) // perturbing would change its history + if (dynamic_cast(model) != nullptr) + throw std::runtime_error("Homogenized tangent computation not implemented for J2Plasticity models."); + } + Log::logger().info("# Homogenized tangent {}: {} solves", consistent ? "from the materials' tangents" : "by perturbation", n_str); + + const bool linear = matmanager->all_linear; + VectorXd stress; + if (!linear) + stress = get_homogenized_stress(); disableMixedBC(); homogenized_tangent.resize(n_str, n_str); for (int j = 0; j < n_str; ++j) { - vector unit(n_str, 0.0); - unit[j] = 1.0; - matmanager->set_gradient(unit); - std::fill_n(v_u, u.size(), 0.0); + vector gradient(n_str, 0.0); + if (linear) { // a unit gradient, from rest + gradient[j] = 1.0; + std::fill_n(v_u, u.size(), 0.0); + } else { // a small step from the gradient of the step + gradient = g0; + gradient[j] += pert_param; + } + matmanager->set_gradient(gradient); solve(); - homogenized_tangent.col(j) = get_homogenized_stress(); + + if (linear) + homogenized_tangent.col(j) = get_homogenized_stress(); + else + homogenized_tangent.col(j) = (get_homogenized_stress() - stress) / pert_param; } + if (!consistent) + homogenized_tangent = 0.5 * (homogenized_tangent + homogenized_tangent.transpose()).eval(); - // Back to the materials and the converged step + // Back to the step matmanager = materials; + matmanager->set_gradient(g0); std::copy(u.begin(), u.end(), v_u); - return true; + return homogenized_tangent; } -template bool Solver<1, 3>::consistent_homogenized_tangent(); -template bool Solver<3, 6>::consistent_homogenized_tangent(); -template bool Solver<3, 9>::consistent_homogenized_tangent(); +template MatrixXd Solver<1, 3>::get_homogenized_tangent(double); +template MatrixXd Solver<3, 6>::get_homogenized_tangent(double); +template MatrixXd Solver<3, 9>::get_homogenized_tangent(double); From 31409dbaeec7daa33f0127843e913367e1c816f6 Mon Sep 17 00:00:00 2001 From: Sanath Date: Tue, 6 Oct 2026 22:04:52 +0200 Subject: [PATCH 4/8] relaxed checks in get_homogenized_tangent --- include/matmodel.h | 5 +++-- include/solver.h | 3 +-- src/neml2_plugin.cpp | 9 +++++---- src/plugin_material.cpp | 5 ----- src/tangent.cpp | 27 +++++++++++---------------- 5 files changed, 20 insertions(+), 29 deletions(-) diff --git a/include/matmodel.h b/include/matmodel.h index 1cc9ba4..18b88a5 100644 --- a/include/matmodel.h +++ b/include/matmodel.h @@ -36,10 +36,11 @@ class Matmodel { for (int i = 0; i < n_gp; ++i) get_tangent(n_str * i, mat_index, element_idx, Tangent(tangent + i * n_str * n_str)); } - // A nonlinear model that has a get_tangent says so here + // A nonlinear model with a tangent: a batched one gives it in evaluate_batch, + // any other writes a get_tangent and says so here virtual bool has_tangent() const { - return false; + return wants_batch(); } // Accessors for internal Gauss point data (populated after getStrainStress call) diff --git a/include/solver.h b/include/solver.h index a614111..92ea030 100644 --- a/include/solver.h +++ b/include/solver.h @@ -698,8 +698,7 @@ void Solver::postprocess(Reader &reader, int load_idx, int time_ << std::setprecision(12) << homogenized_tangent << '\n'; Log::logger().info("{}", output.str()); } - const Matrix rows = homogenized_tangent; // as the file stores it - reader.writeData("homogenized_tangent", load_idx, time_idx, rows.data(), dims, 2); + reader.writeData("homogenized_tangent", load_idx, time_idx, homogenized_tangent.data(), dims, 2); } } diff --git a/src/neml2_plugin.cpp b/src/neml2_plugin.cpp index 8e9b24b..8d3538f 100644 --- a/src/neml2_plugin.cpp +++ b/src/neml2_plugin.cpp @@ -118,13 +118,14 @@ FANSPluginModel *fans_plugin_load(const char *config, char *msg, size_t msglen) throw std::runtime_error("the artifact has no input '" + gradient + "' (the \"gradient\")"); if (fi < 0) throw std::runtime_error("the artifact has no output '" + flux + "' (the \"flux\")"); - h->gradient = make_var(gradient, in_shape[gi]); - h->flux = make_var(flux, out_shape[fi]); + h->gradient = make_var(gradient, in_shape[gi]); + h->flux = make_var(flux, out_shape[fi]); + if (h->flux.size != h->gradient.size) + throw std::runtime_error("the gradient '" + gradient + "' and the flux '" + flux + "' differ in size"); + std::vector tangent_shape = out_shape[fi]; tangent_shape.insert(tangent_shape.end(), in_shape[gi].begin(), in_shape[gi].end()); h->tangent = make_var("tangent", tangent_shape); - if (h->flux.size != h->gradient.size) - throw std::runtime_error("the gradient '" + gradient + "' and the flux '" + flux + "' differ in size"); for (size_t i = 0; i < in.size(); ++i) { const std::string &name = in[i]; diff --git a/src/plugin_material.cpp b/src/plugin_material.cpp index 2cd7518..b56f984 100644 --- a/src/plugin_material.cpp +++ b/src/plugin_material.cpp @@ -201,11 +201,6 @@ class PluginModel : public Base { return true; } - bool has_tangent() const override - { - return true; - } - // "linear": true in the material properties; FANS trusts it bool is_linear() const override { diff --git a/src/tangent.cpp b/src/tangent.cpp index 266ef65..d42980d 100644 --- a/src/tangent.cpp +++ b/src/tangent.cpp @@ -5,11 +5,11 @@ // Its homogenized stiffness is the consistent homogenized tangent, and FANS // finds it as for any linear material: a solve per unit macroscopic gradient, // whose mean flux is a column. -// Otherwise the materials themselves are perturbed, by a finite difference. +// Otherwise the materials themselves are perturbed, by a finite difference. That +// is not for materials with internal variables: every perturbed solve advances them. #include "solver.h" #include "LargeStrainMechModel.h" -#include "material_models/small_strain/J2Plasticity.h" namespace { @@ -82,10 +82,9 @@ MatrixXd Solver::get_homogenized_tangent(double pert_param) // Does every material have a tangent? Linear models with element stiffnesses are // left as they are: perturbing them is exact, and faster. bool consistent = !matmanager->all_stiffness; - for (int phase = 0; phase < reader.n_mat; ++phase) { - const MaterialInfo &info = matmanager->get_info(phase); - consistent = consistent && (info.is_linear || info.model->has_tangent()); - } + for (int phase = 0; phase < reader.n_mat; ++phase) + if (!matmanager->get_info(phase).is_linear && !matmanager->get_info(phase).model->has_tangent()) + consistent = false; // If so, the frozen tangents take the place of the materials std::unique_ptr> linearised; @@ -94,12 +93,10 @@ MatrixXd Solver::get_homogenized_tangent(double pert_param) collect_tangents(*this, frozen->C); linearised = std::make_unique>(frozen, reader.n_mat); // owns `frozen` matmanager = linearised.get(); - } else { - for (auto *model : matmanager->models) // perturbing would change its history - if (dynamic_cast(model) != nullptr) - throw std::runtime_error("Homogenized tangent computation not implemented for J2Plasticity models."); + Log::logger().info("# Homogenized tangent from the materials' tangents"); + } else if (!matmanager->all_linear) { + Log::logger().warn("# Homogenized tangent by perturbation: For materials with internal variables, this is WRONG!"); } - Log::logger().info("# Homogenized tangent {}: {} solves", consistent ? "from the materials' tangents" : "by perturbation", n_str); const bool linear = matmanager->all_linear; VectorXd stress; @@ -110,10 +107,10 @@ MatrixXd Solver::get_homogenized_tangent(double pert_param) homogenized_tangent.resize(n_str, n_str); for (int j = 0; j < n_str; ++j) { vector gradient(n_str, 0.0); - if (linear) { // a unit gradient, from rest + if (linear) { gradient[j] = 1.0; std::fill_n(v_u, u.size(), 0.0); - } else { // a small step from the gradient of the step + } else { gradient = g0; gradient[j] += pert_param; } @@ -125,12 +122,10 @@ MatrixXd Solver::get_homogenized_tangent(double pert_param) else homogenized_tangent.col(j) = (get_homogenized_stress() - stress) / pert_param; } - if (!consistent) - homogenized_tangent = 0.5 * (homogenized_tangent + homogenized_tangent.transpose()).eval(); + homogenized_tangent = 0.5 * (homogenized_tangent + homogenized_tangent.transpose()).eval(); // Back to the step matmanager = materials; - matmanager->set_gradient(g0); std::copy(u.begin(), u.end(), v_u); return homogenized_tangent; } From d8dcb2e84fa7044387d8617981f8bbf3459dd890 Mon Sep 17 00:00:00 2001 From: Sanath Date: Wed, 7 Oct 2026 23:11:12 +0200 Subject: [PATCH 5/8] cleanup --- src/neml2_plugin.cpp | 11 +++-------- src/plugin_material.cpp | 19 ++++++++----------- src/tangent.cpp | 25 ++++++++++--------------- 3 files changed, 21 insertions(+), 34 deletions(-) diff --git a/src/neml2_plugin.cpp b/src/neml2_plugin.cpp index 8d3538f..bdefb40 100644 --- a/src/neml2_plugin.cpp +++ b/src/neml2_plugin.cpp @@ -203,14 +203,9 @@ int fans_plugin_evaluate(FANSPluginModel *m, size_t n_points, double t_old, doub inputs[h.name] = m->to_neml2(new_history.narrow(1, h.offset, h.size), h); // initial guess } - std::map outputs; - if (tangent) { // the same evaluation, with its Jacobian - neml2::aoti::VariablePairJacobian jacobian; - std::tie(outputs, jacobian) = m->model->jacobian(inputs); - m->to_fans(jacobian.at(m->flux.name).at(m->gradient.name), m->tangent, view(tangent, m->tangent.size)); - } else { - outputs = m->model->forward(inputs); - } + const auto outputs = m->model->forward(inputs); + if (tangent) // from the Jacobian of the same evaluation + m->to_fans(m->model->jacobian(inputs).second.at(m->flux.name).at(m->gradient.name), m->tangent, view(tangent, m->tangent.size)); m->to_fans(outputs.at(m->flux.name), m->flux, view(flux, m->flux.size)); for (const Var &h : m->history) diff --git a/src/plugin_material.cpp b/src/plugin_material.cpp index b56f984..d90b5e5 100644 --- a/src/plugin_material.cpp +++ b/src/plugin_material.cpp @@ -50,15 +50,13 @@ Plugin load_plugin(const string &name) throw std::runtime_error("Could not load the '" + name + "' material plugin: " + dlerror() + "\nBuild FANS with -DFANS_" + name + "=ON."); - auto sym = [&](auto &f, const char *n) { - if (!(f = reinterpret_cast>(dlsym(h, n)))) - throw std::runtime_error(lib + " is missing '" + n + "'"); - }; Plugin p; - sym(p.load, "fans_plugin_load"); - sym(p.set_table, "fans_plugin_set_table"); - sym(p.evaluate, "fans_plugin_evaluate"); - sym(p.free_model, "fans_plugin_free"); + p.load = reinterpret_cast(dlsym(h, "fans_plugin_load")); + p.set_table = reinterpret_cast(dlsym(h, "fans_plugin_set_table")); + p.evaluate = reinterpret_cast(dlsym(h, "fans_plugin_evaluate")); + p.free_model = reinterpret_cast(dlsym(h, "fans_plugin_free")); + if (!p.load || !p.set_table || !p.evaluate || !p.free_model) + throw std::runtime_error(lib + " lacks a function of fans_plugin.h"); return p; } @@ -181,11 +179,10 @@ class PluginModel : public Base { return {flux_cache + size_t(element_idx) * n_gp * n_str, n_gp * n_str}; } - // Called per Gauss point: the first call brings the fluxes of them all + // The flux at a Gauss point, from the last evaluate_batch void get_sigma(int i, int, ptrdiff_t element_idx) override { - if (i == 0) - sigma = fluxes(element_idx); + sigma.segment(i, n_str) = fluxes(element_idx).segment(i, n_str); } // The element's residual straight from its fluxes: the gradient and the loop over diff --git a/src/tangent.cpp b/src/tangent.cpp index d42980d..7c77bcf 100644 --- a/src/tangent.cpp +++ b/src/tangent.cpp @@ -13,7 +13,7 @@ namespace { -// The B matrices of the problem +// The model base of the problem: thermal, small or large strain template using Kinematics = std::conditional_t>; @@ -44,14 +44,10 @@ class FrozenTangent : public Kinematics { template void collect_tangents(Solver &solver, vector &C) { - const auto &materials = *solver.matmanager; - const ptrdiff_t n_slab = solver.n_y * solver.n_z * howmany; // nodal values of a slab of the grid - const size_t per_elem = C.size() / (solver.local_n0 * solver.n_y * solver.n_z); - double *u = solver.v_u; + const auto &materials = *solver.matmanager; + const size_t per_elem = C.size() / (solver.local_n0 * solver.n_y * solver.n_z); // The native materials element by element, as Solver::get_homogenized_stress goes through them - MPI_Sendrecv(u, n_slab, MPI_DOUBLE, (solver.world_rank + solver.world_size - 1) % solver.world_size, 0, - u + solver.local_n0 * n_slab, n_slab, MPI_DOUBLE, (solver.world_rank + 1) % solver.world_size, 0, solver.communicator, MPI_STATUS_IGNORE); Matrix ue; solver.template iterateCubes<0>([&](ptrdiff_t *idx, ptrdiff_t *) { const MaterialInfo &info = materials.get_info(solver.ms[idx[0]]); @@ -59,11 +55,11 @@ void collect_tangents(Solver &solver, vector &C) return; for (int i = 0; i < 8; ++i) for (int j = 0; j < howmany; ++j) - ue(howmany * i + j) = u[howmany * idx[i] + j]; + ue(howmany * i + j) = solver.v_u[howmany * idx[i] + j]; info.model->getTangent(&C[idx[0] * per_elem], ue, info.local_mat_id, idx[0]); }); - // The batched materials in batches, from the nodal values of the step, which batch_ue still holds + // The batched materials in batches for (size_t m = 0; m < materials.models.size(); ++m) if (materials.models[m]->wants_batch()) materials.models[m]->evaluate_batch(solver.batch_elems[m], solver.ms, solver.batch_ue.data(), solver.batch_gp_stress.data(), C.data()); @@ -74,7 +70,9 @@ void collect_tangents(Solver &solver, vector &C) template MatrixXd Solver::get_homogenized_tangent(double pert_param) { - // The step, to come back to + // The step: its stress (evaluating it also readies the nodal values the tangents are + // collected from), and what to come back to + const VectorXd stress = get_homogenized_stress(); MaterialManager *materials = matmanager; const vector g0 = matmanager->models[0]->macroscale_loading; const vector u(v_u, v_u + local_n0 * n_y * n_z * howmany); @@ -82,7 +80,7 @@ MatrixXd Solver::get_homogenized_tangent(double pert_param) // Does every material have a tangent? Linear models with element stiffnesses are // left as they are: perturbing them is exact, and faster. bool consistent = !matmanager->all_stiffness; - for (int phase = 0; phase < reader.n_mat; ++phase) + for (int phase = 0; phase < matmanager->get_num_phases(); ++phase) if (!matmanager->get_info(phase).is_linear && !matmanager->get_info(phase).model->has_tangent()) consistent = false; @@ -91,7 +89,7 @@ MatrixXd Solver::get_homogenized_tangent(double pert_param) if (consistent) { auto *frozen = new FrozenTangent(reader, local_n0 * n_y * n_z); collect_tangents(*this, frozen->C); - linearised = std::make_unique>(frozen, reader.n_mat); // owns `frozen` + linearised = std::make_unique>(frozen, matmanager->get_num_phases()); // owns `frozen` matmanager = linearised.get(); Log::logger().info("# Homogenized tangent from the materials' tangents"); } else if (!matmanager->all_linear) { @@ -99,9 +97,6 @@ MatrixXd Solver::get_homogenized_tangent(double pert_param) } const bool linear = matmanager->all_linear; - VectorXd stress; - if (!linear) - stress = get_homogenized_stress(); disableMixedBC(); homogenized_tangent.resize(n_str, n_str); From 337e17efdf4ee40daca1a3a1d254f0dba01c1d99 Mon Sep 17 00:00:00 2001 From: Sanath Date: Fri, 9 Oct 2026 21:44:07 +0200 Subject: [PATCH 6/8] long overdue simplifications --- include/MaterialManager.h | 4 - .../large_strain/CompressibleNeoHookean.h | 20 +++ include/material_models/thermal/GBDiffusion.h | 4 +- include/matmodel.h | 14 +-- include/reader.h | 27 ++-- include/solver.h | 29 +++-- src/plugin_material.cpp | 7 +- src/reader.cpp | 117 +++++------------- 8 files changed, 84 insertions(+), 138 deletions(-) diff --git a/include/MaterialManager.h b/include/MaterialManager.h index 855de88..a927a83 100644 --- a/include/MaterialManager.h +++ b/include/MaterialManager.h @@ -228,10 +228,6 @@ class MaterialManager { { return phase_to_info[phase_id]; } - inline size_t get_num_models() const - { - return models.size(); - } // Get number of unique material models inline int get_num_phases() const { return n_phases; diff --git a/include/material_models/large_strain/CompressibleNeoHookean.h b/include/material_models/large_strain/CompressibleNeoHookean.h index 824c3b2..f65a9f5 100644 --- a/include/material_models/large_strain/CompressibleNeoHookean.h +++ b/include/material_models/large_strain/CompressibleNeoHookean.h @@ -47,6 +47,26 @@ class CompressibleNeoHookean : public LargeStrainMechModel { return lambda[mat_index] * logJ * C_inv + mu[mat_index] * (Matrix3d::Identity() - C_inv); } + bool has_tangent() const override + { + return true; + } + + // A = dP/dF of P = mu F + (lambda log(J) - mu) G with G = F^{-T}: + // A_iJkL = mu d_ik d_JL + lambda G_iJ G_kL + (mu - lambda log(J)) G_iL G_kJ + void get_tangent(int i_eps, int mat_index, ptrdiff_t element_idx, Tangent A) override + { + const Matrix3d F = extract_F(i_eps); + const Matrix3d G = F.inverse().transpose(); + const double c = mu[mat_index] - lambda[mat_index] * log(F.determinant()); + + for (int i = 0; i < 3; ++i) + for (int J = 0; J < 3; ++J) + for (int k = 0; k < 3; ++k) + for (int L = 0; L < 3; ++L) + A(3 * i + J, 3 * k + L) = mu[mat_index] * (i == k) * (J == L) + lambda[mat_index] * G(i, J) * G(k, L) + c * G(i, L) * G(k, J); + } + // At F = I the tangent is that of linear isotropic elasticity Matrix get_reference_stiffness() override { diff --git a/include/material_models/thermal/GBDiffusion.h b/include/material_models/thermal/GBDiffusion.h index 2021e84..8f0b824 100644 --- a/include/material_models/thermal/GBDiffusion.h +++ b/include/material_models/thermal/GBDiffusion.h @@ -74,8 +74,8 @@ class GBDiffusion : public ThermalModel, public LinearModel<1, 3> { const int num_GB = static_cast(boundary_count); n_mat = num_crystals + num_GB; - const vector grain_rot_matrices = Reader::ReadData(reader.ms_filename, reader.MSGroup() + "rotation_matrices"); - const vector GB_normals = Reader::ReadData(reader.ms_filename, reader.MSGroup() + "GB_normals"); + const vector grain_rot_matrices = Reader::ReadData(reader.ms_filename, reader.ms_group + "rotation_matrices"); + const vector GB_normals = Reader::ReadData(reader.ms_filename, reader.ms_group + "GB_normals"); if (grain_rot_matrices.size() < 9 * size_t(num_crystals) || GB_normals.size() < 3 * size_t(n_mat)) throw std::runtime_error("rotation_matrices or GB_normals has too few entries"); diff --git a/include/matmodel.h b/include/matmodel.h index 18b88a5..83f0cf4 100644 --- a/include/matmodel.h +++ b/include/matmodel.h @@ -12,13 +12,9 @@ class Matmodel { public: static constexpr int num_str = n_str; // length of strain and stress - int verbosity; //!< output verbosity - int n_mat; //!< Number of Materials - int n_gp; //!< Number of Gauss points (computed from FE_type) - string FE_type; //!< Finite element type: "HEX8", "HEX8R", "BBAR" - - double *strain; //!< Gradient - double *stress; //!< Flux + int n_mat; //!< Number of Materials + int n_gp; //!< Number of Gauss points (computed from FE_type) + string FE_type; //!< Finite element type: "HEX8", "HEX8R", "BBAR" Matmodel(const Reader &reader); @@ -52,10 +48,6 @@ class Matmodel { { return sigma.data(); } - inline int get_n_gp() const - { - return n_gp; - } virtual void postprocess(Solver &solver, Reader &reader, int load_idx, int time_idx) {}; diff --git a/include/reader.h b/include/reader.h index 2295ad5..e2db91b 100644 --- a/include/reader.h +++ b/include/reader.h @@ -25,9 +25,9 @@ class Reader { ~Reader(); // contents of input file: - char ms_filename[4096]{}; // Name of Micro-structure hdf5 file - char ms_datasetname[4096]{}; // Absolute path of Micro-structure in hdf5 file - char results_prefix[4096]{}; + string ms_filename; // Name of Micro-structure hdf5 file + string ms_datasetname; // Absolute path of Micro-structure in hdf5 file + string ms_group; // Its group, where the datasets that go with it are: "/sphere/32x32x32/" for "/sphere/32x32x32/ms" int n_mat; json inputJson; // Complete input JSON (for MaterialManager) json materialProperties; @@ -45,9 +45,8 @@ class Reader { string strain_type{"small"}; // "small" (default) or "large" string FE_type; // "HEX8" (default), "HEX8R", or "BBAR" vector resultsToWrite; - char results_filename[4096]{}; // Output HDF5 filename - char dataset_name[8192]{}; // Base path for results in HDF5 file - hid_t results_file_id = -1; // Open HDF5 file handle for results + string dataset_name; // Base path for results in HDF5 file + hid_t results_file_id = -1; // Open HDF5 file handle for results // contents of microstructure file: vector dims; @@ -66,12 +65,10 @@ class Reader { ptrdiff_t local_n1; ptrdiff_t local_1_start; - // void Setup(ptrdiff_t howmany); void ReadInputFile(const std::string &input_fn); void ReadMS(int hm); void FreeMS(); void ComputeVolumeFractions(); - // void ReadHDF5(char file_name[], char dset_name[]); void safe_create_group(hid_t file, const char *const name); void OpenResultsFile(const char *output_fn); // Open results file once void CloseResultsFile(); // Explicitly close results file @@ -94,10 +91,6 @@ class Reader { static std::vector DataShape(const string &file, const string &dset_name); template void ReadSlab(T *data, const std::vector &extra_dims, const string &file, const string &dset_name) const; - - // The microstructure's group, e.g. "/sphere/32x32x32/" for "/sphere/32x32x32/ms", - // where the datasets that go with it are - string MSGroup() const; }; // Opens a dataset for reading; the caller closes it and file_id @@ -262,9 +255,8 @@ void Reader::writeData(const char *fieldName, int load_idx, int time_idx, const if (std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) { return; } - char name[5096]; - snprintf(name, sizeof(name), "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, fieldName); - WriteData(data, name, shape, rank); + const string name = dataset_name + "/load" + to_string(load_idx) + "/time_step" + to_string(time_idx) + "/" + fieldName; + WriteData(data, name.c_str(), shape, rank); } template @@ -273,9 +265,8 @@ void Reader::writeSlab(const char *fieldName, int load_idx, int time_idx, const if (std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) { return; } - char name[5096]; - snprintf(name, sizeof(name), "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, fieldName); - WriteSlab(data, extra_dims, name); + const string name = dataset_name + "/load" + to_string(load_idx) + "/time_step" + to_string(time_idx) + "/" + fieldName; + WriteSlab(data, extra_dims, name.c_str()); } #endif diff --git a/include/solver.h b/include/solver.h index 92ea030..a779714 100644 --- a/include/solver.h +++ b/include/solver.h @@ -49,6 +49,7 @@ class Solver : private MixedBCController { template void iterateCubes(F f); + void update_ghost_layer(double *u); //!< the next rank's first layer of nodal values, behind this rank's own void solve(); void extrapolateDisplacement(); //!< Linear extrapolation for next time step @@ -227,6 +228,15 @@ void Solver::CreateFFTWPlans(double *in, fftw_complex *transform new (&rhat) Map((std::complex *) transformed, local_n1 * n_x * (n_z / 2 + 1) * howmany); } +// The elements of a rank's last layer reach into the first layer of nodes of the next rank +template +void Solver::update_ghost_layer(double *u) +{ + const int layer = n_y * n_z * howmany; + MPI_Sendrecv(u, layer, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0, + u + local_n0 * layer, layer, MPI_DOUBLE, (world_rank + 1) % world_size, 0, communicator, MPI_STATUS_IGNORE); +} + // TODO: possibly circumvent the padding problem by accessing r as a matrix? template template @@ -241,10 +251,7 @@ void Solver::compute_residual_basic(RealArray &r_matrix, RealArr r[i] = 0; } - // int MPI_Sendrecv(void *sendbuf, int sendcount, MPI_Datatype sendtype, int dest, int sendtag, void *recvbuf, - // int recvcount, MPI_Datatype recvtype, int source, int recvtag, MPI_Comm comm, MPI_Status *status) - MPI_Sendrecv(u, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0, - u + local_n0 * n_y * n_z * howmany, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + 1) % world_size, 0, communicator, MPI_STATUS_IGNORE); + update_ghost_layer(u); Matrix ue; @@ -494,8 +501,7 @@ void Solver::postprocess(Reader &reader, int load_idx, int time_ vector phase_strain_average(n_mat, VectorXd::Zero(n_str)); vector phase_counts(n_mat, 0); - MPI_Sendrecv(v_u, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0, - v_u + local_n0 * n_y * n_z * howmany, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + 1) % world_size, 0, communicator, MPI_STATUS_IGNORE); + update_ghost_layer(v_u); Matrix ue; int phase_id; @@ -705,17 +711,14 @@ void Solver::postprocess(Reader &reader, int load_idx, int time_ template VectorXd Solver::get_homogenized_stress() { - - VectorXd strain = VectorXd::Zero(local_n0 * n_y * n_z * n_str); - VectorXd stress = VectorXd::Zero(local_n0 * n_y * n_z * n_str); homogenized_stress = VectorXd::Zero(n_str); - MPI_Sendrecv(v_u, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0, - v_u + local_n0 * n_y * n_z * howmany, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + 1) % world_size, 0, communicator, MPI_STATUS_IGNORE); + update_ghost_layer(v_u); if (matmanager->any_batched) evaluate_batched_stress(v_u); Matrix ue; + Matrix strain, stress; // of one element int phase_id; iterateCubes<0>([&](ptrdiff_t *idx, ptrdiff_t *idxPadding) { for (int i = 0; i < 8; ++i) { @@ -726,8 +729,8 @@ VectorXd Solver::get_homogenized_stress() phase_id = ms[idx[0]]; const MaterialInfo &info = matmanager->get_info(phase_id); - info.model->getStrainStress(strain.segment(n_str * idx[0], n_str).data(), stress.segment(n_str * idx[0], n_str).data(), ue, info.local_mat_id, idx[0]); - homogenized_stress += stress.segment(n_str * idx[0], n_str); + info.model->getStrainStress(strain.data(), stress.data(), ue, info.local_mat_id, idx[0]); + homogenized_stress += stress; }); MPI_Allreduce(MPI_IN_PLACE, homogenized_stress.data(), n_str, MPI_DOUBLE, MPI_SUM, communicator); diff --git a/src/plugin_material.cpp b/src/plugin_material.cpp index d90b5e5..6a0b211 100644 --- a/src/plugin_material.cpp +++ b/src/plugin_material.cpp @@ -70,7 +70,7 @@ struct Field { // absolute path; per voxel if it is [Z][Y][X][...] on the grid, else per phase [n_phase][...] Field read_field(const Reader &reader, const string &name, const string &dset, int size) { - const string path = dset[0] == '/' ? dset : reader.MSGroup() + dset; + const string path = dset[0] == '/' ? dset : reader.ms_group + dset; const vector shape = Reader::DataShape(reader.ms_filename, path); const bool per_voxel = shape.size() >= 3 && shape[0] == hsize_t(reader.dims[2]) && shape[1] == hsize_t(reader.dims[1]) && shape[2] == hsize_t(reader.dims[0]); @@ -257,7 +257,7 @@ class PluginModel : public Base { if (!need && !need_gp) return; const size_t n_pts = size_t(solver.local_n0 * solver.n_y * solver.n_z) * n_gp; - const string dir = string(reader.dataset_name) + "/load" + std::to_string(load_idx) + "/time_step" + std::to_string(time_idx) + "/"; + const string dir = reader.dataset_name + "/load" + std::to_string(load_idx) + "/time_step" + std::to_string(time_idx) + "/"; size_t offset = 0; for (auto &[name, size] : history_vars) { @@ -336,8 +336,7 @@ void Solver::evaluate_batched_stress(double *u) batch_gp_stress.resize(size_t(n_elem) * models[0]->n_gp * n_str); } - MPI_Sendrecv(u, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0, - u + local_n0 * n_y * n_z * howmany, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + 1) % world_size, 0, communicator, MPI_STATUS_IGNORE); + update_ghost_layer(u); iterateCubes<0>([&](ptrdiff_t *idx, ptrdiff_t *idxPadding) { if (matmanager->get_info(ms[idx[0]]).model->wants_batch()) diff --git a/src/reader.cpp b/src/reader.cpp index 6be704c..1014865 100644 --- a/src/reader.cpp +++ b/src/reader.cpp @@ -1,15 +1,8 @@ #include "general.h" #include "reader.h" -#include "H5Cpp.h" -#include "fftw3-mpi.h" -#include "hdf5.h" -#include "mpi.h" #include -#include "H5FDmpi.h" -#include "H5FDmpio.h" - Reader::Reader(const MPI_Comm &comm) : communicator(comm) { @@ -19,47 +12,25 @@ Reader::Reader(const MPI_Comm &comm) void Reader::ComputeVolumeFractions() { - unsigned short local_max = 0; - unsigned short local_min = USHRT_MAX; - size_t local_size = local_n0 * dims[1] * dims[2]; - - // Find the local maximum and minimum material indices - for (size_t i = 0; i < local_size; i++) { - unsigned short val = static_cast(ms[i]); - if (val > local_max) { - local_max = val; - } - if (val < local_min) { - local_min = val; - } - } - - // Find the global maximum and minimum material indices - unsigned short global_max, global_min; - MPI_Allreduce(&local_max, &global_max, 1, MPI_UNSIGNED_SHORT, MPI_MAX, communicator); - MPI_Allreduce(&local_min, &global_min, 1, MPI_UNSIGNED_SHORT, MPI_MIN, communicator); - - if (global_min != 0) + const size_t local_size = local_n0 * dims[1] * dims[2]; + + // The phases are 0 ... n_mat - 1 + int highest = *std::max_element(ms, ms + local_size); + MPI_Allreduce(MPI_IN_PLACE, &highest, 1, MPI_INT, MPI_MAX, communicator); + n_mat = highest + 1; + + // Voxels of each phase, summed over all ranks + std::vector voxels(n_mat, 0); + for (size_t i = 0; i < local_size; i++) + voxels[ms[i]]++; + MPI_Allreduce(MPI_IN_PLACE, voxels.data(), n_mat, MPI_LONG, MPI_SUM, communicator); + if (voxels[0] == 0) throw std::invalid_argument("Microstructure phase IDs must start at 0"); - // Calculate total number of materials - n_mat = global_max - global_min + 1; - - Log::logger().info("# Number of materials: {} (from {} to {})", n_mat, global_min, global_max); + Log::logger().info("# Number of materials: {} (from 0 to {})", n_mat, highest); Log::logger().info("# Volume fractions"); - - // Voxels of each phase, summed over all ranks at once - std::vector vol_frac(n_mat, 0); - for (size_t i = 0; i < local_size; i++) { - unsigned short val = static_cast(ms[i]); - int index = val - global_min; // Adjust index to start from 0 - vol_frac[index]++; - } - MPI_Allreduce(MPI_IN_PLACE, vol_frac.data(), n_mat, MPI_LONG, MPI_SUM, communicator); - for (int i = 0; i < n_mat; i++) - Log::logger().info("# material {:4} vol. frac. {:10.4f}% ", static_cast(i) + global_min, - 100. * vol_frac[i] / (double(dims[0]) * dims[1] * dims[2])); + Log::logger().info("# material {:4} vol. frac. {:10.4f}% ", i, 100. * voxels[i] / (double(dims[0]) * dims[1] * dims[2])); } void Reader ::ReadInputFile(const std::string &input_fn) @@ -71,23 +42,17 @@ void Reader ::ReadInputFile(const std::string &input_fn) inputJson = j; // Store complete input JSON for MaterialManager microstructure = j["microstructure"]; - std::snprintf(ms_filename, sizeof(ms_filename), "%s", microstructure["filepath"].get().c_str()); - // dataset name handling - const auto tmp_str = microstructure["datasetname"].get(); - if (tmp_str.empty()) + ms_filename = microstructure["filepath"].get(); + ms_datasetname = microstructure["datasetname"].get(); + if (ms_datasetname.empty()) throw std::invalid_argument("datasetname must not be empty and must refer to a valid HDF5 path"); - // Ensure absolute HDF5 path, leading slash - std::snprintf(ms_datasetname, sizeof(ms_datasetname), "%s%s", tmp_str.front() == '/' ? "" : "/", tmp_str.c_str()); + if (ms_datasetname.front() != '/') + ms_datasetname = "/" + ms_datasetname; // an absolute HDF5 path L = microstructure["L"].get>(); - if (j.contains("results_prefix")) { - std::snprintf(results_prefix, sizeof(results_prefix), "%s", j["results_prefix"].get().c_str()); - } else { - strcpy(results_prefix, ""); - } - - // Construct dataset_name as "_results/" - std::snprintf(dataset_name, sizeof(dataset_name), "%s_results/%s", ms_datasetname, results_prefix); + // The group of the microstructure, and where the results go: "_results/" + ms_group = ms_datasetname.substr(0, ms_datasetname.find_last_of('/') + 1); + dataset_name = ms_datasetname + "_results/" + j.value("results_prefix", ""); errorParameters = j["error_parameters"]; TOL = errorParameters["tolerance"].get(); @@ -105,25 +70,13 @@ void Reader ::ReadInputFile(const std::string &input_fn) problemType = j["problem_type"].get(); method = j["method"].get(); - // Parse strain_type (optional, defaults to "small") - if (j.contains("strain_type")) { - strain_type = j["strain_type"].get(); - if (strain_type != "small" && strain_type != "large") { - throw std::invalid_argument("strain_type must be either 'small' or 'large'"); - } - } else { - strain_type = "small"; // Default to small strain - } + strain_type = j.value("strain_type", "small"); + if (strain_type != "small" && strain_type != "large") + throw std::invalid_argument("strain_type must be either 'small' or 'large'"); - // Parse FE_type (optional, defaults to "HEX8") - if (j.contains("FE_type")) { - FE_type = j["FE_type"].get(); - if (FE_type != "HEX8" && FE_type != "HEX8R" && FE_type != "BBAR") { - throw std::invalid_argument("FE_type must be one of: 'HEX8', 'HEX8R', or 'BBAR'"); - } - } else { - FE_type = "HEX8"; // Default to full integration - } + FE_type = j.value("FE_type", "HEX8"); // full integration + if (FE_type != "HEX8" && FE_type != "HEX8R" && FE_type != "BBAR") + throw std::invalid_argument("FE_type must be one of: 'HEX8', 'HEX8R', or 'BBAR'"); resultsToWrite = j["results"].get>(); // Read the results_to_write field @@ -203,7 +156,7 @@ void Reader ::ReadMS(int hm) hid_t space = H5Dget_space(dset_id); hsize_t file_dims[3]; if (H5Sget_simple_extent_ndims(space) != 3) - throw std::runtime_error("Microstructure dataset " + string(ms_datasetname) + " must be 3D"); + throw std::runtime_error("Microstructure dataset " + ms_datasetname + " must be 3D"); H5Sget_simple_extent_dims(space, file_dims, nullptr); H5Sclose(space); @@ -282,12 +235,6 @@ bool h5_is_zyx(hid_t dset_id) return permute_order.empty() || permute_order[0] == 'z' || permute_order[0] == 'Z'; } -string Reader::MSGroup() const -{ - const string path(ms_datasetname); - return path.substr(0, path.find_last_of('/') + 1); -} - std::vector Reader::DataShape(const string &file, const string &dset_name) { hid_t file_id; @@ -311,10 +258,9 @@ void Reader::FreeMS() void Reader::OpenResultsFile(const char *output_fn) { - std::snprintf(results_filename, sizeof(results_filename), "%s", output_fn); hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS); H5Pset_fapl_mpio(plist_id, communicator, MPI_INFO_NULL); - results_file_id = H5Fcreate(results_filename, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id); + results_file_id = H5Fcreate(output_fn, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id); H5Pclose(plist_id); if (results_file_id < 0) { @@ -328,7 +274,6 @@ void Reader::CloseResultsFile() H5Fclose(results_file_id); results_file_id = -1; } - results_filename[0] = '\0'; } Reader::~Reader() From 6412ea328c09fb79ce9237048e260182b0374b23 Mon Sep 17 00:00:00 2001 From: Sanath Date: Sat, 10 Oct 2026 10:27:45 +0200 Subject: [PATCH 7/8] compact displacement computation --- include/solver.h | 113 +++++++++++++++-------------------------------- 1 file changed, 35 insertions(+), 78 deletions(-) diff --git a/include/solver.h b/include/solver.h index a779714..df2fe36 100644 --- a/include/solver.h +++ b/include/solver.h @@ -486,11 +486,11 @@ void Solver::postprocess(Reader &reader, int load_idx, int time_ // Determine if we need to compute stress/strain at all bool need_compute = need_stress || need_stress_gp || need_strain || need_strain_gp || need_global_avg || need_phase_avg; - // Conditional allocation - VectorXd *strain_elem = need_strain ? new VectorXd(local_n0 * n_y * n_z * n_str) : nullptr; - VectorXd *stress_elem = need_stress ? new VectorXd(local_n0 * n_y * n_z * n_str) : nullptr; - VectorXd *strain_gp = need_strain_gp ? new VectorXd(local_n0 * n_y * n_z * n_gp * n_str) : nullptr; - VectorXd *stress_gp = need_stress_gp ? new VectorXd(local_n0 * n_y * n_z * n_gp * n_str) : nullptr; + // The requested fields; the others stay empty + VectorXd strain_elem(need_strain ? local_n0 * n_y * n_z * n_str : 0); + VectorXd stress_elem(need_stress ? local_n0 * n_y * n_z * n_str : 0); + VectorXd strain_gp(need_strain_gp ? local_n0 * n_y * n_z * n_gp * n_str : 0); + VectorXd stress_gp(need_stress_gp ? local_n0 * n_y * n_z * n_gp * n_str : 0); VectorXd stress_average = VectorXd::Zero(n_str); VectorXd strain_average = VectorXd::Zero(n_str); @@ -527,23 +527,23 @@ void Solver::postprocess(Reader &reader, int load_idx, int time_ // Store element averages if requested if (need_stress) { for (int c = 0; c < n_str; ++c) { - (*stress_elem)(idx[0] * n_str + c) = elem_stress_avg[c]; + stress_elem(idx[0] * n_str + c) = elem_stress_avg[c]; } } if (need_strain) { for (int c = 0; c < n_str; ++c) { - (*strain_elem)(idx[0] * n_str + c) = elem_strain_avg[c]; + strain_elem(idx[0] * n_str + c) = elem_strain_avg[c]; } } // Copy all GP data if requested if (need_stress_gp) { - memcpy(&(*stress_gp)(idx[0] * n_gp * n_str), + memcpy(&stress_gp(idx[0] * n_gp * n_str), info.model->get_sigma_data(), n_gp * n_str * sizeof(double)); } if (need_strain_gp) { - memcpy(&(*strain_gp)(idx[0] * n_gp * n_str), + memcpy(&strain_gp(idx[0] * n_gp * n_str), info.model->get_eps_data(), n_gp * n_str * sizeof(double)); } @@ -589,65 +589,32 @@ void Solver::postprocess(Reader &reader, int load_idx, int time_ homogenized_stress = stress_average; homogenized_strain = strain_average; - /* ====================================================================== * - * u_total = g0·X + ũ (vector or scalar, decided at compile time) - * ====================================================================== */ - const vector &ml = matmanager->models[0]->macroscale_loading; - const double dx = reader.l_e[0]; - const double dy = reader.l_e[1]; - const double dz = reader.l_e[2]; - const double Lx2 = reader.L[0] / 2.0; - const double Ly2 = reader.L[1] / 2.0; - const double Lz2 = reader.L[2] / 2.0; - constexpr double rs2 = 0.7071067811865475; // 1.0 / std::sqrt(2.0) - VectorXd u_total(local_n0 * n_y * n_z * howmany); - /* ---------- single sweep ------------------------------------------------- */ + // u_total = u + G X, with G the macroscale gradient of the nodal values: the temperature + // gradient, the strain (ml holds its Mandel components) or F - I + const vector &ml = matmanager->models[0]->macroscale_loading; + constexpr double rs2 = 0.7071067811865475; // 1.0 / std::sqrt(2.0) + Matrix G; + if constexpr (n_str == 3) + G << ml[0], ml[1], ml[2]; + else if constexpr (n_str == 6) + G << ml[0], ml[3] * rs2, ml[4] * rs2, + ml[3] * rs2, ml[1], ml[5] * rs2, + ml[4] * rs2, ml[5] * rs2, ml[2]; + else + G << ml[0] - 1.0, ml[1], ml[2], + ml[3], ml[4] - 1.0, ml[5], + ml[6], ml[7], ml[8] - 1.0; + + VectorXd u_total(local_n0 * n_y * n_z * howmany); ptrdiff_t n = 0; for (ptrdiff_t ix = 0; ix < local_n0; ++ix) { - const double x = (local_0_start + ix) * dx - Lx2; + const double x = (local_0_start + ix) * reader.l_e[0] - reader.L[0] / 2.0; for (ptrdiff_t iy = 0; iy < n_y; ++iy) { - const double y = iy * dy - Ly2; + const double y = iy * reader.l_e[1] - reader.L[1] / 2.0; for (ptrdiff_t iz = 0; iz < n_z; ++iz, ++n) { - const double z = iz * dz - Lz2; - if (howmany == 3) { /* ===== mechanics (vector) ===== */ - if constexpr (n_str == 6) { - /* Small strain: ml holds 6 Mandel components of the mean strain */ - const double g11 = ml[0]; - const double g22 = ml[1]; - const double g33 = ml[2]; - const double g12 = ml[3] * rs2; - const double g13 = ml[4] * rs2; - const double g23 = ml[5] * rs2; - const double ux = g11 * x + g12 * y + g13 * z; - const double uy = g12 * x + g22 * y + g23 * z; - const double uz = g13 * x + g23 * y + g33 * z; - const ptrdiff_t b = 3 * n; - u_total[b] = v_u[b] + ux; - u_total[b + 1] = v_u[b + 1] + uy; - u_total[b + 2] = v_u[b + 2] + uz; - } else if constexpr (n_str == 9) { - /* Large strain: ml holds 9 components of the mean deformation gradient F */ - const double F11 = ml[0]; - const double F12 = ml[1]; - const double F13 = ml[2]; - const double F21 = ml[3]; - const double F22 = ml[4]; - const double F23 = ml[5]; - const double F31 = ml[6]; - const double F32 = ml[7]; - const double F33 = ml[8]; - // u = (F - I) * X - const double ux = (F11 - 1.0) * x + F12 * y + F13 * z; - const double uy = F21 * x + (F22 - 1.0) * y + F23 * z; - const double uz = F31 * x + F32 * y + (F33 - 1.0) * z; - const ptrdiff_t b = 3 * n; - u_total[b] = v_u[b] + ux; - u_total[b + 1] = v_u[b + 1] + uy; - u_total[b + 2] = v_u[b + 2] + uz; - } - } else { /* ===== scalar (howmany==1) ==== */ - u_total[n] = v_u[n] + (ml[0] * x + ml[1] * y + ml[2] * z); - } + const double z = iz * reader.l_e[2] - reader.L[2] / 2.0; + for (int i = 0; i < howmany; ++i) + u_total[howmany * n + i] = v_u[howmany * n + i] + (G(i, 0) * x + G(i, 1) * y + G(i, 2) * z); } } } @@ -674,26 +641,16 @@ void Solver::postprocess(Reader &reader, int load_idx, int time_ reader.writeSlab("residual", load_idx, time_idx, v_r, {howmany}); if (need_strain) - reader.writeSlab("strain", load_idx, time_idx, strain_elem->data(), {n_str}); + reader.writeSlab("strain", load_idx, time_idx, strain_elem.data(), {n_str}); if (need_stress) - reader.writeSlab("stress", load_idx, time_idx, stress_elem->data(), {n_str}); + reader.writeSlab("stress", load_idx, time_idx, stress_elem.data(), {n_str}); if (need_strain_gp) - reader.writeSlab("strain_gp", load_idx, time_idx, strain_gp->data(), {n_gp, n_str}); + reader.writeSlab("strain_gp", load_idx, time_idx, strain_gp.data(), {n_gp, n_str}); if (need_stress_gp) - reader.writeSlab("stress_gp", load_idx, time_idx, stress_gp->data(), {n_gp, n_str}); + reader.writeSlab("stress_gp", load_idx, time_idx, stress_gp.data(), {n_gp, n_str}); matmanager->postprocess(*this, reader, load_idx, time_idx); - // Cleanup - if (strain_elem) - delete strain_elem; - if (stress_elem) - delete stress_elem; - if (strain_gp) - delete strain_gp; - if (stress_gp) - delete stress_gp; - // Compute homogenized tangent only if requested if (find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), "homogenized_tangent") != reader.resultsToWrite.end()) { homogenized_tangent = get_homogenized_tangent(1e-6); From 46ce0f6bf57ee32505d5f479e677d4748682dc33 Mon Sep 17 00:00:00 2001 From: Sanath Date: Sat, 10 Oct 2026 11:02:51 +0200 Subject: [PATCH 8/8] add changelog entry --- CHANGELOG.md | 1 + docs/NEML2.md | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 64ea413..4f0d076 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## latest +- Add consistent homogenized tangent from the materials' own tangents (NEML2 and native models) [#163](https://github.com/DataAnalyticsEngineering/FANS/pull/163) - Add NEML2 material support [#162](https://github.com/DataAnalyticsEngineering/FANS/pull/162) - Account for grain orientations in GBDiffusion material model [#161](https://github.com/DataAnalyticsEngineering/FANS/pull/161) diff --git a/docs/NEML2.md b/docs/NEML2.md index c5ead69..b4169b1 100644 --- a/docs/NEML2.md +++ b/docs/NEML2.md @@ -30,7 +30,7 @@ This gives `test/FANS_neml2`. Without pixi, configure with `-DFANS_NEML2=ON` in pixi run -e dev-neml2 neml2-compile model.i --model model --device cpu cuda --dtype float64 --output-dir compiled_models -d state/S:forces/E ``` -This compiles the block `[model]` of `model.i` for both devices into `compiled_models/model`. `--load a.py` (repeatable) imports a Python file the model needs. `-d flux:gradient` also compiles the model's tangent. +This compiles the block `[model]` of `model.i` for both devices into `compiled_models/model`. `--load a.py` (repeatable) imports a Python file the model needs. `-d flux:gradient` also compiles the model's tangent, which the `homogenized_tangent` result needs. ## Input file