Hello Prof. Bangerth,
I have attached a test case.
Best,
Yiyang
On Friday, October 27, 2017 at 5:12:49 PM UTC-5, Wolfgang Bangerth wrote:
>
> On 10/27/2017 02:13 PM, Yiyang Zhang wrote:
> >
> > Yes I think I am setting them in a way that is exactly same throughout
> > all processes.
> > Since I can check the n_active_cells() before the refinement, and also
> > the refine_flags and coarsen_flags for each process. They are the same
> > for each process. But after the refinement, the number of active_cells
> > suddenly becomes different.
> >
> > I attached the refine_mesh() code.
>
> Can you fabricate a complete testcase that we can run that shows this
> phenomenon?
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bang...@colostate.edu
> <javascript:>
> www: http://www.math.colostate.edu/~bangerth/
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/usr/local/Cellar/open-mpi/2.1.1/bin/mpirun -np 2 ./tria
Process No. 0 ; distributed size = 10240
Process No. 1 ; distributed size = 10240
Process No. 0. Refine Count = 8194
Process No. 0. Coarsen Count = 2048
Process No. 1. Refine Count = 8194
Process No. 1. Coarsen Count = 2048
Process No. 0. Pre_refine, n_active_cells = 20480
Process No. 1. Pre_refine, n_active_cells = 20480
Process No. 0. Post_refine, n_active_cells = 43772
Process No. 1. Post_refine, n_active_cells = 43772
Process No. 1 ; distributed size = 21886
Process No. 0 ; distributed size = 21886
Process No. 1. Refine Count = 6144
Process No. 1. Coarsen Count = 1442
Process No. 0. Refine Count = 6144
Process No. 0. Coarsen Count = 1442
Process No. 1. Pre_refine, n_active_cells = 43772
Process No. 0. Pre_refine, n_active_cells = 43772
Process No. 0. Post_refine, n_active_cells = 61367
Process No. 1. Post_refine, n_active_cells = 61364
Active cells not synced.
Active cells not synced.
Process No. 0 ; distributed size = 30684
Process No. 1 ; distributed size = 30682
libc++abi.dylib: terminating with uncaught exception of type
dealii::StandardExceptions::ExcDimensionMismatch:
--------------------------------------------------------
An error occurred in line <792> of file
</tmp/dealii-20170729-98705-mhpw56/dealii-8.5.0/source/lac/trilinos_vector.cc>
in function
void dealii::TrilinosWrappers::Vector::reinit(const
dealii::TrilinosWrappers::VectorBase &, const bool, const bool)
The violated condition was:
size() == v.size()
Additional information:
Dimension 61366 not equal to 52139.
--------------------------------------------------------
libc++abi.dylib: terminating with uncaught exception of type
dealii::StandardExceptions::ExcDimensionMismatch:
--------------------------------------------------------
An error occurred in line <792> of file
</tmp/dealii-20170729-98705-mhpw56/dealii-8.5.0/source/lac/trilinos_vector.cc>
in function
void dealii::TrilinosWrappers::Vector::reinit(const
dealii::TrilinosWrappers::VectorBase &, const bool, const bool)
The violated condition was:
size() == v.size()
Additional information:
Dimension 61366 not equal to 52139.
#ifndef SHARED_TRIA_TEST
#define SHARED_TRIA_TEST
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/table_handler.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparsity_tools.h> //need SparsityTools::distribute_sparsity_pattern()
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/trilinos_block_vector.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <algorithm>
#include <math.h>
//elliptic integral of first kind
#include <boost/math/special_functions/ellint_1.hpp>
//elliptic integral of third kind
#include <boost/math/special_functions/ellint_3.hpp>
namespace TEST{
using namespace dealii;
const double eta = 1.0;
const double lambda = 1.0;
const double PI =3.1415926535897932384626;
template<int dim>
class UnrelaxedLoop : public Function<dim> {
private:
double radius;
//RO is the radius of the loop
const double loop_magnitude(const Point<dim>& p, const double R0 = 1.0) const {
const double rho = sqrt(p(0)*p(0) + p(1)*p(1));
double dist2;
if (dim == 2) {
dist2 = pow(rho - R0, 2);
}
else {
dist2 = pow(rho - R0, 2) + p(2)*p(2);
}
return 1.0 - exp(-dist2);/*divided by 1 length^2*/
}
const double loop_solid_angle(const Point<dim>& p, const double R0 = 1.0 /*radius of the loop*/) const {
double val;
const double dist = std::max(p.distance(Point<dim>()),1e-10);
const double lam = dist / R0;
double theta;
if (dim == 3) {
theta = acos(p(2) / dist);
}
else {
theta = PI / 2.0;
}
if (abs(theta - PI / 2.0) < 1e-6) {
if (lam < 1) return 2.0*PI;
else return 0.0;
}
else {
const double coef = 1.0 + lam*lam + 2.0*lam*sin(theta);
const double n = 2.0*sin(theta) / (1 + sin(theta));
const double m = n*pow(1.0 + lam, 2) / coef;
const double k = sqrt(4.0*lam*sin(theta) / coef);
const double ellipticPi_m = boost::math::ellint_3(k, m);
const double ellipticPi_n = boost::math::ellint_3(k, n);
const double ellipticK = boost::math::ellint_1(k);
return 2.0*PI - 2.0*cos(theta) / sqrt(coef)*(
(1.0 + lam*lam - 2.0*lam*sin(theta)) / ((1.0 + lam)*(1.0 + sin(theta)))*ellipticPi_m
- (1.0 - lam) / (1.0 + sin(theta))*ellipticPi_n
+ 2.0*lam / (1.0 + lam)*ellipticK);
}
}
public:
UnrelaxedLoop(const double r = 1.0) :Function<dim>(2) {
radius = r;
}
virtual double value(const Point<dim>& p, const unsigned int component = 0) const {
Assert(component < this->n_components, ExcIndexRange(component, 0, this->n_components));
double val;
const double mag = loop_magnitude(p, radius);
const double arg = loop_solid_angle(p, radius);
if (component == 0) {
val = eta*mag*cos(arg / 2.0);
}
else {
val = eta*mag*sin(arg / 2.0);
}
return val;
}
virtual void vector_value(const Point<dim>& p, Vector<double>& values) const {
for (unsigned int c = 0; c < this->n_components; ++c) {
values(c) = this->value(p, c);
}
return;
}
};
template<int dim>
class TestFunc : public Function<dim> {
private:
double radius;
//RO is the radius of the loop
const double loop_magnitude(const Point<dim>& p, const double R0 = 1.0) const {
const double rho = sqrt(p(0)*p(0) + p(1)*p(1));
double dist2;
if (dim == 2) {
dist2 = pow(rho - R0, 2);
}
else {
dist2 = pow(rho - R0, 2) + p(2)*p(2);
}
return 1.0 - exp(-dist2);/*divided by 1 length^2*/
}
public:
TestFunc(const double r = 1.0) :Function<dim>(2) {
radius = r;
}
virtual double value(const Point<dim>& p, const unsigned int component = 0) const {
Assert(component < this->n_components, ExcIndexRange(component, 0, this->n_components));
double val;
const double dist = p.distance(Point<dim>());
const double mag = loop_magnitude(p, radius);
if(component == 0){
val = mag * p(0)/dist;
}
if(component == 1){
val = mag * p(1)/dist;
}
return val;
}
virtual void vector_value(const Point<dim>& p, Vector<double>& values) const {
for (unsigned int c = 0; c < this->n_components; ++c) {
values(c) = this->value(p, c);
}
return;
}
};
template<int dim>
class SharedTria{
public:
SharedTria();
~SharedTria();
void run();
private:
void make_mesh(const double size, const int refine);
void make_ball_mesh_2d(const double radius, const int refine);
void setup_system();
void init_sol(const double radius);
void refine_mesh();
parallel::shared::Triangulation<dim> tria;
DoFHandler<dim> dof_handler;
FESystem<dim> fe;
SphericalManifold<dim> spherical_mfd;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
ConstraintMatrix constraints;
TrilinosWrappers::MPI::Vector locally_relevant_sln;
MPI_Comm mpi_communicator;
};
template<int dim>
SharedTria<dim>::SharedTria()
:
dof_handler(tria),
mpi_communicator(MPI_COMM_WORLD),
fe(FE_Q<dim>(1), 2),
spherical_mfd(),
tria(MPI_COMM_WORLD, typename Triangulation<dim>::MeshSmoothing
(Triangulation<dim>::smoothing_on_refinement | Triangulation<dim>::smoothing_on_coarsening) )
{}
template<int dim>
SharedTria<dim>::~SharedTria(){
dof_handler.clear();
}
template<int dim>
void SharedTria<dim>::init_sol(const double radius) {
TrilinosWrappers::MPI::Vector distributed_sln(this->locally_owned_dofs, this->mpi_communicator);
VectorTools::interpolate(
this->dof_handler,
UnrelaxedLoop<dim>(radius),
distributed_sln);
//ConstraintMatrix::distribute() can only be used on vectors without ghose elements.
this->constraints.distribute(distributed_sln);
locally_relevant_sln = distributed_sln;
return;
}
template<int dim>
void SharedTria<dim>::make_mesh(const double size, const int refine){
GridGenerator::hyper_cube(tria, -size, size);
tria.refine_global(refine);
return;
}
template<int dim>
void SharedTria<dim>::make_ball_mesh_2d(const double radius, const int refine){
//adapt the strategy from Step-6
const double isCenter = 1e-3; //some small value
GridGenerator::hyper_ball(tria, Point<dim>(), radius);
//set manifold id to get spherical manifold description
types::boundary_id MFD_ID_SPH = 1;
tria.set_all_manifold_ids_on_boundary(MFD_ID_SPH);
tria.set_manifold(MFD_ID_SPH, spherical_mfd);
const Point<dim> mesh_center;
const double core_radius = 1.0 / 100.0 * radius; //This is essentially the central resolution.
const double inner_radius = 1.0 / 50.0 * radius;
// Step 1: Shrink the inner cell
//We cannot subdivide the inner four cells (or their faces) onto concentric rings.
for (typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active();
cell != tria.end(); ++cell) {
if (mesh_center.distance(cell->center()) < isCenter)
//cell center almost on the mesh_center
{
for (unsigned int v = 0;
v < GeometryInfo<dim>::vertices_per_cell;
++v)
cell->vertex(v) *= core_radius / mesh_center.distance(cell->vertex(v));
//the vertices of this cell now have a radius of core_radius.
}
}
// Step 2: Refine all cells except the central one
for (typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active();
cell != tria.end(); ++cell) {
if (mesh_center.distance(cell->center()) >= core_radius/*isCenter*/)
cell->set_refine_flag();
}
tria.execute_coarsening_and_refinement();
// Step 3: Resize the inner children of the outer cells
for (typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active();
cell != tria.end(); ++cell) {
for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
{
const double dist = mesh_center.distance(cell->vertex(v));
if (dist > core_radius*1.0001 && dist < 0.9999*radius)
cell->vertex(v) *= inner_radius / dist;
}
}
// Step 4: Apply curved manifold description
for (typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active();
cell != tria.end(); ++cell) {
bool is_in_inner_circle = false;
for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
if (mesh_center.distance(cell->vertex(v)) < inner_radius /* modified when step-4 moved up. default: inner_radius*/)
{
is_in_inner_circle = true;
break;
}
if (is_in_inner_circle == false)
cell->set_all_manifold_ids(MFD_ID_SPH);
}
tria.refine_global(refine);
return;
}
template<int dim>
void SharedTria<dim>::setup_system(){
this->dof_handler.distribute_dofs(this->fe);
this->locally_owned_dofs = this->dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(this->dof_handler, this->locally_relevant_dofs);
this->constraints.clear();
this->constraints.reinit(this->locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(this->dof_handler, this->constraints);
this->constraints.close();
locally_relevant_sln.reinit(this->locally_owned_dofs, this->locally_relevant_dofs, this->mpi_communicator);
return;
}
template<int dim>
void SharedTria<dim>::refine_mesh(){
const auto n_active_cells = this->tria.n_active_cells();
const auto max_nac = Utilities::MPI::max(n_active_cells, this->mpi_communicator);
const auto min_nac = Utilities::MPI::min(n_active_cells, this->mpi_communicator);
if(max_nac != min_nac){
std::cerr << "Active cells not synced." << std::endl;
}
Vector<float> estimated_error_per_cell(this->tria.n_active_cells());
KellyErrorEstimator<dim>::estimate(this->dof_handler,
QGauss<dim - 1>(4), //What are these parameters?
typename FunctionMap<dim>::type(),
locally_relevant_sln,
estimated_error_per_cell,
ComponentMask(),
0,
MultithreadInfo::n_threads(),
Utilities::MPI::this_mpi_process(this->mpi_communicator));
IndexSet local_active_cells(this->tria.n_active_cells());
local_active_cells.clear();
for(auto cell = this->tria.begin_active(); cell != this->tria.end(); ++cell){
if(cell->is_locally_owned()){
local_active_cells.add_index(cell->active_cell_index());
}
}
TrilinosWrappers::MPI::Vector distributed_error_per_cell(local_active_cells, this->mpi_communicator);
for(auto cell = this->tria.begin_active(); cell != this->tria.end(); ++cell){
if(cell->is_locally_owned()){
distributed_error_per_cell(cell->active_cell_index()) = estimated_error_per_cell(cell->active_cell_index());
}
}
std::cout<<"Process No. "<<Utilities::MPI::this_mpi_process(this->mpi_communicator) <<" ; distributed size = " << distributed_error_per_cell.local_size()<<std::endl;
distributed_error_per_cell.compress(VectorOperation::insert);
estimated_error_per_cell = distributed_error_per_cell;
GridRefinement::refine_and_coarsen_fixed_number(this->tria,
estimated_error_per_cell,
0.4,
0.1,
60000
);
int refine_count = 0;
int coarsen_count = 0;
for(auto cell = this->tria.begin_active(); cell != this->tria.end(); ++cell){
if(cell->refine_flag_set()) refine_count++;
if(cell->coarsen_flag_set()) coarsen_count++;
}
std::cout<<"Process No. "<<Utilities::MPI::this_mpi_process(this->mpi_communicator) <<". Refine Count = "<<refine_count<<std::endl;
std::cout<<"Process No. "<<Utilities::MPI::this_mpi_process(this->mpi_communicator) <<". Coarsen Count = "<<coarsen_count<<std::endl;
this->tria.prepare_coarsening_and_refinement();
std::cout<<"Process No. "<<Utilities::MPI::this_mpi_process(this->mpi_communicator) <<". Pre_refine, n_active_cells = "<<this->tria.n_active_cells()<<std::endl;
this->tria.execute_coarsening_and_refinement();
std::cout<<"Process No. "<<Utilities::MPI::this_mpi_process(this->mpi_communicator) <<". Post_refine, n_active_cells = "<<this->tria.n_active_cells()<<std::endl;
setup_system();
return;
}
template<int dim>
void SharedTria<dim>::run(){
//with make_mesh() function, the program also collapse, although the number of active cells are the same for both MPI processes.
//make_mesh(10.0, 5);
//with make_ball_mesh_2d() function, we can have a case, where the refine_flag and coarsen_flag is same for both MPI processes
//but the number of active cells are different after the mesh refinement.
//I am using the complicated function UnrelaxedLoop().
//-- It does depend on the function one uses, also depends on how many times we refine the mesh with refine_global().
//However, if one switches to the function TestFunc() -- a simpler function, one will find cases where the refine_flag and coarsen_flag
//are different for different MPI processes, and thus different number of active cells after mesh refinement. I think this is also not
//right, since the behaviour of the mesh should be exactly the same for different processes.
make_ball_mesh_2d(10.0, 5);
setup_system();
init_sol(5.0);
for(int i=0; i < 10;++i){
refine_mesh();
init_sol(5.0);
}
return;
}
}
#endif
#include "shared_tria_test.h"
int main(int argc, char *argv[]){
using namespace dealii;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, /*max_num_threads*/1);
TEST::SharedTria<2> shared;
shared.run();
return 0;
}
