Added support for large strain in pyFANS

CHANGELOG.md

@@ -6,6 +6,7 @@
 
 ## v0.6.1
 
+- Add pyFANS support for large strain, fixed pyFANS MPI, added preliminary pyFANS state setting&loading and added logging [#124](https://github.com/DataAnalyticsEngineering/FANS/pull/124)
 - Adds support to write integration point data to disk [#117](https://github.com/DataAnalyticsEngineering/FANS/pull/117)
 
 ## v0.6.0

CMakeLists.txt

@@ -112,6 +112,9 @@ find_package(nlohmann_json REQUIRED)
 # TARGETS
 # ##############################################################################
 
+set(FANS_VERBOSITY 5 CACHE STRING "Set verbosity level: 0 none, 5 all")
+add_definitions(-DVERBOSITY=${FANS_VERBOSITY})
+
 option(FANS_BUILD_STATIC "Build static library" OFF)
 if (FANS_BUILD_STATIC)
     add_library(FANS_FANS STATIC)
@@ -150,6 +153,7 @@ target_include_directories(FANS_FANS PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SO
 # FANS.h header that consuming projects can use
 set_property(TARGET FANS_FANS PROPERTY PUBLIC_HEADER
         include/general.h
+        include/logging.h
         include/matmodel.h
         include/MaterialManager.h
         include/reader.h
@@ -184,6 +188,7 @@ target_include_directories(FANS_main PRIVATE "${PROJECT_BINARY_DIR}/include")
 
 target_sources(FANS_FANS PRIVATE
         src/reader.cpp
+        src/logging.cpp
 )
 
 target_sources(FANS_main PRIVATE
@@ -216,7 +221,7 @@ target_compile_definitions(FANS_FANS PUBLIC ${FFTW3_DEFINITIONS})
 target_link_libraries(FANS_FANS PUBLIC Eigen3::Eigen)
 target_link_libraries(FANS_FANS PUBLIC nlohmann_json::nlohmann_json)
 
-target_include_directories(FANS_main PRIVATE ${HDF5_INCLUDE_DIRS})
+target_include_directories(FANS_main PRIVATE ${HDF5_INCLUDE_DIRS} ${FFTW3_INCLUDE_DIRS})
 target_link_libraries(FANS_main PRIVATE FANS::FANS)
 
 # ##############################################################################

cmake/modules/FindFFTW3.cmake

@@ -72,13 +72,25 @@ macro(find_specific_libraries KIND PARALLEL)
     # adding target properties to the different cases
     ##   MPI
     if(PARALLEL STREQUAL "MPI")
-      if(MPI_C_LIBRARIES)
-        set_target_properties(fftw3::${kind}::mpi PROPERTIES
-          IMPORTED_LOCATION "${FFTW3_${KIND}_${PARALLEL}_LIBRARY}"
-          INTERFACE_INCLUDE_DIRECTORIES "${FFTW3_INCLUDE_DIR_PARALLEL}"
-          IMPORTED_LINK_INTERFACE_LIBRARIES ${MPI_C_LIBRARIES})
+      set_target_properties(fftw3::${kind}::mpi PROPERTIES
+              IMPORTED_LOCATION "${FFTW3_${KIND}_${PARALLEL}_LIBRARY}"
+              INTERFACE_INCLUDE_DIRECTORIES "${FFTW3_INCLUDE_DIR_PARALLEL}"
+      )
+
+      if(MPI_C_LIBRARIES AND NOT MPI_C_LIBRARIES STREQUAL "")
+        set_property(TARGET fftw3::${kind}::mpi PROPERTY
+                INTERFACE_LINK_LIBRARIES "${MPI_C_LIBRARIES}"
+        )
       endif()
     endif()
+    #if(PARALLEL STREQUAL "MPI")
+    #  if(MPI_C_LIBRARIES)
+    #    set_target_properties(fftw3::${kind}::mpi PROPERTIES
+    #      IMPORTED_LOCATION "${FFTW3_${KIND}_${PARALLEL}_LIBRARY}"
+    #      INTERFACE_INCLUDE_DIRECTORIES "${FFTW3_INCLUDE_DIR_PARALLEL}"
+    #      IMPORTED_LINK_INTERFACE_LIBRARIES ${MPI_C_LIBRARIES})
+    #  endif()
+    #endif()
     ##   OpenMP
     if(PARALLEL STREQUAL "OPENMP")
       if(OPENMP_C_FLAGS)

include/MaterialManager.h

@@ -2,6 +2,7 @@
 #define MATERIALMANAGER_H
 
 #include "general.h"
+#include "logging.h"
 #include "matmodel.h"
 
 template <int howmany, int n_str>
@@ -110,51 +111,49 @@ class MaterialManager {
         compute_reference_stiffness(reader);
 
         // Print detailed information about material configuration for logging
-        if (reader.world_rank == 0) {
-            printf("\n# MaterialManager initialized:\n");
-            printf("#   Number of material models: %zu\n", models.size());
-            printf("#   Number of phases: %d\n#\n", n_phases);
-
-            for (size_t i = 0; i < mats.size(); ++i) {
-                const auto &mg = mats[i];
-                printf("# Material Model %zu: %s\n", i + 1, mg["matmodel"].get<string>().c_str());
-
-                // Print phases
-                auto phases = mg["phases"].get<vector<int>>();
-                printf("#   Phases: [");
-                for (size_t j = 0; j < phases.size(); ++j) {
-                    printf("%d", phases[j]);
-                    if (j < phases.size() - 1)
-                        printf(", ");
-                }
-                printf("]\n");
-
-                // Print material properties
-                printf("#   Material properties:\n");
-                const auto &props = mg["material_properties"];
-                for (auto it = props.begin(); it != props.end(); ++it) {
-                    printf("#     %s: ", it.key().c_str());
-                    if (it.value().is_array()) {
-                        printf("[");
-                        for (size_t k = 0; k < it.value().size(); ++k) {
-                            if (it.value()[k].is_number()) {
-                                printf("%.5g", it.value()[k].get<double>());
-                            } else if (it.value()[k].is_string()) {
-                                printf("\"%s\"", it.value()[k].get<string>().c_str());
-                            }
-                            if (k < it.value().size() - 1)
-                                printf(", ");
+        Log::io->info() << "\n# MaterialManager initialized:\n";
+        Log::io->info() << "#   Number of material models: " << models.size() << "\n";
+        Log::io->info() << "#   Number of phases: " << n_phases << "\n#\n";
+
+        for (size_t i = 0; i < mats.size(); ++i) {
+            const auto &mg = mats[i];
+            Log::io->info() << "# Material Model " << i + 1 << ": " << mg["matmodel"].get<string>() << "\n";
+
+            // Print phases
+            auto phases = mg["phases"].get<vector<int>>();
+            Log::io->info() << "#   Phases: [";
+            for (size_t j = 0; j < phases.size(); ++j) {
+                Log::io->info(true) << phases[j];
+                if (j < phases.size() - 1)
+                    Log::io->info(true) << ", ";
+            }
+            Log::io->info(true) << "]\n";
+
+            // Print material properties
+            Log::io->info() << "#   Material properties:\n";
+            const auto &props = mg["material_properties"];
+            for (auto it = props.begin(); it != props.end(); ++it) {
+                Log::io->info() << "#     " << it.key() << ": ";
+                if (it.value().is_array()) {
+                    Log::io->info(true) << "[";
+                    for (size_t k = 0; k < it.value().size(); ++k) {
+                        if (it.value()[k].is_number()) {
+                            Log::io->info(true) << std::setprecision(5) << it.value()[k].get<double>() << std::defaultfloat;
+                        } else if (it.value()[k].is_string()) {
+                            Log::io->info(true) << "\"" << it.value()[k].get<string>() << "\"";
                         }
-                        printf("]");
-                    } else if (it.value().is_number()) {
-                        printf("%.5g", it.value().get<double>());
-                    } else if (it.value().is_string()) {
-                        printf("\"%s\"", it.value().get<string>().c_str());
+                        if (k < it.value().size() - 1)
+                            Log::io->info(true) << ", ";
                     }
-                    printf("\n");
+                    Log::io->info(true) << "]";
+                } else if (it.value().is_number()) {
+                    Log::io->info(true) << std::setprecision(5) << it.value().get<double>() << std::defaultfloat;
+                } else if (it.value().is_string()) {
+                    Log::io->info(true) << "\"" << it.value().get<string>() << "\"";
                 }
-                printf("#\n");
+                Log::io->info(true) << "\n";
             }
+            Log::io->info() << "#\n";
         }
     }
 
@@ -213,7 +212,7 @@ class MaterialManager {
     }
 
     // Set macroscale loading gradient for all models
-    void set_gradient(vector<double> g0)
+    void set_gradient(const vector<double> &g0)
     {
         for (auto *model : models) {
             model->setGradient(g0);

include/general.h

@@ -58,6 +58,4 @@ inline void FANS_free(V *p)
 }
 #endif // FANS_MALLOC_H
 
-#define VERBOSITY 0
-
 // #define EIGEN_RUNTIME_NO_MALLOC

include/logging.h

@@ -0,0 +1,106 @@
+#ifndef LOGGING_H
+#define LOGGING_H
+
+#include <chrono>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <memory>
+#include "mpi.h"
+
+namespace Log {
+
+enum Level {
+    Error,
+    Info,
+    Warn,
+    Debug,
+    All
+};
+
+[[maybe_unused]] const std::string level_to_string(Level level);
+
+extern Level active_level;
+
+class Logger;
+class MPI_TraceSync {
+  public:
+    explicit MPI_TraceSync(Logger &log, bool append);
+    ~MPI_TraceSync();
+
+    template <class T>
+    MPI_TraceSync &operator<<(const T &v);
+    MPI_TraceSync &operator<<(std::ostream &(*m)(std::ostream &) );
+
+  private:
+    Logger            &_log;
+    bool               _append;
+    std::ostringstream _buffer;
+};
+
+template <class T>
+MPI_TraceSync &MPI_TraceSync::operator<<(const T &v)
+{
+    _buffer << v;
+    return *this;
+}
+
+class Logger {
+  public:
+    friend class MPI_TraceSync;
+    explicit Logger(std::string prefix, int comm_rank, int comm_size, const MPI_Comm &comm);
+
+    /// error > info > warn > debug > trace
+    std::ostream &error(bool append = false);
+    /// error > info > warn > debug > trace
+    std::ostream &info(bool append = false);
+    /// error > info > warn > debug > trace
+    std::ostream &warn(bool append = false);
+    /// error > info > warn > debug > trace
+    std::ostream &debug(bool append = false);
+    /// error > info > warn > debug > trace
+    MPI_TraceSync trace(bool append = false);
+    /// progress bar
+    void progress(const std::string &prefix, int step, int max);
+
+  private:
+    std::ostream &trace_impl(bool append = false);
+    /// starting time
+    std::chrono::steady_clock::time_point _start_time;
+    /// what the logger should always print first
+    const std::string _prefix;
+    /// empty stream to write nothing
+    std::ofstream _nullstr;
+    /// MPI comm rank
+    int _comm_rank;
+    /// MPI comm size
+    int _comm_size;
+    /// communicator
+    MPI_Comm _comm;
+};
+
+/**
+ * Creates all loggers and sets the level
+ * */
+void init(int comm_rank, int comm_size, const MPI_Comm &comm);
+
+/**
+ * Frees all memory
+ * */
+void finalize();
+
+/**
+ * Set activate rank
+ */
+[[maybe_unused]] void setActiveRank(int rank);
+
+/// logger with prefix [GENERAL]
+extern std::unique_ptr<Logger> general;
+/// logger with prefix [SOLVER]
+extern std::unique_ptr<Logger> solver;
+/// logger with prefix [IO]
+extern std::unique_ptr<Logger> io;
+
+}; // namespace Log
+
+#endif

include/matmodel.h

@@ -27,7 +27,7 @@ class Matmodel {
     Matrix<double, howmany * 8, howmany * 8> Compute_Reference_ElementStiffness(const Matrix<double, n_str, n_str> &kapparef_mat);
     Matrix<double, howmany * 8, 1>          &element_residual(Matrix<double, howmany * 8, 1> &ue, int mat_index, ptrdiff_t element_idx);
     void                                     getStrainStress(double *strain, double *stress, Matrix<double, howmany * 8, 1> &ue, int mat_index, ptrdiff_t element_idx);
-    void                                     setGradient(vector<double> _g0);
+    void                                     setGradient(const vector<double> &_g0);
 
     // Accessors for internal Gauss point data (populated after getStrainStress call)
     inline const double *get_eps_data() const
@@ -225,11 +225,13 @@ void Matmodel<howmany, n_str>::getStrainStress(double *strain, double *stress, M
 }
 
 template <int howmany, int n_str>
-void Matmodel<howmany, n_str>::setGradient(vector<double> _g0)
+void Matmodel<howmany, n_str>::setGradient(const vector<double> &_g0)
 {
+    // _g0 is smaller than g0, can stay in cache
+    // if j loop inner, then not using cache-lines well
     macroscale_loading = _g0;
-    for (int i = 0; i < n_str; i++) {
-        for (int j = 0; j < n_gp; ++j) {
+    for (int j = 0; j < n_gp; ++j) {
+        for (int i = 0; i < n_str; i++) {
             g0(n_str * j + i, 0) = _g0[i];
         }
     }
@@ -307,6 +309,11 @@ template <int howmany, int n_str>
 class LinearModel {
   public:
     Matrix<double, howmany * 8, howmany * 8> *phase_stiffness;
+
+    ~LinearModel()
+    {
+        delete[] phase_stiffness;
+    }
 };
 
 #endif // MATMODEL_H

include/reader.h

@@ -15,12 +15,13 @@ using namespace std;
 
 class Reader {
   public:
-    // Default constructor
-    Reader() = default;
-
     // Destructor to free allocated memory
     ~Reader();
 
+    // MPI controlls
+    bool     force_single_rank;
+    MPI_Comm communicator;
+
     // 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
@@ -64,7 +65,9 @@ class Reader {
 
     // void Setup(ptrdiff_t howmany);
     void ReadInputFile(char input_fn[]);
+    void ReadJson(json j);
     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);