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
<https://lh3.googleusercontent.com/-2TadguSBFTw/WeYDTK6Do_I/AAAAAAAAAD4/h7iKKtK4G8o_96qDO84gDU4niUsxGuQZwCLcBGAs/s1600/heat-equation-mpi.gif>
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写道：
>
> 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.
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> <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 [email protected].
For more options, visit https://groups.google.com/d/optout.
