AMR , how to pass solution vector to refined mesh

[Jaekwang Kim]

Hi, all. 

I wonder if it is possible to pass solution vector to refined-meshed....

I am solving non-linear problem with iterative method. 

I'd like to use coarse mesh solution as a initial solution to find accurate solution fast. 

However, my present code sets up system again on every refinement cycle. 

So, I lose all solution that I had in coarse mesh. 

I think it is really helpful for me if I knew a way to pass the previous solution to as a initial solution for next iterative method. 

I think it might be complicated... first vector size of solution might different when mesh is refined....

is there anyone who has idea on this ?

Thanks, 

Jaekwang Kim 

template <int dim>

void nonlinear<dim>::run ()

{

    unsigned int cycle_control=4;

    

    for (unsigned int cycle=0; cycle<cycle_control; ++cycle)

    {

        std::cout << "Cycle " << cycle << ':' << std::endl;

        

        if (cycle == 0)

        {

            

            GridGenerator::hyper_cube (triangulation, 0.0, 1.0);

            triangulation.refine_global (1);

            

            for (typename Triangulation<dim>::active_cell_iterator

                 cell=triangulation.begin_active();

                 cell!=triangulation.end(); ++cell)

            {

                for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)

                {

                    if (cell->face(f)->at_boundary() == true)

                    {

                        if ( std::abs(cell->face(f)->center()[0]-0.5)< 1e-12 )

                        {

                            cell->face(f)->set_boundary_id(10);

                        }

                        

                        // Boundary faces

                        static const double tol = 1e-12;

                        if ( std::abs(cell->face(f)->center()[0]-0.0)< tol )

                        {

                            cell->face(f)->set_boundary_id(1);

                        }

                        if ( std::abs(cell->face(f)->center()[0]-1.0)< tol )

                        {

                            cell->face(f)->set_boundary_id(1);

                        }

                        if ( std::abs(cell->face(f)->center()[1]-0.0)< tol )

                        {

                            cell->face(f)->set_boundary_id(1);

                        }

                        if ( std::abs(cell->face(f)->center()[1]-1.0)< tol )

                        {

                            cell->face(f)->set_boundary_id(1);

                        }

                    }

                }

            }

            

            

            std::ofstream out ("grid-1.ucd");

            GridOut grid_out;

            grid_out.write_ucd (triangulation, out);

            

            

        }

        

        else

            

            refine_grid ();

        

        std::cout << "   Number of active cells:       "<< triangulation.n_active_cells() << std::endl;

        

        setup_system ();

        int iteration=0;

     

        assemble_system ();

        solve ();   

        

        Vector<double> difference(dof_handler.n_dofs());

        Vector<double> temp_vector(dof_handler.n_dofs());

        Vector<double> temp_difference(dof_handler.n_dofs());

        

        difference=solution; // std::cout << " ||u_k-u_{k-1}||" << difference.l2_norm() << std::endl;

        previous_solution = solution ;

                

        double omega = 0.5; //relxation control number

        

        int success_step=0 ;

        

        do{

            

            assemble_system();

            solve();

            

            

            difference = solution;

            difference -= previous_solution; //This is temporary step

            

          

            temp_vector=solution;

            temp_vector.add(omega,difference);

            

            temp_difference = temp_vector;

            temp_difference -= solution; //Calculate different again

            

            if (temp_difference.l2_norm()< difference.l2_norm())  //temp_difference = temp_vector - solution , difference=solution-previous_solution

            {

                

                success_step+=1;

                

                difference=temp_vector;

                difference-=solution ;

                solution=temp_vector;

            }

            

            

            previous_solution=solution ;

            

            iteration+=1;

            

            std::ofstream outerror2("iter_error.dat",std::ios::app);

            outerror2 << degree << " " << triangulation.n_active_cells() << " " << iteration << " " << difference.l2_norm() << " " << error << std::endl;

            

        }while (difference.l2_norm() > 0.00001);

        

        std::cout << "   ||u_k-u_{k-1}|| = " << difference.l2_norm() <<  "   Iteration Number: " << iteration << std::endl;

        

        error_evaluation();

        

        std::ofstream outerror("error.dat",std::ios::app);

        outerror << degree << " " << triangulation.n_active_cells() << " " << iteration << " " << error << std::endl;

        

        output_results ();

        

    }

    

}

