Hi all,
I have to solve the continuity equation in 1D using Implicit Runge Kutta
method.
The formula I made is this: d(U_t+1)/dt = [-grad(\Gamma_e)+R_i]/exp(U_t)
\Gamma_e ,R_i : Load the value calculated by another function.
U_t : Solution calculated at the previous time.
When implicit_runge_kutta.evolve_one_time_step(..) and
evaluate_diffusion_energy(..) is repeated, the value(U_t) I load changes.
U_t, the value that should be divided, becomes 0 and an error occurs.
I do not know why this value changes. I just want to get a fixed value.
template<int dim>
Vector<double> Electron<dim>::evaluate_diffusion_energy (const Vector<double>
&y)
{
Vector<double> tmp(dof_handler.n_dofs());
tmp = 0.;
system_matrix.vmult(tmp,y); //tmp=(A-B-C)*y
const QGauss<dim> quadrature_formula(fe_degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values |update_gradients|
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector<Tensor<1, dim>> grad_V(n_q_points);
std::vector<Tensor<1, dim>> grad_energyflux(n_q_points);
Vector<double> cell_source(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell->get_dof_indices(local_dof_indices);
cell_source = 0.;
fe_values.reinit(cell);
fe_values.get_function_gradients(Potential,grad_V);
fe_values.get_function_gradients(EnergyFlux,grad_energyflux);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
double E = coeffi_E(local_dof_indices,grad_V[q_point]);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_source(i) += fe_values.JxW(q_point)*fe_values.shape_value(i,
q_point)* // phi_i(x)
(E-grad_energyflux[i][0])/exp(initial_EnergyDensity[local_dof_indices[i]]);
}
constraints_2.distribute_local_to_global(cell_source,
local_dof_indices,
tmp); //(A-B-C)*y+E
}
Vector<double> value(dof_handler.n_dofs());
inverse_mass_matrix.vmult(value,tmp); //value=M^-1*tmp
return value;
}
I look forward to hearing from you, thank you.
Best regards,
sora
--
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].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/f7655830-f6f6-4cbc-a592-254232d93590%40googlegroups.com.
