From 7b2f70cc93d0f8f41455d4b93731ddaa403936bf Mon Sep 17 00:00:00 2001 From: Jiaqi Zhang Date: Mon, 14 Sep 2020 11:30:02 -0400 Subject: [PATCH] fix the error estimator --- examples/step-50/step-50.cc | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/examples/step-50/step-50.cc b/examples/step-50/step-50.cc index 7d3b4cc544..576c9d8a73 100644 --- a/examples/step-50/step-50.cc +++ b/examples/step-50/step-50.cc @@ -439,7 +439,7 @@ private: MatrixFreeActiveMatrix mf_system_matrix; VectorType solution; VectorType right_hand_side; - Vector estimated_error_per_cell; + Vector estimated_error_square_per_cell; MGLevelObject mg_matrix; MGLevelObject mg_interface_in; @@ -1279,11 +1279,11 @@ void LaplaceProblem::estimate() const Coefficient coefficient; - estimated_error_per_cell.reinit(triangulation.n_active_cells()); + estimated_error_square_per_cell.reinit(triangulation.n_active_cells()); using Iterator = typename DoFHandler::active_cell_iterator; - // Assembler for cell residual $h \| f + \epsilon \triangle u \|_K$ + // Assembler for cell residual $h^2 \| f + \epsilon \triangle u \|_K^2$ auto cell_worker = [&](const Iterator & cell, ScratchData &scratch_data, CopyData & copy_data) { @@ -1307,7 +1307,8 @@ void LaplaceProblem::estimate() residual_norm_square += residual * residual * fe_values.JxW(k); } - copy_data.value = cell->diameter() * std::sqrt(residual_norm_square); + copy_data.value = + cell->diameter() * cell->diameter() * residual_norm_square; }; // Assembler for face term $\sum_F h_F \| \jump{\epsilon \nabla u \cdot n} @@ -1356,17 +1357,17 @@ void LaplaceProblem::estimate() } const double h = cell->face(f)->measure(); - copy_data_face.values[0] = 0.5 * h * std::sqrt(jump_norm_square); + copy_data_face.values[0] = 0.5 * h * jump_norm_square; copy_data_face.values[1] = copy_data_face.values[0]; }; auto copier = [&](const CopyData ©_data) { if (copy_data.cell_index != numbers::invalid_unsigned_int) - estimated_error_per_cell[copy_data.cell_index] += copy_data.value; + estimated_error_square_per_cell[copy_data.cell_index] += copy_data.value; for (auto &cdf : copy_data.face_data) for (unsigned int j = 0; j < 2; ++j) - estimated_error_per_cell[cdf.cell_indices[j]] += cdf.values[j]; + estimated_error_square_per_cell[cdf.cell_indices[j]] += cdf.values[j]; }; const unsigned int n_gauss_points = degree + 1; @@ -1410,7 +1411,7 @@ void LaplaceProblem::refine_grid() const double refinement_fraction = 1. / (std::pow(2.0, dim) - 1.); parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number( - triangulation, estimated_error_per_cell, refinement_fraction, 0.0); + triangulation, estimated_error_square_per_cell, refinement_fraction, 0.0); triangulation.execute_coarsening_and_refinement(); } @@ -1445,9 +1446,9 @@ void LaplaceProblem::output_results(const unsigned int cycle) level(cell->active_cell_index()) = cell->level(); data_out.add_data_vector(level, "level"); - if (estimated_error_per_cell.size() > 0) - data_out.add_data_vector(estimated_error_per_cell, - "estimated_error_per_cell"); + if (estimated_error_square_per_cell.size() > 0) + data_out.add_data_vector(estimated_error_square_per_cell, + "estimated_error_square_per_cell"); data_out.build_patches(); -- 2.39.5