[Wolfgang Bangerth]

On 03/15/2017 02:37 PM, Jaekwang Kim wrote:
>
> I wonder if it is possible to pass solution vector to refined-meshed....
>
> I am solving non-linear problem with iterative method.
>
> I'd like to use coarse mesh solution as a initial solution to find accurate
> solution fast.
>
> However, my present code sets up system again on every refinement cycle.
> So, I lose all solution that I had in coarse mesh.
>
> I think it is really helpful for me if I knew a way to pass the previous
> solution to as a initial solution for next iterative method.
>
> I think it might be complicated... first vector size of solution might
> different when mesh is refined....
>
> is there anyone who has idea on this ?

Take a look at step-15 -- it does exactly this, using the SolutionTransfer class.

Best
W.

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

[Jaekwang Kim]

Thank you, Dr. Bangerth!

Step-15 tutorial was really helpful for me and I tried to used solution transfer class. 

I refined grid as follow. 

Because I am using Picard iteration, I use it as a previous_solution vector (in first refined mesh) when it is interpolated into refined meshes. 

template <int dim>

void nonlinear<dim>::refine_grid ()

{

    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());

    KellyErrorEstimator<dim>::estimate (dof_handler,

                                        QGauss<dim-1>(3),

                                        typename FunctionMap<dim>::type(),

                                        solution,

                                        estimated_error_per_cell);

    GridRefinement::refine_and_coarsen_fixed_number (triangulation,

                                                     estimated_error_per_cell,

                                                     0.3, 0.0);

    

    /* You need to tranfer solution (n-1 step) domain into previous solution (n step)*/

   

    triangulation.prepare_coarsening_and_refinement ();

    SolutionTransfer<dim> solution_transfer(dof_handler);

    solution_transfer.prepare_for_coarsening_and_refinement(solution);

    triangulation.execute_coarsening_and_refinement(); // I think n_dof() has increased here. You need to check this if necessary

    dof_handler.distribute_dofs(fe);

    

   // Vector<double> tmp(dof_handler.n_dofs()); //tmp - n step dof

    

    previous_solution.reinit (dof_handler.n_dofs());

    

    solution_transfer.interpolate(solution, previous_solution);  // interpolate (input, output)

    

    solution.reinit (dof_handler.n_dofs());

    

}

Also My setup_system function is is as follow 

template <int dim>

void nonlinear<dim>::setup_system (const unsigned int refinement_cycle)

{

    if (refinement_cycle==0)

    { dof_handler.distribute_dofs (fe);

    

    std::cout << "   Number of degrees of freedom: "

    << dof_handler.n_dofs()

    << std::endl;

    

    solution.reinit (dof_handler.n_dofs());

    previous_solution.reinit (dof_handler.n_dofs());

    system_rhs.reinit (dof_handler.n_dofs());

    

    

    for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)//////////

    {                                                //////

        previous_solution(i)=1;

    }

    

    

    DynamicSparsityPattern dsp(dof_handler.n_dofs());

    DoFTools::make_sparsity_pattern (dof_handler, dsp);

    sparsity_pattern.copy_from(dsp);

    

    system_matrix.reinit (sparsity_pattern);

    

    

    

    constraints.clear ();

    DoFTools::make_hanging_node_constraints (dof_handler,

                                             constraints);

    

    VectorTools::interpolate_boundary_values (dof_handler,1,

                                              BoundaryValues<dim>(),

                                              constraints);

    

    constraints.close ();

        

    }

    else

    {

        

        //You don't need to make previous solution. Or solution vector , (it will be done on refinement step)

        //What you should do?  set boundary condition again.

        DynamicSparsityPattern dsp(dof_handler.n_dofs());

        DoFTools::make_sparsity_pattern (dof_handler, dsp);

        sparsity_pattern.copy_from(dsp);

        

        system_matrix.reinit (sparsity_pattern);

        system_rhs.reinit (dof_handler.n_dofs());

        

        

        constraints.clear ();

        DoFTools::make_hanging_node_constraints (dof_handler,

                                                 constraints);

        VectorTools::interpolate_boundary_values (dof_handler,1,

                                                  BoundaryValues<dim>(),

                                                  constraints);

        constraints.close ();

        

        std::cout << "Set up system finished" << std::endl;

    }

    

}

