From 75e0547c2a27272c34519eafc555ac0ad7b58129 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Fri, 17 Apr 2015 11:59:26 -0400 Subject: [PATCH] step-50 debugging work --- examples/step-50/step-50.cc | 314 +++++++++++++++++++++++++++++++----- 1 file changed, 277 insertions(+), 37 deletions(-) diff --git a/examples/step-50/step-50.cc b/examples/step-50/step-50.cc index b172983285..cd18fec6a0 100644 --- a/examples/step-50/step-50.cc +++ b/examples/step-50/step-50.cc @@ -105,6 +105,8 @@ #include #include #include +#include +#include // The last step is as in all // previous programs: @@ -113,6 +115,152 @@ namespace Step50 using namespace dealii; + + + + // send to proc 0 and compute&print checksum + void check(const ConstraintMatrix & cm) + { + unsigned int myid=Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + std::string rank0data; + { + const IndexSet & is = cm.get_local_lines(); + unsigned int num_lines = cm.n_constraints(); + std::ostringstream oss; + { + boost::archive::text_oarchive oa(oss); + oa << num_lines; + for (unsigned int n=0;n > * v + = cm.get_constraint_entries(idx); + oa << idx; + unsigned int s = v->size(); + oa << s; + for (unsigned int i=0;isize();++i) + oa << (*v)[i].first << (*v)[i].second; + double inhom = cm.get_inhomogeneity(idx); + oa << inhom; + //std::cout << "check() idx: " << idx << " v: " << v->size() << " inhom: " << inhom << std::endl; + } + } + } +// std::cout << "check(): myid=" << myid +// << " data length: " << oss.str().size() +// << " lines: " << num_lines +// << std::endl; + std::string data = oss.str(); + char * ptr = (char*)data.c_str(); + if (myid!=0) + MPI_Send(ptr,data.size()+1,MPI_CHAR, 0,7001,MPI_COMM_WORLD); + else + rank0data = data; + } + + MPI_Barrier(MPI_COMM_WORLD); + if (myid==0) + { + ConstraintMatrix all; + unsigned int numprocs = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); + std::vector receive; + for (unsigned int p=0;p> num_lines; + std::cout << "check(): got " + << " lines: " << num_lines + << " from " << status.MPI_SOURCE + << std::endl; + + for (unsigned int line=0;line > d; + double inhom; + unsigned int ll; + + ia >> idx >> ll; + for (unsigned int i=0;i> idx >> val; + d.push_back(std::make_pair(idx, val)); + } + ia >> inhom; +// std::cout << "check() idx: " << idx << " d: " << d.size() << " inhom: " << inhom << std::endl; + if (!all.is_constrained(idx)) + { + all.add_line(idx); + all.add_entries(idx,d); + all.set_inhomogeneity(idx, inhom); + } + else + { + const std::vector > * v + = all.get_constraint_entries(idx); + if (v->size() != d.size()) + std::cout << "check(): DIFFERENT size " << idx << std::endl; + else + { + for (unsigned int i=0;iLaplaceProblem class template} // This main class is basically the same @@ -455,17 +603,6 @@ namespace Step50 - // Here we output not only the - // degrees of freedom on the finest - // level, but also in the - // multilevel structure - deallog << "Number of degrees of freedom: " - << mg_dof_handler.n_dofs(); - - for (unsigned int l=0; l(buffers2[index].c_str()); MPI_Isend(ptr2, buffers2[index].length(), @@ -858,18 +1020,33 @@ namespace Step50 std::istringstream iss(std::string(ptr,len)); IndexSet reply_ref_edge; - IndexSet reply_ref_edge_bdry; + //IndexSet reply_ref_edge_bdry; reply_ref_edge.block_read(iss); - reply_ref_edge_bdry.block_read(iss); + IndexSet reply_bdry; + reply_bdry.block_read(iss); + //reply_ref_edge_bdry.block_read(iss); unsigned int c1 = mg_constrained_dofs.refinement_edge_indices[l].n_elements(); - unsigned int c2 = mg_constrained_dofs.refinement_edge_boundary_indices[l].n_elements(); + //unsigned int c2 = mg_constrained_dofs.refinement_edge_boundary_indices[l].n_elements(); mg_constrained_dofs.refinement_edge_indices[l].add_indices(reply_ref_edge); - mg_constrained_dofs.refinement_edge_boundary_indices[l].add_indices(reply_ref_edge_bdry); - - std::cout << "new " << mg_constrained_dofs.refinement_edge_indices[l].n_elements()-c1 - << " and " << mg_constrained_dofs.refinement_edge_boundary_indices[l].n_elements()-c2 + //mg_constrained_dofs.refinement_edge_boundary_indices[l].add_indices(reply_ref_edge_bdry); + + unsigned int c3 = mg_constrained_dofs.boundary_indices[l].size(); + for (unsigned int n=0;nc3) + std::cout << "HEY: boundary indices transfer did something!" << std::endl; + } } @@ -903,9 +1080,13 @@ namespace Step50 for (unsigned int level=0; level to_refine; @@ -1265,20 +1446,60 @@ namespace Step50 triangulation.prepare_coarsening_and_refinement (); - MPI_Barrier(MPI_COMM_WORLD); - typename DoFHandler::active_cell_iterator - cell = mg_dof_handler.begin_active(), - endc = mg_dof_handler.end(); - for (; cell!=endc; ++cell) - if (cell->is_locally_owned() && cell->refine_flag_set()) - { - std::cout << cell->id() << std::endl; - } - MPI_Barrier(MPI_COMM_WORLD); - } + if (Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)>1) + { + // write refine.n + std::ofstream f; + int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + if (myid==0) + f.open(("refine." + Utilities::int_to_string(n)).c_str()); + { + // local contribution + std::ostringstream oss; + typename DoFHandler::active_cell_iterator + cell = mg_dof_handler.begin_active(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) + { + if (cell->is_locally_owned() && cell->refine_flag_set()) + oss << cell->id() << std::endl; + + cell->clear_coarsen_flag(); + + } + + + if (myid!=0) + { + std::string s= oss.str(); + + MPI_Send((void*)s.c_str(), s.size(), MPI_CHAR, 0, 0, MPI_COMM_WORLD); + } + else + f << oss.str(); + } + + if (myid==0) + for (unsigned int i=0;i receive; + receive.resize(len+1); + + char *ptr = &receive[0]; + MPI_Recv(ptr, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG, + MPI_COMM_WORLD, &status); + f << ptr; + } + } + } + triangulation.execute_coarsening_and_refinement (); } @@ -1340,6 +1561,9 @@ namespace Step50 ".visit"); std::ofstream visit_master (visit_master_filename.c_str()); data_out.write_visit_record (visit_master, filenames); + + std::cout << "wrote " << pvtu_master_filename << std::endl; + } } @@ -1357,15 +1581,16 @@ namespace Step50 template void LaplaceProblem::run () { - for (unsigned int cycle=0; cycle<5; ++cycle) + for (unsigned int cycle=0; cycle<7; ++cycle) { deallog << "Cycle " << cycle << ':' << std::endl; + std::cout << "Cycle " << cycle << ':' << std::endl; if (cycle == 0) { GridGenerator::hyper_cube (triangulation); - triangulation.refine_global (3); + triangulation.refine_global (4); } else refine_grid (); @@ -1425,9 +1650,24 @@ namespace Step50 - output_results (cycle); solve (); + output_results (cycle); + + MPI_Barrier(MPI_COMM_WORLD); + TrilinosWrappers::MPI::Vector temp = solution; + system_matrix.residual(temp,solution,system_rhs); + constraints.set_zero(temp); + deallog << "residual " << temp.l2_norm() << std::endl; + if (temp.l2_norm()>1e-5) + { + for (unsigned int i=0;i0.06) + std::cout << "RESIDUAL at " << i + << " is " << temp[i] << std::endl; + + } + Assert(temp.l2_norm()<1e-5, ExcMessage("Error: residual is large for some reason!")); } } @@ -1447,7 +1687,7 @@ int main (int argc, char *argv[]) using namespace dealii; using namespace Step50; - LaplaceProblem<2> laplace_problem(1); + LaplaceProblem<2> laplace_problem(1/*degree*/); laplace_problem.run (); } catch (std::exception &exc) -- 2.39.5