Solving time dependent heat equation with MPI (PETsc)

[Mark Ma]

Dear experienced deal.II users and developers,

I want to solve a heat equation in the time domain with distributed memory using MPI, but the results are incorrect. In order to do so, I reference tutorial step-23 for time updating method and step-40 for implementing MPI. May I ask whether my boundary condition is right or not? Should we do compress() after apply_boundary_values()? Thanks in advance!

With best regards,
Mark





In PETsc or Trilinos, how could we set the boundary conditions and initial conditions?


Since old_solution_T would be read in other places (not shown here), it is initialized with ghost cells; while old_solutuon_T_cal is only used for written, it does not have ghost cells. see code as following,
//
old_solution_T_cal.reinit (locally_owned_dofs,mpi_communicator);
old_solution_T.reinit (locally_owned_dofs,locally_relevant_dofs,mpi_communicator);
//
..........


  template <int dim>
  void ThermalDiffusion<dim>::run ()
  {


    setup_system();
    assemble_system();

   
    VectorTools::project (dof_handler, constraints, QGauss<dim>(degree),
                          InitialValues_T<dim>(),
                          old_solution_T_cal);

    old_solution_T = old_solution_T_cal;

    LA::MPI::Vector tmp (locally_owned_dofs,mpi_communicator);
    LA::MPI::Vector forcing_terms (locally_owned_dofs,mpi_communicator);

    for (timestep_number=1, time=time_step;
         time<= global_simulation_end_time;
         time+=time_step, ++timestep_number)
      {
        pcout << "Time step " << timestep_number
                  << " at t=" << time
                  << std::endl;


//-----------------------------------
//run to solve T
//-----------------------------------

  //
  //time dependent

        //assign right hand side
        mass_matrix_T.vmult (system_rhs, old_solution_T_cal);

        laplace_matrix_T.vmult (tmp, old_solution_T_cal);
        system_rhs.add (-time_step * (1-theta), tmp);

        assemble_rhs_T (time);
        forcing_terms = dynamic_rhs_T;
        forcing_terms *= time_step * theta;

        assemble_rhs_T (time - time_step);
        tmp = dynamic_rhs_T;
        forcing_terms.add (time_step*(1-theta),tmp);
        system_rhs.add (1,forcing_terms);

        //assign system matrix
        system_matrix.copy_from (mass_matrix_T);
        system_matrix.add (time_step * theta, laplace_matrix_T);



        {
          BoundaryValues_Temperature<dim> boundary_values_function_T;
          //boundary_values_function.set_time (time);

          std::map<types::global_dof_index,double> boundary_values_T;

          VectorTools::interpolate_boundary_values (dof_handler,
                                                    BOUNDARY_NUM,
                                                    boundary_values_function_T,
                                                    boundary_values_T);

          MatrixTools::apply_boundary_values (boundary_values_T,
                                              system_matrix,
                                              solution_T,
                                              system_rhs,
                                              false);


        }
        solve_T ();

        if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 32)
        {
            TimerOutput::Scope t(computing_timer,"output");
            output_results ();
        }

        computing_timer.print_summary ();
        computing_timer.reset();
       
        pcout << std::endl;


        old_solution_T = solution_T;


      }
  }
}

//definition of solve_T()

  template <int dim>
  void ThermalDiffusion<dim>::solve_T ()
  {
    TimerOutput::Scope t(computing_timer,"solve_T");
    LA::MPI::Vector completely_distributed_solution (locally_owned_dofs,mpi_communicator);
    SolverControl   solver_control (dof_handler.n_dofs(),1e-12);

#ifdef USE_PETSC_LA
    LA::SolverCG solver(solver_control,mpi_communicator);
#else
    LA::SolverCG solver(solver_control);
#endif

    LA::MPI::PreconditionAMG preconditioner;
    LA::MPI::PreconditionAMG::AdditionalData    data;

#ifdef USE_PETSC_LA
    data.symmetric_operator = true;
#else
    /* Trilinos defaults are good */
#endif

    preconditioner.initialize(system_matrix,data);

    solver.solve (system_matrix,completely_distributed_solution,system_rhs,preconditioner);

   

    pcout << "  Solved in " << solver_control.last_step()
          << "  iterations."<< std::endl;   

    constraints.distribute (solution_T);
   
    solution_T = completely_distributed_solution;

    pcout << "success...solve_T()..." <<std::endl;
  }

