Best way to transfer old solutions in time stepping over mesh refinement?

[Maxi Miller]

I am trying to implement example 26 using Newton's method and Sacado, in order to get a better feeling for the method. As mentioned before, I had problems with retaining the old solution when remeshing. Thus I rewrote my remeshing sequence as

template <int dim>
    void HeatEquation<dim>::refine_mesh ()
    {
        Vector<float> estimated_error_per_cell (triangulation.n_active_cells());

        KellyErrorEstimator<dim>::estimate (dof_handler,
                                            QGauss<dim-1>(fe.degree+1),
                                            typename FunctionMap<dim>::type(),
                                            present_solution,
                                            estimated_error_per_cell);

        GridRefinement::refine_and_coarsen_fixed_number (triangulation,
                                                         estimated_error_per_cell,
                                                         0.6, 0.4);

        if (triangulation.n_levels() > max_grid_level)
            for (auto cell = triangulation.begin_active(max_grid_level); cell != triangulation.end(); ++cell)
                cell->clear_refine_flag ();
        for (auto cell = triangulation.begin_active(min_grid_level); cell != triangulation.end_active(min_grid_level); ++cell)
            cell->clear_coarsen_flag ();

        triangulation.prepare_coarsening_and_refinement ();

        SolutionTransfer<dim> solution_transfer(dof_handler);
        solution_transfer.prepare_for_coarsening_and_refinement(present_solution);

        triangulation.execute_coarsening_and_refinement();

        dof_handler.distribute_dofs(fe);

        Vector<double> tmp(dof_handler.n_dofs());
        solution_transfer.interpolate(present_solution, tmp);
        present_solution = tmp;
        solution_transfer.interpolate(old_solution, tmp);
        old_solution = tmp;

        set_boundary_values ();


        constraint_matrix.clear();

        DoFTools::make_hanging_node_constraints(dof_handler,
                                                constraint_matrix);
        constraint_matrix.close();

        constraint_matrix.distribute (present_solution);
        constraint_matrix.distribute(old_solution);

        setup_system (false);
    }

and then I use the same grid both for the new and the old solution for the fe_values.grad_values(). Is that approach correct? I assume it still introduces some errors, but I am not sure if they are responsible for my (wrong) results:


It appears as if the heat is never spreading out, but also the maximum temperature is way below the example results, even though I used exactly the same boundary/Input values.

[Jean-Paul Pelteret]

Dear Maxi,

I'm going to put my foot down and say that, particularly as of late, you're not giving your problems sufficient attention before posting them to the forum. We have a set of guidelines

https://groups.google.com/forum/#!topic/dealii/GRZMUTLIm2I

(which, to the best of my recollection, I've mentioned to you before) that includes a minimum set of expectations that you should consider before asking a question. Given the time that elapsed between this post and your last (as well as those preceding it), I judge that you are not giving your issues due consideration. Everyone runs into these sorts of issues when developing a code, and it takes time and patience, and sometimes a little playing around, to think about what the problem might be. 

Please respect the time of the group moderators, developers and other community members. We all want to help you succeed in solving your problem, but the google group is not meant to be a one-stop shop for trouble-shooting and debugging your code.

Good luck in resolving your issue.

Jean-Paul

[Wolfgang Bangerth]

On 08/29/2017 02:44 AM, 'Maxi Miller' via deal.II User Group wrote:
> and then I use the same grid both for the new and the old solution for the
> fe_values.grad_values(). Is that approach correct? I assume it still
> introduces some errors,

For the record, yes, that's what we typically do -- transfer the old solution
to the new mesh and then continue working only on the new mesh. There are
alternative approaches whereby you keep the old solution on the old mesh, but
it's cumbersome to work with two meshes.

As for the reason why your program doesn't seem to work -- Jean-Paul already
replied to that.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/