From 49d7a9fa39dce753793a58831eacd797496373c5 Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Thu, 14 Sep 2023 13:46:01 +0200 Subject: [PATCH 1/3] addresses Carsten review --- app/form/SeasQDDiscreteGreenOperator.cpp | 277 +++++++++++------------ app/form/SeasQDDiscreteGreenOperator.h | 19 +- app/tandem/SEAS.cpp | 32 +-- 3 files changed, 141 insertions(+), 187 deletions(-) diff --git a/app/form/SeasQDDiscreteGreenOperator.cpp b/app/form/SeasQDDiscreteGreenOperator.cpp index 6f3a3bd2..f5140462 100644 --- a/app/form/SeasQDDiscreteGreenOperator.cpp +++ b/app/form/SeasQDDiscreteGreenOperator.cpp @@ -6,6 +6,7 @@ #include "util/Stopwatch.h" #include +#include namespace fs = std::filesystem; namespace tndm { @@ -14,24 +15,40 @@ SeasQDDiscreteGreenOperator::SeasQDDiscreteGreenOperator( std::unique_ptr dgop, std::unique_ptr adapter, std::unique_ptr friction, LocalSimplexMesh const& mesh, - bool matrix_free, MGConfig const& mg_config, std::string prefix) + std::optional prefix, + bool matrix_free, MGConfig const& mg_config) : base(std::move(dgop), std::move(adapter), std::move(friction), matrix_free, mg_config) { int rank; MPI_Comm_rank(base::comm(),&rank); // if prefix is not empty, set filenames and mark checkpoint_enabled_ = true - if (!prefix.empty()) { - if (rank == 0) { std::cout << "Using GF checkpoint path: " << prefix << std::endl; } - fs::path pckp(prefix); + if (prefix.has_value()) { + std::string sprefix = prefix.value_or(""); + if (rank == 0) { std::cout << "Using GF checkpoint path: " << sprefix << std::endl; } + auto pckp = fs::path(sprefix); if (rank == 0) { bool exists = fs::exists(pckp); if (!exists) { bool ret = fs::create_directories(pckp); - if (!ret) std::cout << "--> Failed to create directory!" << std::endl; + if (!ret) throw std::runtime_error("Failed to create directory!"); } } - prepend_checkpoint_path(prefix); + auto prepend_checkpoint_path = [&](std::string pre_path) { + fs::path pckpOp(pre_path); + pckpOp /= gf_operator_filename_; + gf_operator_filename_ = pckpOp; + + fs::path pckpVec(pre_path); + pckpVec /= gf_traction_filename_; + gf_traction_filename_ = pckpVec; + + fs::path pckpFacet(pre_path); + pckpFacet /= gf_facet_filename_; + gf_facet_filename_ = pckpFacet; + }; + + prepend_checkpoint_path(sprefix); checkpoint_enabled_ = true; } if (!checkpoint_enabled_) { @@ -53,36 +70,7 @@ void SeasQDDiscreteGreenOperator::set_boundary( } } - std::tuple SeasQDDiscreteGreenOperator::get_checkpoint_filenames(void) { - return std::make_tuple(gf_operator_filename_, gf_traction_filename_); - } - double SeasQDDiscreteGreenOperator::get_checkpoint_time_interval(void) { - return checkpoint_every_nmins_; - } - - void SeasQDDiscreteGreenOperator::set_checkpoint_filenames(std::string mat_fname, std::string vec_fname) { - gf_operator_filename_ = mat_fname; - gf_traction_filename_ = vec_fname; - } - - void SeasQDDiscreteGreenOperator::prepend_checkpoint_path(std::string pre_path) { - fs::path pckpOp(pre_path); - pckpOp /= gf_operator_filename_; - gf_operator_filename_ = pckpOp; - - fs::path pckpVec(pre_path); - pckpVec /= gf_traction_filename_; - gf_traction_filename_ = pckpVec; - - fs::path pckpFacet(pre_path); - pckpFacet /= gf_facet_filename_; - gf_facet_filename_ = pckpFacet; - } - - void SeasQDDiscreteGreenOperator::set_checkpoint_time_interval(double t) { - checkpoint_every_nmins_ = t; - } void SeasQDDiscreteGreenOperator::update_internal_state(double time, BlockVector const& state, @@ -229,9 +217,9 @@ void SeasQDDiscreteGreenOperator::compute_boundary_traction() { CHKERRTHROW(VecCopy(base::traction_.vec(), t_boundary_->vec())); } -void SeasQDDiscreteGreenOperator::create_discrete_greens_function() { +PetscInt SeasQDDiscreteGreenOperator::create_discrete_greens_function() { auto slip_block_size = base::friction().slip_block_size(); - PetscInt M,N; + PetscInt M,N,n_gf; PetscInt num_local_elements = base::adapter().num_local_elements(); PetscInt m_bs = base::adapter().traction_block_size(); PetscInt n_bs = 1; @@ -257,11 +245,11 @@ void SeasQDDiscreteGreenOperator::create_discrete_greens_function() { S_ = std::make_unique(slip_block_size, num_local_elements, comm); t_boundary_ = std::make_unique(m_bs, num_local_elements, comm); - current_gf_ = 0; - CHKERRTHROW(VecGetSize(S_->vec(), &n_gf_)); + CHKERRTHROW(VecGetSize(S_->vec(), &n_gf)); + return n_gf; } -void SeasQDDiscreteGreenOperator::write_discrete_greens_operator(LocalSimplexMesh const& mesh) { +void SeasQDDiscreteGreenOperator::write_discrete_greens_operator(LocalSimplexMesh const& mesh, PetscInt current_gf, PetscInt n_gf) { PetscViewer v; PetscLogDouble t0,t1; int commsize; @@ -272,18 +260,18 @@ void SeasQDDiscreteGreenOperator::write_discrete_greens_operator(LocalSimplexMes CHKERRTHROW(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)G_),gf_operator_filename_.c_str(),FILE_MODE_WRITE,&v)); { - PetscInt _commsize = (PetscInt)commsize; - CHKERRTHROW(PetscViewerBinaryWrite(v,&_commsize,1,PETSC_INT)); + PetscInt commsize_checkpoint = (PetscInt)commsize; + CHKERRTHROW(PetscViewerBinaryWrite(v,&commsize_checkpoint,1,PETSC_INT)); } - CHKERRTHROW(PetscViewerBinaryWrite(v,¤t_gf_,1,PETSC_INT)); + CHKERRTHROW(PetscViewerBinaryWrite(v,¤t_gf,1,PETSC_INT)); CHKERRTHROW(MatView(G_,v)); CHKERRTHROW(PetscViewerDestroy(&v)); CHKERRTHROW(PetscTime(&t1)); CHKERRTHROW(PetscPrintf(PetscObjectComm((PetscObject)G_),"write_discrete_greens_operator():matrix %1.2e (sec)\n",(double)(t1 - t0))); - CHKERRTHROW(PetscPrintf(PetscObjectComm((PetscObject)G_)," status: computed %D / pending %D\n",current_gf_, n_gf_ - current_gf_)); + CHKERRTHROW(PetscPrintf(PetscObjectComm((PetscObject)G_)," status: computed %D / pending %D\n",current_gf, n_gf - current_gf)); write_facet_labels_IS(mesh); } @@ -293,22 +281,25 @@ void SeasQDDiscreteGreenOperator :: write_facet_labels_IS(LocalSimplexMesh const& mesh, IS is) +void SeasQDDiscreteGreenOperator::create_permutation_redundant_IS(LocalSimplexMesh const& mesh, IS is) { MPI_Comm comm = base::comm(); - std::size_t bndNo; - PetscInt nfacets = (PetscInt)base::adapter().fault_map().local_size(); - PetscInt facet_size = (PetscInt)(DomainDimension); - PetscInt *_idx = NULL,is_len; + auto const& fault_map = base::adapter().fault_map(); + PetscInt nfacets = (PetscInt) fault_map.local_size(); + constexpr PetscInt facet_size = (PetscInt)(DomainDimension); + PetscInt *idx_ = NULL,is_len; PetscInt *fault_map_index = NULL; int rank; MPI_Comm_rank(comm,&rank); if (rank == 0) { CHKERRTHROW(ISGetSize(is,&is_len)); - CHKERRTHROW(ISGetIndices(is,(const PetscInt**)&_idx)); + CHKERRTHROW(ISGetIndices(is,(const PetscInt**)&idx_)); } MPI_Bcast((void*)&is_len,1,MPIU_INT,0,comm); - if (!_idx) { - CHKERRTHROW(PetscCalloc1(is_len,&_idx)); + if (!idx_) { + CHKERRTHROW(PetscCalloc1(is_len,&idx_)); } - MPI_Bcast((void*)_idx,is_len,MPIU_INT,0,comm); + MPI_Bcast((void*)idx_,is_len,MPIU_INT,0,comm); - CHKERRTHROW(PetscCalloc1((PetscInt)base::adapter().fault_map().local_size(),&fault_map_index)); + CHKERRTHROW(PetscCalloc1(nfacets,&fault_map_index)); - for (bndNo = 0; bndNo < base::adapter().fault_map().local_size(); ++bndNo) { - auto fctNo = base::adapter().fault_map().fctNo(bndNo); - auto facets = mesh.facets()[fctNo]; - PetscInt facet_size = (PetscInt)(DomainDimension); - PetscInt n,map_index = -1; - - for (n=0; n, std::size_t, SimplexHash>; + auto map = map_t(nfacets); + std::size_t map_index = 0; + for (PetscInt n=0; n{}; + for (PetscInt d=0; dsecond; } if (rank != 0) { - CHKERRTHROW(PetscFree(_idx)); + CHKERRTHROW(PetscFree(idx_)); } else { - CHKERRTHROW(ISRestoreIndices(is,(const PetscInt**)&_idx)); + CHKERRTHROW(ISRestoreIndices(is,(const PetscInt**)&idx_)); } { - PetscInt len = (PetscInt)base::adapter().fault_map().local_size(); - CHKERRTHROW(ISCreateGeneral(comm,len,(const PetscInt*)fault_map_index,PETSC_COPY_VALUES,&is_perm_)); + CHKERRTHROW(ISCreateGeneral(comm,nfacets,(const PetscInt*)fault_map_index,PETSC_COPY_VALUES,&is_perm_)); CHKERRTHROW(PetscFree(fault_map_index)); } } @@ -405,7 +394,6 @@ std::tuple SeasQDDiscreteGreenOperator :: create_row_col_permutation_m MPI_Comm comm = base::comm(); PetscInt num_local_elements = base::adapter().num_local_elements(); PetscInt num_global_elements; - std::size_t bndNo; PetscInt nfacets = (PetscInt)base::adapter().fault_map().local_size(); PetscInt M,N,m,n,br,bc,d; Mat Rperm = NULL,Cperm = NULL; @@ -417,8 +405,7 @@ std::tuple SeasQDDiscreteGreenOperator :: create_row_col_permutation_m MPI_Scan(&num_local_elements, &mb_offset, 1, MPIU_INT, MPI_SUM, comm); mb_offset -= num_local_elements; - num_global_elements = num_local_elements; - MPI_Allreduce(MPI_IN_PLACE,&num_global_elements,1,MPIU_INT,MPI_SUM,comm); + MPI_Allreduce(&num_local_elements,&num_global_elements,1,MPIU_INT,MPI_SUM,comm); CHKERRTHROW(MatGetSize(G_,&M,&N)); CHKERRTHROW(MatGetLocalSize(G_,&m,&n)); @@ -426,56 +413,54 @@ std::tuple SeasQDDiscreteGreenOperator :: create_row_col_permutation_m br = M/num_global_elements; bc = N/num_global_elements; + assert(M % num_global_elements == 0); + assert(N % num_global_elements == 0); + // G = Rperm G_ Cperm // Rperm is the transpose of Cperm just with a different block size - if (create_row) { - //CHKERRTHROW(MatCreateDense(comm, m, m, M, M, nullptr, &Rperm)); - CHKERRTHROW(MatCreate(comm,&Rperm)); - CHKERRTHROW(MatSetSizes(Rperm,m,m,M,M)); - CHKERRTHROW(MatSetType(Rperm,MATAIJ)); - CHKERRTHROW(MatSeqAIJSetPreallocation(Rperm,br,NULL)); - CHKERRTHROW(MatMPIAIJSetPreallocation(Rperm,br,NULL,br,NULL)); - - for (bndNo = 0; bndNo < base::adapter().fault_map().local_size(); ++bndNo) { - PetscInt from = mb_offset + bndNo; - PetscInt to = fault_map_index[bndNo]; - for (d=0; dadapter().fault_map().local_size(); ++bndNo) { PetscInt from = mb_offset + bndNo; PetscInt to = fault_map_index[bndNo]; - for (d=0; d const& mesh) { +PetscInt SeasQDDiscreteGreenOperator::load_discrete_greens_operator(LocalSimplexMesh const& mesh, PetscInt n_gf) { PetscViewer v; PetscLogDouble t0,t1; MPI_Comm comm = base::comm(); int commsize; - PetscInt _commsize; + PetscInt commsize_checkpoint; + PetscInt current_gf = 0; if (!G_) { CHKERRTHROW(PetscPrintf(base::comm(),"G_ is NULL. Cannot load the operator. Must call create_discrete_greens_function() first!\n")); @@ -484,23 +469,23 @@ void SeasQDDiscreteGreenOperator :: load_discrete_greens_operator(LocalSimplexMe CHKERRTHROW(PetscTime(&t0)); CHKERRTHROW(PetscViewerBinaryOpen(PetscObjectComm((PetscObject)G_),gf_operator_filename_.c_str(),FILE_MODE_READ,&v)); - CHKERRTHROW(PetscViewerBinaryRead(v,&_commsize,1,NULL,PETSC_INT)); + CHKERRTHROW(PetscViewerBinaryRead(v,&commsize_checkpoint,1,NULL,PETSC_INT)); MPI_Comm_size(comm,&commsize); - if ((PetscInt)commsize != _commsize) { - CHKERRTHROW(PetscPrintf(comm,"GF loaded was created with commsize = %d. Current commsize = %d. Repartitioning required.\n",(int)_commsize,(int)commsize)); + if ((PetscInt)commsize != commsize_checkpoint) { + CHKERRTHROW(PetscPrintf(comm,"GF loaded was created with commsize = %d. Current commsize = %d. Repartitioning required.\n",(int)commsize_checkpoint,(int)commsize)); repartition_gfs_ = true; } else { CHKERRTHROW(PetscPrintf(comm,"GF loaded was created with commsize matching current (%d).\n",(int)commsize)); repartition_gfs_ = false; } - CHKERRTHROW(PetscViewerBinaryRead(v,¤t_gf_,1,NULL,PETSC_INT)); + CHKERRTHROW(PetscViewerBinaryRead(v,¤t_gf,1,NULL,PETSC_INT)); CHKERRTHROW(MatLoad(G_,v)); CHKERRTHROW(PetscViewerDestroy(&v)); CHKERRTHROW(PetscTime(&t1)); CHKERRTHROW(PetscPrintf(PetscObjectComm((PetscObject)G_),"load_discrete_greens_operator() %1.2e (sec)\n",(double)(t1 - t0))); - CHKERRTHROW(PetscPrintf(PetscObjectComm((PetscObject)G_)," status: loaded %D / pending %D\n",current_gf_, n_gf_ - current_gf_)); + CHKERRTHROW(PetscPrintf(PetscObjectComm((PetscObject)G_)," status: loaded %D / pending %D\n",current_gf, n_gf - current_gf)); if (repartition_gfs_) { IS is = load_facet_labels_seq_IS(); @@ -508,17 +493,18 @@ void SeasQDDiscreteGreenOperator :: load_discrete_greens_operator(LocalSimplexMe CHKERRTHROW(ISDestroy(&is)); Mat GCperm; - auto perm_mat = create_row_col_permutation_matrices(PETSC_TRUE, PETSC_TRUE); + auto [Rperm, Cperm] = create_row_col_permutation_matrices(PETSC_TRUE, PETSC_TRUE); - CHKERRTHROW(MatMatMult(G_,std::get<1>(perm_mat) ,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&GCperm)); - CHKERRTHROW(MatMatMult(std::get<0>(perm_mat),GCperm,MAT_REUSE_MATRIX,PETSC_DEFAULT,&G_)); + CHKERRTHROW(MatMatMult(G_, Cperm ,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&GCperm)); + CHKERRTHROW(MatMatMult(Rperm,GCperm,MAT_REUSE_MATRIX,PETSC_DEFAULT,&G_)); CHKERRTHROW(MatDestroy(&GCperm)); - CHKERRTHROW(MatDestroy(&std::get<0>(perm_mat))); - CHKERRTHROW(MatDestroy(&std::get<1>(perm_mat))); + CHKERRTHROW(MatDestroy(&Rperm)); + CHKERRTHROW(MatDestroy(&Cperm)); } + return current_gf; } -void SeasQDDiscreteGreenOperator::partial_assemble_discrete_greens_function(LocalSimplexMesh const& mesh) { +void SeasQDDiscreteGreenOperator::partial_assemble_discrete_greens_function(LocalSimplexMesh const& mesh, PetscInt current_gf, PetscInt n_gf) { auto slip_block_size = base::friction().slip_block_size(); PetscInt num_local_elements = base::adapter().num_local_elements(); @@ -528,8 +514,8 @@ void SeasQDDiscreteGreenOperator::partial_assemble_discrete_greens_function(Loca PetscInt n = num_local_elements * slip_block_size * n_bs; PetscInt mb_offset = 0; PetscInt nb_offset = 0; - PetscInt start = current_gf_; - PetscInt N = n_gf_; + PetscInt start = current_gf; + PetscInt N = n_gf; Stopwatch sw; double solve_time; int rank; @@ -578,12 +564,12 @@ void SeasQDDiscreteGreenOperator::partial_assemble_discrete_greens_function(Loca base::traction_.end_access_readonly(traction_handle); solve_time += sw.stop(); - current_gf_ = i + 1; + current_gf = i + 1; /* checkpoint */ MPI_Bcast(&solve_time,1,MPI_DOUBLE,0,comm); if (solve_time/60.0 > checkpoint_every_nmins_) { - write_discrete_greens_operator(mesh); + write_discrete_greens_operator(mesh, current_gf, n_gf); solve_time = 0.0; } } @@ -594,26 +580,25 @@ void SeasQDDiscreteGreenOperator::partial_assemble_discrete_greens_function(Loca void SeasQDDiscreteGreenOperator::get_discrete_greens_function(LocalSimplexMesh const& mesh) { PetscBool found = PETSC_FALSE; - PetscInt n_gf_loaded = 0; + PetscInt n_gfloaded = 0, n_gf; // Create an empty dense matrix to store all GFs if (!G_) { - create_discrete_greens_function(); + n_gf = create_discrete_greens_function(); } // If a checkpoint file is found, load it. Record the number of assembled GFs found in file. CHKERRTHROW(PetscTestFile(gf_operator_filename_.c_str(),'r',&found)); if (found) { - load_discrete_greens_operator(mesh); - n_gf_loaded = current_gf_; + n_gfloaded = load_discrete_greens_operator(mesh, n_gf); } - // Assemble as many GFs as possible in the range [current_gf_, n_gf_) - partial_assemble_discrete_greens_function(mesh); + // Assemble as many GFs as possible in the range [current_gf, n_gf) + partial_assemble_discrete_greens_function(mesh, n_gfloaded, n_gf); // Write out the operator whenever the fully assembled operator was not loaded from file - if (n_gf_loaded != n_gf_) { - write_discrete_greens_operator(mesh); + if (n_gfloaded != n_gf) { + write_discrete_greens_operator(mesh, n_gfloaded, n_gf); } } @@ -644,9 +629,9 @@ void SeasQDDiscreteGreenOperator::load_discrete_greens_traction() { Vec tmp; CHKERRTHROW(VecDuplicate(base::traction_.vec(),&tmp)); - auto perm_mat = create_row_col_permutation_matrices(true, false); - CHKERRTHROW(MatMult(std::get<0>(perm_mat),base::traction_.vec(),tmp)); - CHKERRTHROW(MatDestroy(&std::get<0>(perm_mat))); + auto [Rperm, Cperm] = create_row_col_permutation_matrices(true, false); + CHKERRTHROW(MatMult(Rperm,base::traction_.vec(),tmp)); + CHKERRTHROW(MatDestroy(&Rperm)); CHKERRTHROW(VecCopy(tmp,base::traction_.vec())); CHKERRTHROW(VecDestroy(&tmp)); } diff --git a/app/form/SeasQDDiscreteGreenOperator.h b/app/form/SeasQDDiscreteGreenOperator.h index e0f6dc99..db8b7759 100644 --- a/app/form/SeasQDDiscreteGreenOperator.h +++ b/app/form/SeasQDDiscreteGreenOperator.h @@ -25,8 +25,8 @@ class SeasQDDiscreteGreenOperator : public SeasQDOperator { SeasQDDiscreteGreenOperator(std::unique_ptr dgop, std::unique_ptr adapter, std::unique_ptr friction, -LocalSimplexMesh const& mesh, - bool matrix_free = false, MGConfig const& mg_config = MGConfig(), std::string prefix = ""); +LocalSimplexMesh const& mesh, std::optional prefix, + bool matrix_free = false, MGConfig const& mg_config = MGConfig()); ~SeasQDDiscreteGreenOperator(); void set_boundary(std::unique_ptr fun) override; @@ -49,11 +49,6 @@ LocalSimplexMesh const& mesh, bool state_changed_since_last_rhs, bool require_traction, bool require_displacement); - std::tuple get_checkpoint_filenames(void); - double get_checkpoint_time_interval(void); - void set_checkpoint_filenames(std::string, std::string); - void set_checkpoint_time_interval(double); - void prepend_checkpoint_path(std::string); protected: std::string gf_operator_filename_ = "gf_mat.bin"; @@ -66,10 +61,10 @@ LocalSimplexMesh const& mesh, private: void compute_discrete_greens_function(); void compute_boundary_traction(); - void create_discrete_greens_function(); - void partial_assemble_discrete_greens_function(LocalSimplexMesh const& mesh); - void write_discrete_greens_operator(LocalSimplexMesh const& mesh); - void load_discrete_greens_operator(LocalSimplexMesh const& mesh); + PetscInt create_discrete_greens_function(); + void partial_assemble_discrete_greens_function(LocalSimplexMesh const& mesh, PetscInt current_gf_, PetscInt n_gf_); + void write_discrete_greens_operator(LocalSimplexMesh const& mesh, PetscInt current_gf_, PetscInt n_gf_); + PetscInt load_discrete_greens_operator(LocalSimplexMesh const& mesh, PetscInt n_gf_); // all logic associated with matix craetion, loading / partial assembly is done here void get_discrete_greens_function(LocalSimplexMesh const& mesh); void write_discrete_greens_traction(); @@ -82,8 +77,6 @@ LocalSimplexMesh const& mesh, std::tuple create_row_col_permutation_matrices(bool create_row, bool create_col); bool checkpoint_enabled_ = false; - PetscInt current_gf_ = 0; - PetscInt n_gf_ = 0; Mat G_ = nullptr; std::unique_ptr S_; std::unique_ptr t_boundary_; diff --git a/app/tandem/SEAS.cpp b/app/tandem/SEAS.cpp index ee0b471e..9444478f 100644 --- a/app/tandem/SEAS.cpp +++ b/app/tandem/SEAS.cpp @@ -140,42 +140,18 @@ template <> struct operator_specifics : public qd_operator_specifics {}; - -template struct qd_green_operator_specifics { - using monitor_t = seas::MonitorQD; +template <> struct operator_specifics : public qd_operator_specifics { static auto make(LocalSimplexMesh const& mesh, Config const& cfg, seas::ContextBase& ctx) { - std::string prefix; - if (cfg.gf_checkpoint_prefix) { prefix = *(cfg.gf_checkpoint_prefix); } - else { prefix = ""; } - auto seasop = std::make_shared(std::move(ctx.dg()), std::move(ctx.adapter()), - std::move(ctx.friction()), mesh, cfg.matrix_free, - MGConfig(cfg.mg_coarse_level, cfg.mg_strategy), prefix); + auto seasop = std::make_shared(std::move(ctx.dg()), std::move(ctx.adapter()), + std::move(ctx.friction()), mesh, cfg.gf_checkpoint_prefix, cfg.matrix_free, + MGConfig(cfg.mg_coarse_level, cfg.mg_strategy)); ctx.setup_seasop(*seasop); seasop->warmup(); return seasop; } - - static std::optional cfl_time_step(T const&) { return std::nullopt; } - - template - static auto displacement(double time, TimeSolver const& ts, T& seasop) { - seasop.update_internal_state(time, ts.state(0), true, false, true); - return seasop.displacement(); - } - - template static auto state(double, TimeSolver const& ts, T& seasop) { - return seasop.friction().raw_state(ts.state(0)); - } - - static void print_profile(T const&) {} }; - -template <> -struct operator_specifics - : public qd_green_operator_specifics {}; - template <> struct operator_specifics { using monitor_t = seas::MonitorFD; From 8459461a158ab9f581004871ca65d34c4dcdfda3 Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Fri, 15 Sep 2023 10:29:11 +0200 Subject: [PATCH 2/3] fix 2 bugs --- app/form/SeasQDDiscreteGreenOperator.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/app/form/SeasQDDiscreteGreenOperator.cpp b/app/form/SeasQDDiscreteGreenOperator.cpp index f5140462..4dbd81cc 100644 --- a/app/form/SeasQDDiscreteGreenOperator.cpp +++ b/app/form/SeasQDDiscreteGreenOperator.cpp @@ -303,6 +303,7 @@ void SeasQDDiscreteGreenOperator :: write_facet_labels_IS(LocalSimplexMesh SeasQDDiscreteGreenOperator :: create_row_col_permutation_m // Rperm is the transpose of Cperm just with a different block size auto create_permutation_matrix = [this,&comm, mb_offset,&fault_map_index](PetscInt m, PetscInt M, PetscInt block_size, bool swap) { - Mat perm = NULL, Rperm = NULL; + Mat perm = NULL; CHKERRTHROW(MatCreate(comm,&perm)); CHKERRTHROW(MatSetSizes(perm,m,m,M,M)); CHKERRTHROW(MatSetType(perm,MATAIJ)); - CHKERRTHROW(MatSeqAIJSetPreallocation(Rperm,block_size,NULL)); - CHKERRTHROW(MatMPIAIJSetPreallocation(Rperm,block_size,NULL,block_size,NULL)); + CHKERRTHROW(MatSeqAIJSetPreallocation(perm,block_size,NULL)); + CHKERRTHROW(MatMPIAIJSetPreallocation(perm,block_size,NULL,block_size,NULL)); for (std::size_t bndNo = 0; bndNo < this->adapter().fault_map().local_size(); ++bndNo) { PetscInt from = mb_offset + bndNo; @@ -598,7 +599,7 @@ void SeasQDDiscreteGreenOperator::get_discrete_greens_function(LocalSimplexMesh< // Write out the operator whenever the fully assembled operator was not loaded from file if (n_gfloaded != n_gf) { - write_discrete_greens_operator(mesh, n_gfloaded, n_gf); + write_discrete_greens_operator(mesh, n_gf, n_gf); } } From e2da7970f14869d5f800835db5982a18d43e540f Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Fri, 15 Sep 2023 10:56:41 +0200 Subject: [PATCH 3/3] add gf_checkpoint_every_nmins parameter --- app/form/SeasQDDiscreteGreenOperator.cpp | 4 ++++ app/form/SeasQDDiscreteGreenOperator.h | 6 ++++-- app/tandem/SEAS.cpp | 2 +- app/tandem/SeasConfig.cpp | 7 ++++--- app/tandem/SeasConfig.h | 1 + 5 files changed, 14 insertions(+), 6 deletions(-) diff --git a/app/form/SeasQDDiscreteGreenOperator.cpp b/app/form/SeasQDDiscreteGreenOperator.cpp index 4dbd81cc..1eb28e28 100644 --- a/app/form/SeasQDDiscreteGreenOperator.cpp +++ b/app/form/SeasQDDiscreteGreenOperator.cpp @@ -16,6 +16,7 @@ SeasQDDiscreteGreenOperator::SeasQDDiscreteGreenOperator( std::unique_ptr friction, LocalSimplexMesh const& mesh, std::optional prefix, + double gf_checkpoint_every_nmins, bool matrix_free, MGConfig const& mg_config) : base(std::move(dgop), std::move(adapter), std::move(friction), matrix_free, mg_config) { @@ -23,6 +24,9 @@ SeasQDDiscreteGreenOperator::SeasQDDiscreteGreenOperator( MPI_Comm_rank(base::comm(),&rank); // if prefix is not empty, set filenames and mark checkpoint_enabled_ = true + + checkpoint_every_nmins_ = gf_checkpoint_every_nmins; + if (prefix.has_value()) { std::string sprefix = prefix.value_or(""); if (rank == 0) { std::cout << "Using GF checkpoint path: " << sprefix << std::endl; } diff --git a/app/form/SeasQDDiscreteGreenOperator.h b/app/form/SeasQDDiscreteGreenOperator.h index db8b7759..3ae14149 100644 --- a/app/form/SeasQDDiscreteGreenOperator.h +++ b/app/form/SeasQDDiscreteGreenOperator.h @@ -25,7 +25,9 @@ class SeasQDDiscreteGreenOperator : public SeasQDOperator { SeasQDDiscreteGreenOperator(std::unique_ptr dgop, std::unique_ptr adapter, std::unique_ptr friction, -LocalSimplexMesh const& mesh, std::optional prefix, + LocalSimplexMesh const& mesh, + std::optional prefix, + double gf_checkpoint_every_nmins, bool matrix_free = false, MGConfig const& mg_config = MGConfig()); ~SeasQDDiscreteGreenOperator(); @@ -54,7 +56,7 @@ LocalSimplexMesh const& mesh, std::optional prefix std::string gf_operator_filename_ = "gf_mat.bin"; std::string gf_traction_filename_ = "gf_vec.bin"; std::string gf_facet_filename_ = "gf_facet_labels.bin"; - double checkpoint_every_nmins_ = 30.0; + double checkpoint_every_nmins_; void update_traction(double time, BlockVector const& state); diff --git a/app/tandem/SEAS.cpp b/app/tandem/SEAS.cpp index 9444478f..b19bc416 100644 --- a/app/tandem/SEAS.cpp +++ b/app/tandem/SEAS.cpp @@ -144,7 +144,7 @@ template <> struct operator_specifics : public qd_o static auto make(LocalSimplexMesh const& mesh, Config const& cfg, seas::ContextBase& ctx) { auto seasop = std::make_shared(std::move(ctx.dg()), std::move(ctx.adapter()), - std::move(ctx.friction()), mesh, cfg.gf_checkpoint_prefix, cfg.matrix_free, + std::move(ctx.friction()), mesh, cfg.gf_checkpoint_prefix, cfg.gf_checkpoint_every_nmins, cfg.matrix_free, MGConfig(cfg.mg_coarse_level, cfg.mg_strategy)); ctx.setup_seasop(*seasop); seasop->warmup(); diff --git a/app/tandem/SeasConfig.cpp b/app/tandem/SeasConfig.cpp index b250b838..776853c6 100644 --- a/app/tandem/SeasConfig.cpp +++ b/app/tandem/SeasConfig.cpp @@ -117,11 +117,12 @@ void setConfigSchema(TableSchema& schema, .default_value(false) .help("Assert that boundary is a linear function of time (i.e. boundary(x, t) = f(x) t)."); - schema.add_value("gf_checkpoint_prefix", &Config::gf_checkpoint_prefix) .help("Path where Green's function operator and RHS will be checkpointed"); - - + schema.add_value("gf_checkpoint_every_nmins", &Config::gf_checkpoint_every_nmins) + .default_value(30.0) + .help("time interval, in minutes, at which the Green's function operator data is saved to disk"); + schema.add_value("matrix_free", &Config::matrix_free) .default_value(false) .help("Use matrix-free operators"); diff --git a/app/tandem/SeasConfig.h b/app/tandem/SeasConfig.h index e0856b81..17735333 100644 --- a/app/tandem/SeasConfig.h +++ b/app/tandem/SeasConfig.h @@ -76,6 +76,7 @@ struct Config { unsigned mg_coarse_level; std::optional gf_checkpoint_prefix; + double gf_checkpoint_every_nmins; std::optional> generate_mesh; std::optional fault_output;