[Wolfgang Bangerth]

On 10/10/2017 08:40 AM, Mark Ma wrote:
>
> I want to solve a heat equation in the time domain with distributed memory
> using MPI, but the results are incorrect. In order to do so, I reference
> tutorial step-23 for time updating method and step-40 for implementing MPI.
> May I ask whether my boundary condition is right or not? Should we do
> compress() after apply_boundary_values()? Thanks in advance!

Jack -- how exactly is your solution wrong when you look at it? Do the
boundary values look wrong? Are they correct if you run your MPI program with
just one MPI process?

In general, using MatrixTools::apply_boundary_values() is not the way to go
with MPI programs. Rather, use a ConstraintMatrix and incorporate the boundary
values into the same object as you do with hanging node constraints. That's
what all of the parallel programs do, if I recall correctly.

Best
W.

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

[Mark]

Hi Prof. Wolfgang Bangerth,

I really appreciate for your quick reply. Now I do as what you said, using VectorTools in setup_system(). It seems working, but there is still one question that I always find very large values at points near boundaries, see attached figures. Is there a way to solve this? Thanks for advance!

Have a nice day,
Mark


Code used as follows,

constraints.clear();
constraints.reinit(locally_relevant_dofs);
VectorTools::interpolate_boundary_values(dof_handler, BOUNDARY_NUM, BoundaryValues_T<dim>(),constraints);
constraints.close();

Initial value projected by

    VectorTools::project (dof_handler, constraints, QGauss<dim>(degree),
                          InitialValues_T<dim>(),
                          old_solution_T_cal);

After 30us, Temperature distribution

在 2017年10月10日星期二 UTC+2下午5:48:57，Wolfgang Bangerth写道：

[Wolfgang Bangerth]

On 10/10/2017 11:40 AM, Mark wrote:
>
> I really appreciate for your quick reply. Now I do as what you said,
> using VectorTools in setup_system(). It seems working, but there is
> still one question that I always find very large values at points near
> boundaries, see attached figures. Is there a way to solve this?

There are so many possible reasons for this that I can't even speculate.
Play with what you have to find out what the common theme is -- how do
things change if you try different boundary conditions, initial
conditions, right hand sides, time steps, etc. What affects the problem
and what doesn't? Can you make the problem simpler (e.g., use zero
boundary values and zero initial conditions)? Does this only happen in
parallel or also if you run the program on a single processor?

This is sort of the mental list that I go through when debugging these
kinds of problems. I try to isolate what it is that affects the problem
and make the problem as simple as possible. Whatever is left must then
at least be related to the cause.

[Mark Ma]

Hi Prof. Wolfgang Bangerth,

Thanks for your useful advice. Finally I solved this problem.

Best,

Mark

                           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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/5hC7jODg-7k/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Wolfgang Bangerth]

On 10/12/2017 03:02 AM, Mark Ma wrote:
>
> Thanks for your useful advice. Finally I solved this problem.

So what was the problem? Maybe we can learn from your example?

[Mark Ma]

It is actually my stupid mistake in solvercontrol,

I used

SolverControl solver_control (dof_handler.n_dofs(),1e-12)

This works when the geometry size in the order of 1, but fails at 1e-6, or even 1e-9. Here I actually want to make a micrometer or nanometer size structure.

later, I changed the control into

SolverControl solver_control (5*system_rhs.size(),1e-12*system_rhs.l2_norm())

this works well for structure size of um or nm. I think previous setting may lead to a loss of precision so that the results are always incorrect.

Thanks for your ideas and have a nice day,
Mark

[Wolfgang Bangerth]

On 10/13/2017 02:06 AM, Mark Ma wrote:
>
> later, I changed the control into
> |
> |