I think this two function do everything I needed... but when I run my code, I run into error message when I assemble system in second time. 

what might be the problem.....?

Always thank you!!

Jaekwang Kim 

An error occurred in line <1668> of file </Users/kimjaekwang/dealii-8.4.1/include/deal.II/lac/constraint_matrix.templates.h> in function

    void dealii::internals::dealiiSparseMatrix::add_value(const LocalType, const size_type, const size_type, SparseMatrixIterator &) [SparseMatrixIterator = dealii::SparseMatrixIterators::Iterator<double, false>, LocalType = double]

The violated condition was: 

    matrix_values->column() == column

The name and call sequence of the exception was:

    typename SparseMatrix<typename SparseMatrixIterator::MatrixType::value_type>::ExcInvalidIndex(row, column)

Additional Information: 

You are trying to access the matrix entry with index <0,30>, but this entry does not exist in the sparsity pattern of this matrix.

The most common cause for this problem is that you used a method to build the sparsity pattern that did not (completely) take into account all of the entries you will later try to write into. An example would be building a sparsity pattern that does not include the entries you will write into due to constraints on degrees of freedom such as hanging nodes or periodic boundary conditions. In such cases, building the sparsity pattern will succeed, but you will get errors such as the current one at one point or other when trying to write into the entries of the matrix. 

[Wolfgang Bangerth]

On 03/16/2017 08:38 PM, Jaekwang Kim wrote:
>
> I think this two function do everything I needed... but when I run my code, I
> run into error message when I assemble system in second time.
> what might be the problem.....?

As the error message says, you are trying to write into a matrix entry that
does not exist. Did you re-initialize the sparsity pattern and matrix after
you refined the mesh?

[Jaekwang Kim]

After the first cycle, I am doing

template <int dim>

void nonlinear<dim>::refine_grid ()

{

Vector<float> estimated_error_per_cell (triangulation.n_active_cells());

    KellyErrorEstimator<dim>::estimate (dof_handler,

                                        QGauss<dim-1>(3),

                                        typename FunctionMap<dim>::type(),

                                        solution,

                                        estimated_error_per_cell);

    GridRefinement::refine_and_coarsen_fixed_number (triangulation,

                                                     estimated_error_per_cell,

                                                     0.3, 0.0);

    

   

    triangulation.prepare_coarsening_and_refinement ();

    SolutionTransfer<dim> solution_transfer(dof_handler);

    solution_transfer.prepare_for_coarsening_and_refinement(solution);

    triangulation.execute_coarsening_and_refinement(); // I think n_dof() has increased here. You need to check this if necessary

    dof_handler.distribute_dofs(fe);

    

    

    previous_solution.reinit (dof_handler.n_dofs());

    

    solution_transfer.interpolate(solution, previous_solution);  

    

        DynamicSparsityPattern dsp(dof_handler.n_dofs());

        DoFTools::make_sparsity_pattern (dof_handler, dsp);

        sparsity_pattern.copy_from(dsp);

        

        solution.reinit (dof_handler.n_dofs());

        system_matrix.reinit (sparsity_pattern);

        system_rhs.reinit (dof_handler.n_dofs());

}