> SolverControlsolver_control (5*system_rhs.size(),1e-12*system_rhs.l2_norm())

> |
> |
> this works well for structure size of um or nm. I think previous setting may
> lead to a loss of precision so that the results are always incorrect.

Yes, indeed -- using a tolerance relative to the size of the rhs vector is
important.

[Lucas Campos]

Dear Bangerth,

When you mention 

In general, using MatrixTools::apply_boundary_values() is not the way to go 
with MPI programs. Rather, use a ConstraintMatrix and incorporate the boundary 
values into the same object as you do with hanging node constraints.

This is the way to go due to correctness, or in the sense of scalability?

Bests,

Lucas

[Wolfgang Bangerth]

On 10/13/2017 08:39 AM, Lucas Campos wrote:
>
> In general, using MatrixTools::apply_boundary_values() is not the way to go
> with MPI programs. Rather, use a ConstraintMatrix and incorporate the
> boundary
> values into the same object as you do with hanging node constraints.
>
>
> This is the way to go due to correctness, or in the sense of scalability?

MatrixTools::apply_boundary_values() needs access to the elements of the
matrix because it wants to modify elements after they have already been
written into the matrix. That is already difficult if the matrix is owned by
PETSc or Trilinos -- we can get access to these elements, but it is not
efficient to do so.

But the bigger issue is that the function wants to access elements not only
for the rows of constrained DoFs, but also for the columns. That means that
you may have to access elements that are actually stored on other processors
-- something that can not be done efficiently. Consequently,
MatrixTools::apply_boundary_values() does not attempt to eliminate columns of
the matrix, and you will end up with a non-symmetric matrix even if your
problem is symmetric.

It is better to use the approach via ConstraintMatrix that deals with entries
before they even get into the matrix (and therefore in particular before
matrix entries are sent to other processors).

[Lucas Campos]

Dear Wolfgang,

Thank you for your explanation. Currently, I am using a code that was not written by me, and uses the MatrixTools::apply_boundary_values() approach. I try to change it to use the ConstraintMatrix one. For that, step-40 seems to be the best starting point, like Mark did. 

Thanks again,

Lucas

                           www: http://www.math.colostate.edu/~bangerth/

[Mark Ma]

Hi Lucas,

I am currently studying of this, if you have some problems or very exciting process,  we could discuss if you would like to share. :)

Best,

Mark

在 2017年10月13日星期五 UTC+2下午4:50:02，Lucas Campos写道：

To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.

[Mark Ma]

Prof. Wolfgang Bangerth,

Now the projection of initial values (rewrite the code by manually assemble the matrix and system_rhs and calculate) run OK, but the time updating of T is not correct, same phenomenon appears. I believe this may arise from the fact that direct using matrix vmult (i.e. 

mass_matrix_T.vmult (system_rhs, old_solution_T_cal);

) instead of assembling and distribute_local_to_global again may ignore eliminating the constraint points in matrix or vector, when using constriantMatrix and interporate_boundary_values to apply the boundary condition, I am now checking of this. 

Is there a simple way to update the RHS of old value using something simple like vmult?

Best,

Mark

//---------------------------------------
//Updating of T in time domain------the old way, working in non-mpi version
//---------------------------------------

  //
  //time dependent

        //assign right hand side
        mass_matrix_T.vmult (system_rhs, old_solution_T_cal);

        laplace_matrix_T.vmult (tmp, old_solution_T_cal);
        system_rhs.add (-time_step * (1-theta), tmp);

        assemble_rhs_T (time);
        forcing_terms = dynamic_rhs_T;
        forcing_terms *= time_step * theta;

        assemble_rhs_T (time - time_step);
        tmp = dynamic_rhs_T;
        forcing_terms.add (time_step*(1-theta),tmp);
        system_rhs.add (1,forcing_terms);

        //assign system matrix
        system_matrix.copy_from (mass_matrix_T);
        system_matrix.add (time_step * theta, laplace_matrix_T);

在 2017年10月13日星期五 UTC+2下午4:43:30，Wolfgang Bangerth写道：