I am sure that I am reinitialize matrixes... 

but I am not sure about whether the redline is doing enough for reinitialization of sparsity pattern....

Jaekwang Kim 

2017년 3월 17일 금요일 오전 8시 30분 21초 UTC-5, Wolfgang Bangerth 님의 말:

[Wolfgang Bangerth]

On 03/17/2017 08:03 AM, Jaekwang Kim wrote:
> After the first cycle, I am doing
>
> template<intdim>
>

> voidnonlinear<dim>::refine_grid ()

>
> {
>
>
> Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
>
> KellyErrorEstimator<dim>::estimate (dof_handler,
>
> QGauss<dim-1>(3),
>

> typenameFunctionMap<dim>::type(),

>
> solution,
>
> estimated_error_per_cell);
>
> GridRefinement::refine_and_coarsen_fixed_number (triangulation,
>
> estimated_error_per_cell,
>
> 0.3, 0.0);
>
>
>
>
>
> triangulation.prepare_coarsening_and_refinement ();
>
> SolutionTransfer<dim> solution_transfer(dof_handler);
>
> solution_transfer.prepare_for_coarsening_and_refinement(solution);
>
> triangulation.execute_coarsening_and_refinement(); // I think n_dof() has
> increased here. You need to check this if necessary
>
> dof_handler.distribute_dofs(fe);
>
>
>
>
>
> previous_solution.reinit (dof_handler.n_dofs());
>
>
>
> solution_transfer.interpolate(solution, previous_solution);
>
>
>
> DynamicSparsityPattern dsp(dof_handler.n_dofs());
>
> DoFTools::make_sparsity_pattern (dof_handler, dsp);
>
> sparsity_pattern.copy_from(dsp);
>
>
>
> solution.reinit (dof_handler.n_dofs());
>
> system_matrix.reinit (sparsity_pattern);
>
> system_rhs.reinit (dof_handler.n_dofs());
>
>
> }
>
>
> I am sure that I am reinitialize matrixes...
> but I am not sure about whether the redline is doing enough for
> reinitialization of sparsity pattern....

It looks like it, but it seems like it's not working, and just seeing the
little code snippet is not enough for me to tell for sure. You'll have to
debug what is going on. My first stab would be: At the beginning of your
assembly function, put something like
Assert (system_matrix.n() == dof_handler.n_dofs(), ExcInternalError());
to verify that indeed your matrix has the right size.

[Jaekwang Kim]

Thanks for your response. 

I tried to assert assemble function as you suggested, but I didn't receive any error message. 

What make me more confuse me is the code is working with global refinement.

if I input 100% cells to be refined.. 

 GridRefinement::refine_and_coarsen_fixed_number (triangulation,

                                                     estimated_error_per_cell,

                                                     1.0, 0.0);

if I input 90 cells to be refined, it goes at least one cycle, 

but same or less than 80 cell, I keep receiving error message that... 

Might AMR be related to my program?  

 

An error occurred in line <1668> of file </Users/kimjaekwang/dealii-8.4.1/include/deal.II/lac/constraint_matrix.templates.h> in function

    void dealii::internals::dealiiSparseMatrix::add_value(const LocalType, const size_type, const size_type, SparseMatrixIterator &) [SparseMatrixIterator = dealii::SparseMatrixIterators::Iterator<double, false>, LocalType = double]

The violated condition was: 

    matrix_values->column() == column

The name and call sequence of the exception was:

    typename SparseMatrix<typename SparseMatrixIterator::MatrixType::value_type>::ExcInvalidIndex(row, column)

Additional Information: 

You are trying to access the matrix entry with index <1,305>, but this entry does not exist in the sparsity pattern of this matrix.