[Mark Ma]

Prof. W. Bangerth,

Please find attached source code, where the projection part is ok now, but time updating of solution part results is incorrect I believe (see boundary results). This confused me a few weeks and still did not find a solution. 

Best,

Mark

在 2017年10月13日星期五 UTC+2下午2:55:04，Wolfgang Bangerth写道：

[Wolfgang Bangerth]

On 10/15/2017 12:41 PM, Mark Ma wrote:
>
> Now the projection of initial values (rewrite the code by manually assemble
> the matrix and system_rhs and calculate) run OK, but the time updating of T is
> not correct, same phenomenon appears. I believe this may arise from the fact
> that direct using matrix vmult (i.e.
> |
> mass_matrix_T.vmult (system_rhs,old_solution_T_cal);
> |
>
> ) instead of assembling and distribute_local_to_global again may ignore
> eliminating the constraint points in matrix or vector, when using
> constriantMatrix and interporate_boundary_values to apply the boundary
> condition, I am now checking of this.

So when you visualize the solution, the error is at the boundary but it looks
correct in the interior?

> Is there a simple way to update the RHS of old value using something simple
> like vmult?

What is the equation you are trying to implement?

[Mark Ma]

On 10/15/2017 12:41 PM, Mark Ma wrote:
>
> Now the projection of initial values (rewrite the code by manually assemble
> the matrix and system_rhs and calculate) run OK, but the time updating of T is
> not correct, same phenomenon appears. I believe this may arise from the fact
> that direct using matrix vmult (i.e.
> |
> mass_matrix_T.vmult (system_rhs,old_solution_T_cal);
> |
>
> ) instead of assembling and distribute_local_to_global again may ignore
> eliminating the constraint points in matrix or vector, when using
> constriantMatrix and interporate_boundary_values to apply the boundary
> condition, I am now checking of this.

So when you visualize the solution, the error is at the boundary but it looks
correct in the interior?

Yes, it is. Thereafter, the error at boundaries does propagate and then interfere with the interior.

> Is there a simple way to update the RHS of old value using something simple
> like vmult?

What is the equation you are trying to implement?

I am trying just the Heat equation,

rho * C * dT/dt = nabla *(kappba nalbaT) + f(x,y,t);

Boundary condition is,

T(x,y,z,t) = 300K, at boundaries.

[Mark Ma]

Projection of initial value, time 000

time 001

time 002

在 2017年10月10日星期二 UTC+2下午4:40:34，Mark Ma写道：

[Wolfgang Bangerth]

> So when you visualize the solution, the error is at the boundary but it looks
> correct in the interior?
>
> Yes, it is. Thereafter, the error at boundaries does propagate and then
> interfere with the interior.

From the pictures you posted in a follow-up, it looks like the boundary
values are actually correct, but that the problem starts one layer of cells
into the domain. Do I see this correctly?

> > Is there a simple way to update the RHS of old value using something
> simple
> > like vmult?
>
> What is the equation you are trying to implement?
>
> I am trying just the Heat equation,
> rho * C * dT/dt = nabla *(kappba nalbaT) + f(x,y,t);
> Boundary condition is,
> T(x,y,z,t) = 300K, at boundaries.

I meant the "update the RHS of old value" equation.

[Mark Ma]

在 2017年10月16日星期一 UTC+2下午4:31:25，Wolfgang Bangerth写道：

>     So when you visualize the solution, the error is at the boundary but it looks
>     correct in the interior?
>
> Yes, it is. Thereafter, the error at boundaries does propagate and then
> interfere with the interior.

 From the pictures you posted in a follow-up, it looks like the boundary
values are actually correct, but that the problem starts one layer of cells
into the domain. Do I see this correctly?

Yes, exactly. And if I use the BC values as 0, then it seems OK now (see figs below). What if I want to apply some value other than 0, this problem seems annoying.

 