The most common cause for this problem is that you used a method to build the sparsity pattern that did not (completely) take into account all of the entries you will later try to write into. An example would be building a sparsity pattern that does not include the entries you will write into due to constraints on degrees of freedom such as hanging nodes or periodic boundary conditions. In such cases, building the sparsity pattern will succeed, but you will get errors such as the current one at one point or other when trying to write into the entries of the matrix.

template <int dim>

void nonlinear<dim>::refine_grid ()

{

    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());

    KellyErrorEstimator<dim>::estimate (dof_handler,

                                        QGauss<dim-1>(3),

                                        typename FunctionMap<dim>::type(),

                                        solution,

                                        estimated_error_per_cell);

    GridRefinement::refine_and_coarsen_fixed_number (triangulation,

                                                     estimated_error_per_cell,

                                                     1.0, 0.0);

    

    /* You need to tranfer solution (n-1 step) domain into previous solution (n step)*/

   

    triangulation.prepare_coarsening_and_refinement ();

    SolutionTransfer<dim> solution_transfer(dof_handler);

    solution_transfer.prepare_for_coarsening_and_refinement(solution);

    triangulation.execute_coarsening_and_refinement(); // I think n_dof() has increased here. You need to check this if necessary

    dof_handler.distribute_dofs(fe);

    

   // Vector<double> tmp(dof_handler.n_dofs()); //tmp - n step dof

    

    previous_solution.reinit (dof_handler.n_dofs());

    

    solution_transfer.interpolate(solution, previous_solution);  // interpolate (input, output)

    

    

    DynamicSparsityPattern dsp(dof_handler.n_dofs());

    DoFTools::make_sparsity_pattern (dof_handler, dsp);

    sparsity_pattern.copy_from(dsp);

    

    solution.reinit (dof_handler.n_dofs());

    system_matrix.reinit (sparsity_pattern);

    system_rhs.reinit (dof_handler.n_dofs());

    

    constraints.clear ();

[Wolfgang Bangerth]

On 03/20/2017 10:58 AM, Jaekwang Kim wrote:
>
> I tried to assert assemble function as you suggested, but I didn't
> receive any error message.
>
> What make me more confuse me is the code is working with global refinement.

Then are you building the sparsity pattern taking into account hanging
node constraints? Take a look how that is done in step-6, for example.

[Jaekwang Kim]

Thank you ! 

You just pointed out exact problem that my code had!

I was not taking into constraints when I make new sparsity pattern! 

  DoFTools::make_sparsity_pattern(dof_handler,

                                    dsp,

                                    constraints,

                                    /*keep_constrained_dofs = */ false);

2017년 3월 20일 월요일 오후 1시 3분 50초 UTC-4, Wolfgang Bangerth 님의 말:

[Jaekwang Kim]

Dear, Dr. Bangerth

Adding up to this, may I ask you one more question?

After this,,, I was trying to transfer my block vector solution through Solutiontransfer class. but I am having difficulty in this. 

For example, I have fluid problem vector value solution (u_x, u_y, pressure)