>      > Is there a simple way to update the RHS of old value using something
>     simple
>      > like vmult?
>
>     What is the equation you are trying to implement?
>
> I am trying just the Heat equation,
> rho * C * dT/dt = nabla *(kappba nalbaT) + f(x,y,t);
> Boundary condition is,
> T(x,y,z,t) = 300K, at boundaries.
I meant the "update the RHS of old value" equation.
 

sorry, I mistook your idea. The RHS function is,

f(x,y,z,t) = I0 * exp(-((x-v*t)^2 + (y-v*t)^2 +(z-v*t)^2)/2c^2), which stands for absorption of the laser.

Boundary value = 0, initial values = 0

[Wolfgang Bangerth]

On 10/16/2017 08:55 AM, Mark Ma wrote:
>
> > So when you visualize the solution, the error is at the boundary
> but it looks
> > correct in the interior?
> >
> > Yes, it is. Thereafter, the error at boundaries does propagate and then
> > interfere with the interior.
>
> From the pictures you posted in a follow-up, it looks like the boundary
> values are actually correct, but that the problem starts one layer of cells
> into the domain. Do I see this correctly?
>
>
> Yes, exactly. And if I use the BC values as 0, then it seems OK now (see figs
> below). What if I want to apply some value other than 0, this problem seems
> annoying.

Yes, bugs are highly annoying indeed :-)

It almost looks to me like you're applying both Dirichlet values (via the
ConstraintMatrix) and Neumann boundary values (via a boundary integral). But
then I haven't taken a look at the code, so I can't really say for sure.

[Mark Ma]

在 2017年10月16日星期一 UTC+2下午5:05:20，Wolfgang Bangerth写道：

On 10/16/2017 08:55 AM, Mark Ma wrote:
>
>      >     So when you visualize the solution, the error is at the boundary
>     but it looks
>      >     correct in the interior?
>      >
>      > Yes, it is. Thereafter, the error at boundaries does propagate and then
>      > interfere with the interior.
>
>       From the pictures you posted in a follow-up, it looks like the boundary
>     values are actually correct, but that the problem starts one layer of cells
>     into the domain. Do I see this correctly?
>
>
> Yes, exactly. And if I use the BC values as 0, then it seems OK now (see figs
> below). What if I want to apply some value other than 0, this problem seems
> annoying.

Yes, bugs are highly annoying indeed :-)

It almost looks to me like you're applying both Dirichlet values (via the
ConstraintMatrix) and Neumann boundary values (via a boundary integral). But
then I haven't taken a look at the code, so I can't really say for sure.

I think I only applied Dirichlet values via ConstraintMatrix. There is no boundary integrals, so I believe Neumann BC is not implemented.

But may be you are right since for most of the case, heat equations are solved assuming an adiabatic process, that is dT/dx = 0 at boundaries.

This zero values is automatically satisfied.

I think the problem comes from the way of setting Boundary values using constraint method, which in dealii it makes some kind of constraints at vertex of boundaries, i.e. x22 = x1/2 + x2/2, after elliminating the value x22 as setting x22=0, then x1 = -x2.  So, if all values are zero at initial, this would not be a problem. Otherwise, you would see this problem that at vertex point close to boundary, the value is negative to the boundary. This is what I could think of to my best knowledge of deal.II. Maybe I need clear the constraints before solve() function.

Best,

Mark

[Wolfgang Bangerth]

On 10/16/2017 09:23 AM, Mark Ma wrote:
>
> It almost looks to me like you're applying both Dirichlet values (via
> the ConstraintMatrix) and Neumann boundary values (via a boundary
> integral). But then I haven't taken a look at the code, so I can't
> really say for sure.
>
>
> I think I only applied Dirichlet values via ConstraintMatrix. There
> is no boundary integrals, so I believe Neumann BC is not
> implemented.

My guess came from the fact that a Neumann boundary condition would act
like a boundary heat source. It looks like you have a heat source in the
last layer of cells around the boundary.

I don't know whether that is true, but it's a place I'd try to
investigate in debugging.

> But may be you are right since for most of the case, heat equations
> are solved assuming an adiabatic process, that is dT/dx = 0 at
> boundaries. This zero values is automatically satisfied.

That's a different question. Whether one may *physically* want a
different equation is besides the point -- you have chosen one equation
you'd like to solve, but the numerical solution is wrong. It's
worthwhile trying to understand why that is before thinking about
*other* equations.

> I think the problem comes from the way of setting Boundary values
> using constraint method, which in dealii it makes some kind of
> constraints at vertex of boundaries, i.e. x22 = x1/2 + x2/2,

Hm, that would be the constraint of a hanging node. Do you have hanging
nodes in your problem? If you do, may I suggest to simplify the problem
and see if you have the same problem on uniform meshes?

[Mark Ma]

On 10/16/2017 09:23 AM, Mark Ma wrote:
>
> It almost looks to me like you're applying both Dirichlet values (via
> the ConstraintMatrix) and Neumann boundary values (via a boundary
> integral). But then I haven't taken a look at the code, so I can't
> really say for sure.
>
>
> I think I only applied Dirichlet values via ConstraintMatrix. There
> is no boundary integrals, so I believe Neumann BC is not
> implemented.

My guess came from the fact that a Neumann boundary condition would act
like a boundary heat source. It looks like you have a heat source in the
last layer of cells around the boundary.

I don't know whether that is true, but it's a place I'd try to
investigate in debugging.

Yes, it looks like so. But then I comment the time updating code of rhs source by the following code, this boundary heat still appears(see figs below)

 //
  //time dependent

        //assign right hand side
        mass_matrix_T.vmult (system_rhs, old_solution_T_cal);
        laplace_matrix_T.vmult (tmp, old_solution_T_cal);
        system_rhs.add (-time_step * (1-theta), tmp);

/*

        assemble_rhs_T (time);
        forcing_terms = dynamic_rhs_T;
        forcing_terms *= time_step * theta;

        assemble_rhs_T (time - time_step);
        tmp = dynamic_rhs_T;
        forcing_terms.add (time_step*(1-theta),tmp);
        system_rhs.add (1,forcing_terms);

*/

        //assign system matrix
        system_matrix.copy_from (mass_matrix_T);
        system_matrix.add (time_step * theta, laplace_matrix_T);

I think this problem lies in the time updating of solution using old_solution, since the mass_matrix and laplace_matrix have already eliminated the constraint node, mass_matrix_T.vmult (system_rhs, old_solution_T_cal); is no longer valid for this case. This may result into ZEROs at all constrained nodes, so at boundaries these zero values then propagate into center, as can be seen from the figure below. For time updating of solution using old_solution, maybe tutorial-32 could give some ideas. However to completely understand step-32 may spend some more time for me. 

 

> But may be you are right since for most of the case, heat equations
> are solved assuming an adiabatic process, that is dT/dx = 0 at
> boundaries. This zero values is automatically satisfied.

That's a different question. Whether one may *physically* want a
different equation is besides the point -- you have chosen one equation
you'd like to solve, but the numerical solution is wrong. It's
worthwhile trying to understand why that is before thinking about
*other* equations.


> I think the problem comes from the way of setting Boundary values
> using constraint method, which in dealii it makes some kind of
> constraints at vertex of boundaries, i.e. x22 = x1/2 + x2/2,

Hm, that would be the constraint of a hanging node. Do you have hanging
nodes in your problem? If you do, may I suggest to simplify the problem
and see if you have the same problem on uniform meshes?

there is no handing node in my code since I used a uniformed mesh. And also I tried this, results show no difference.

I am really appreciate your patience to look for this question during your busy schedule.

Thanks and have a nice day,

Mark

-------------------------------------------------

Laboratoire Hubert Curien, UMR CNRS 5516,
Bâtiment F 18 Rue du Professeur Benoît Lauras
42000 Saint-Etienne
FRANCE

[Mark Ma]

source code are attached.

在 2017年10月16日星期一 UTC+2下午9:14:13，Wolfgang Bangerth写道：

[Wolfgang Bangerth]

On 10/16/2017 01:53 PM, Mark Ma wrote:
>
> I think this problem lies in the time updating of solution using
> old_solution, since the mass_matrix and laplace_matrix have already

> eliminated the constraint node, /*mass_matrix_T.vmult
> (system_rhs, old_solution_T_cal*//*);*/ is no longer valid for this
> case.

Yes, this sounds reasonable from your description. Have you looked at
steps 24, 25, 26 to see how it is done there? I could imagine that you
need to think about which degrees of freedom you want eliminated/not
eliminated in the matrices you multiply with. To this end, write out the
weak form of the problem you need to solve in each time step, and think
specifically about the boundary values of the trial and test functions
and whether they should be part of the matrix you use on the right hand
side or not. It is possible that you need to use a different matrix for
the lhs and rhs operations.

[Mark Ma]

Hi Wlofgang,

Good news, finally I resolved this problem by assembling the left_Mass_matrix and right_rhs by the code following,

              local_rho_c_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
                                            phi_i_u *
                                            phi_j_u *
                                            fe_values.JxW (q))
                                            +
                                            time_step_00 * (1-theta) *
                                            (kappa_fun_T.value(fe_values.quadrature_point(q)) *
                                            div_phi_i_u *
                                            div_phi_j_u *
                                            fe_values.JxW (q));

              local_kappa_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
                                            phi_i_u *
                                            phi_j_u *
                                            fe_values.JxW (q))
                                            +
                                            time_step_00 * (theta) *
                                            (kappa_fun_T.value(fe_values.quadrature_point(q)) *
                                            div_phi_i_u *
                                            div_phi_j_u *
                                            fe_values.JxW (q));

.....
--------------------------------------------------------------------------
      constraints.distribute_local_to_global(local_rho_c_T_matrix,
                                             local_dof_indices,
                                             mass_matrix_T);

      constraints.distribute_local_to_global(local_kappa_T_matrix,
                                             local_dof_indices,
                                             laplace_matrix_T);


And in the run() function,

 for (timestep_number=1, time=time_step;
         time<= global_simulation_end_time;
         time+=time_step, ++timestep_number)
      {
        pcout << "Time step " << timestep_number
                  << " at t=" << time
                  << std::endl;

       

        mass_matrix_T.vmult (system_rhs_T, old_solution_T_cal);

        assemble_rhs_T (time);
        forcing_terms = dynamic_rhs_T;
        forcing_terms *= time_step * theta;

        assemble_rhs_T (time - time_step);
        tmp = dynamic_rhs_T;
        forcing_terms.add (time_step*(1-theta),tmp);

        system_rhs_T.add(1,forcing_terms);


        system_matrix_T.copy_from (laplace_matrix_T);


        solve_T_run ();     // solve and store solution to new_solution_T

        //
        old_solution_T = new_solution_T;        // used for output, GHOST CELLS
        old_solution_T_cal = new_solution_T;

This problem seems coming from the assemble_system() that I could not tell why. But confused me a few weeks, it finally turned out from here (previous code),

              local_rho_c_T_matrix(i,j) += (rho_C_fun_T.value(fe_values.quadrature_point(q)) *
                                            phi_i_u *
                                            phi_j_u *
                                            fe_values.JxW (q));

              local_kappa_T_matrix(i,j) += (kappa_fun_T.value(fe_values.quadrature_point(q)) *
                                            div_phi_i_u *
                                            div_phi_j_u *
                                            fe_values.JxW (q));
.....
--------------------------------------------------------------------------
      constraints.distribute_local_to_global(local_rho_c_T_matrix,
                                             local_dof_indices,
                                             mass_matrix_T);

      constraints.distribute_local_to_global(local_kappa_T_matrix,
                                             local_dof_indices,
                                             laplace_matrix_T);

and in run(),

    for (timestep_number=1, time=time_step;
         time<= global_simulation_end_time;
         time+=time_step, ++timestep_number)
      {
        pcout << "Time step " << timestep_number
                  << " at t=" << time
                  << std::endl;

       
        mass_matrix_T.vmult (system_rhs_T, old_solution_T_cal);
        laplace_matrix_T.vmult (tmp,old_solution_T_cal);
        system_rhs_T.add(-time_step * (1-theta),tmp);

        assemble_rhs_T (time);
        forcing_terms = dynamic_rhs_T;
        forcing_terms *= time_step * theta;

        assemble_rhs_T (time - time_step);
        tmp = dynamic_rhs_T;
        forcing_terms.add (time_step*(1-theta),tmp);

        system_rhs_T.add(1,forcing_terms);

                // system_matrix_T
        system_matrix_T.copy_from (mass_matrix_T);
        system_matrix_T.add (time_step * theta, laplace_matrix_T);


        solve_T_run ();     // solve and store solution to new_solution_T

        //
        old_solution_T = new_solution_T;        // used for output, GHOST CELLS
        old_solution_T_cal = new_solution_T;

        output_results (timestep_number);


      }

Finally, the results

Best,
Mark


PS: I think this could be another tutorial show how to implement this, I would like to contribute if you think it is interesting for Deal.II project.

在 2017年10月17日星期二 UTC+2上午2:31:15，Wolfgang Bangerth写道：

[Wolfgang Bangerth]

On 10/17/2017 07:21 AM, Mark Ma wrote:
>
> PS: I think this could be another tutorial show how to implement this, I would
> like to contribute if you think it is interesting for Deal.II project.

Yes, I think that would be great. Or, maybe simpler and with less
documentation requirements: make it a part of the code gallery:
http://dealii.org/code-gallery.html
Want to do that?

Thanks in advance!
Wolfgang

[Mark Ma]

Sure, no problem.

                           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 a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/5hC7jODg-7k/unsubscribe.

To unsubscribe from this group and all its topics, send an email to dealii+unsubscribe@googlegroups.com.

[Wolfgang Bangerth]

On 10/17/2017 07:53 AM, Mark Ma wrote:
> Sure, no problem.

Excellent -- let us know if you need help! We're eager to grow the code
gallery!

Best

[syed ansari]

Hello all,

            I tried to find this code in the deal.ii code gallery but I am unable to locate it. I request the group members to help me to find this code.

Thanks in advance.

Best Regards,

Syed Ansari S.

[Wolfgang Bangerth]

On 1/22/22 2:56 AM, syed ansari wrote:
>             I tried to find this code in the deal.ii code gallery but I am
> unable to locate it. I request the group members to help me to find this code.

I don't think it was ever submitted.

Best
W.

[syed ansari]

Thank you very much Prof. Bangerth for your kind reply.

--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/1e09ca9d-c53d-605b-069d-26d35ea0c25a%40colostate.edu.

[Mark Ma]

Hi,

The code was submitted to the github (https://github.com/dealii/code-gallery/pull/92) but is waiting for checking. Alternatively, you can find the code at my github at:

https://github.com/MarkMa1990/code-gallery/tree/master/Distributed_Moving_Laser_Heating

Have a nice day,

Mark

syed ansari <ansa...@gmail.com> 于2022年1月24日周一 09:35写道：

You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/5hC7jODg-7k/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CANaEQub-80n7v-4wDArz_6EyHfzs_8XzxTZvXc3CncNxbBmMWg%40mail.gmail.com.

[syed ansari]

Thank you very much, Mark, for sharing. This is indeed very helpful for me.

Best Regards,

Syed Ansari S.

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CANwa5smXTP_vK5Tb0jBg_zwnC_hBe6PeHkgGLR%2B2KMJX-d-xyA%40mail.gmail.com.

[Mark Ma]

that's good. feel free to discuss any questions related on this topic, I will get back to you as soon as possible. 

syed ansari <ansa...@gmail.com> 于2022年1月25日周二 14:41写道：

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CANaEQuYhXp1KWb4bxinrS8PqGtKSQdgdNO%2BduHXwDBWUw-kJLw%40mail.gmail.com.

[Wolfgang Bangerth]

> The code was submitted to the github
> (https://github.com/dealii/code-gallery/pull/92

Oh, doh, my apologies for forgetting about it! I hope to get around to
reviewing it soon!