at refie_mesh function, I tried...

    template <int dim>

    void

    StokesProblem<dim>::refine_mesh () //AMR based on VELOCITY GRADIENT

    {

         const FEValuesExtractors::Vector velocities (0);

            

         Vector<float> estimated_error_per_cell (triangulation.n_active_cells());

        

         KellyErrorEstimator<dim>::estimate (dof_handler,QGauss<dim-1>(degree+1),typename FunctionMap<dim>::type(),

         solution,estimated_error_per_cell, fe.component_mask(velocities));

         

         GridRefinement::refine_and_coarsen_fixed_number (triangulation,estimated_error_per_cell,0.3, 0.0);

                

        triangulation.prepare_coarsening_and_refinement ();

       

        SolutionTransfer<dim> solution_transfer(dof_handler); // Here my doc hander might include all (u,v,p) solution 

After this I reinitialized my previous solution, where my current solution will be interpolated, for newly refined dof_handler. 

And When I tried to transfer my solution.... 

        solution_transfer.prepare_for_coarsening_and_refinement(solution.block(0));

        //since this function only takes vector type inputs... I assume that I cannot transfer whole solution (u,v,p) at once? 

        

        triangulation.execute_coarsening_and_refinement();

        

        

        //Transfer solution

        solution_transfer.interpolate(solution.block(0), previous_solution.block(0));

        

        solution_transfer.prepare_for_coarsening_and_refinement(solution.block(1));

        

        solution_transfer.interpolate(solution.block(1), previous_solution.block(1));

but I run into this error. 

I think it is because there two lines are conflicting....  because they don't have same number of dof_handler...

1. SolutionTransfer<dim> solution_transfer(dof_handler);

2. solution_transfer.prepare_for_coarsening_and_refinement(solution.block(0));

will there be other easy way I can transfer vector-valued solution without having this difficulty? 

I always appreciate your help !

Thanks

Jaekwang Kim 

An error occurred in line <253> of file </Users/kimjaekwang/dealii-8.4.1/source/numerics/solution_transfer.cc> in function

    void dealii::SolutionTransfer<2, dealii::Vector<double>, dealii::DoFHandler<2, 2> >::prepare_for_coarsening_and_refinement(const std::vector<VectorType> &) [dim = 2, VectorType = dealii::Vector<double>, DoFHandlerType = dealii::DoFHandler<2, 2>]

The violated condition was: 

    all_in[i].size()==n_dofs_old

The name and call sequence of the exception was:

    ExcDimensionMismatch(all_in[i].size(),n_dofs_old)

Additional Information: 

Dimension 238 not equal to 274

[Jean-Paul Pelteret]

Dear Jaekwang,

Step-31 demonstrates how to use the SolutionTransfer class with BlockVectors. In short, you need to specify the vector type as a template argument to the class. Here's how its done in that tutorial:

SolutionTransfer<dim,TrilinosWrappers::MPI::Vector>

temperature_trans(temperature_dof_handler); // non-block vector

SolutionTransfer<dim,TrilinosWrappers::MPI::BlockVector>

stokes_trans(stokes_dof_handler); // block vector

triangulation.prepare_coarsening_and_refinement();

temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);

stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);

I hope this helps,

Jean-Paul

[Jaekwang Kim]

thank you for your response. 

I am trying to following this... 

but it seems that I have to use the class,  SolutionTransfer<dim,TirilinoisWarppers::MPI::BlockVector> to  transfer block vector.... 

Without using MPI, isn't it possible to assess this class? like... SolutionTransfer<dim,BlockVector> , while it shows error ,

error: use of

      class template 'BlockVector' requires template arguments

        SolutionTransfer<dim,BlockVector> solution_transfer(dof_handler)...

                             ^~~~~~~~~~~

/Users/kimjaekwang/Programs/dealii/include/deal.II/n

For now I don't have mpi in my computer.....If there is no way to do this without using MPI, i have to install..

Thank you. 

2017년 3월 21일 화요일 오전 4시 25분 37초 UTC-4, Jean-Paul Pelteret 님의 말:

[Jean-Paul Pelteret]

Hi Jaekwang,

No, you've misunderstood in that I was just showing that to you as an example. You don't need to use that exact Vector class. You were nearly at the answer with what you've posted - you just need to read the compiler error more carefully to understand what its telling you is wrong. The deal.II BlockVector class is templatised on a number type, so must pass the arguments along with the vector type when defining the vector template argument of the SolutionTransfer class. So, in your case, this would be

SolutionTransfer<dim,BlockVector<double> > solution_transfer(dof_handler);

This I suspect will fix the problem you've shown above.

Regards,

Jean-Paul 

[Jaekwang Kim]

Thank you!!

Just got what've meant. 

and fixed my code properly

Thanks ! 

Jaekwang Kim 

2017년 3월 21일 화요일 오후 12시 55분 7초 UTC-4, Jean-Paul Pelteret 님의 말: