Skip to content

Commit ffe3dfb

Browse files
YiminJinYiminJin
committed
Fix Newton derivatives in GMG solver for compressible models
1 parent 3ef7b3f commit ffe3dfb

File tree

3 files changed

+69
-33
lines changed

3 files changed

+69
-33
lines changed

‎include/aspect/simulator/solver/matrix_free_operators.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -122,6 +122,11 @@ namespace aspect
122122
*/
123123
bool apply_stabilization_free_surface_faces;
124124

125+
/**
126+
* If true, average the Newton factors in each cell.
127+
*/
128+
bool average_newton_factors;
129+
125130
/**
126131
* Table which stores viscosity values for each cell.
127132
*

‎source/simulator/solver/matrix_free_operators.cc‎

Lines changed: 36 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -196,7 +196,8 @@ namespace aspect
196196
VectorizedArray<number> deta_deps_times_sym_grad_u(0.);
197197
VectorizedArray<number> eps_times_sym_grad_u(0.);
198198
VectorizedArray<number> deta_dp_times_p(0.);
199-
if (cell_data->enable_newton_derivatives)
199+
if (cell_data->enable_newton_derivatives
200+
&& cell_data->average_newton_factors)
200201
{
201202
SymmetricTensor<2,dim,VectorizedArray<number>> sym_grad_u;
202203
VectorizedArray<number> val_p;
@@ -248,6 +249,15 @@ namespace aspect
248249
// Add the Newton derivatives if required.
249250
if (cell_data->enable_newton_derivatives)
250251
{
252+
if(!cell_data->average_newton_factors)
253+
{
254+
deta_deps_times_sym_grad_u = cell_data->newton_factor_wrt_strain_rate_table(cell,q)
255+
* sym_grad_u;
256+
deta_dp_times_p = cell_data->newton_factor_wrt_pressure_table(cell,q) * val_p;
257+
if (cell_data->symmetrize_newton_system)
258+
eps_times_sym_grad_u = cell_data->strain_rate_table(cell,q) * sym_grad_u;
259+
}
260+
251261
velocity_terms +=
252262
( cell_data->symmetrize_newton_system ?
253263
( cell_data->strain_rate_table(cell,q) * deta_deps_times_sym_grad_u +
@@ -256,15 +266,32 @@ namespace aspect
256266
+
257267
2. * cell_data->strain_rate_table(cell,q) * deta_dp_times_p;
258268

259-
if (cell_data->enable_prescribed_dilation)
269+
if (cell_data->is_compressible ||
270+
cell_data->enable_prescribed_dilation)
260271
{
261-
pressure_terms += ( ( cell_data->dilation_derivative_wrt_strain_rate_table(cell,q)
262-
* sym_grad_u )
263-
+
264-
( cell_data->dilation_derivative_wrt_pressure_table(cell,q)
265-
* cell_data->pressure_scaling * val_p )
266-
)
267-
* cell_data->pressure_scaling;
272+
const number one_third = 1.0 / 3.0;
273+
if (cell_data->symmetrize_newton_system)
274+
{
275+
velocity_terms -= one_third * cell_data->newton_factor_wrt_strain_rate_table(cell,q) * (div_u * div_u);
276+
for (unsigned int d = 0; d < dim; ++d)
277+
velocity_terms[d][d] -= one_third * div_u * deta_deps_times_sym_grad_u;
278+
}
279+
else
280+
{
281+
for (unsigned int d = 0; d < dim; ++d)
282+
velocity_terms[d][d] -= 2.0 * one_third * div_u * deta_deps_times_sym_grad_u;
283+
}
284+
285+
if (cell_data->enable_prescribed_dilation)
286+
{
287+
pressure_terms += ( ( cell_data->dilation_derivative_wrt_strain_rate_table(cell,q)
288+
* sym_grad_u )
289+
+
290+
( cell_data->dilation_derivative_wrt_pressure_table(cell,q)
291+
* cell_data->pressure_scaling * val_p )
292+
)
293+
* cell_data->pressure_scaling;
294+
}
268295
}
269296
}
270297

‎source/simulator/solver/stokes_matrix_free_local_smoothing.cc‎

Lines changed: 28 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -454,12 +454,14 @@ namespace aspect
454454
active_cell_data.enable_newton_derivatives = (Parameters<dim>::is_defect_correction(this->get_parameters().nonlinear_solver)
455455
&& this->get_newton_handler().parameters.newton_derivative_scaling_factor != 0);
456456
active_cell_data.enable_prescribed_dilation = this->get_parameters().enable_prescribed_dilation;
457+
active_cell_data.average_newton_factors = (this->get_parameters().material_averaging != MaterialModel::MaterialAveraging::none);
457458

458459
// TODO: these are not implemented yet
459460
for (unsigned int level=0; level<n_levels; ++level)
460461
{
461462
level_cell_data[level].enable_newton_derivatives = false;
462463
level_cell_data[level].enable_prescribed_dilation = false;
464+
level_cell_data[level].average_newton_factors = false;
463465
}
464466

465467
FEValues<dim> fe_values (this->get_mapping(),
@@ -556,15 +558,8 @@ namespace aspect
556558

557559
for (unsigned int q=0; q<n_q_points; ++q)
558560
{
559-
// use the correct strain rate for the Jacobian
560-
// when elasticity is enabled use viscoelastic strain rate
561-
// when stabilization is enabled, use the deviatoric strain rate because the SPD factor
562-
// that is computed is only safe for the deviatoric strain rate (see PR #5580 and issue #5555)
563-
SymmetricTensor<2,dim> effective_strain_rate = in.strain_rate[q];
564-
if (elastic_out != nullptr)
565-
effective_strain_rate = elastic_out->viscoelastic_strain_rate[q];
566-
else if ((this->get_newton_handler().parameters.velocity_block_stabilization & Newton::Parameters::Stabilization::PD) != Newton::Parameters::Stabilization::none)
567-
effective_strain_rate = Utilities::Tensors::consistent_deviator(effective_strain_rate);
561+
const SymmetricTensor<2,dim> effective_strain_rate
562+
= (elastic_out == nullptr ? in.strain_rate[q] : elastic_out->viscoelastic_strain_rate[q]);
568563

569564
// use the spd factor when the stabilization is PD or SPD.
570565
const double alpha = (this->get_newton_handler().parameters.velocity_block_stabilization
@@ -579,14 +574,19 @@ namespace aspect
579574
1.0;
580575

581576
active_cell_data.newton_factor_wrt_pressure_table(cell,q)[i]
582-
= derivatives->viscosity_derivative_wrt_pressure[q] *
583-
derivatives->viscosity_derivative_averaging_weights[q] *
584-
newton_derivative_scaling_factor;
577+
= newton_derivative_scaling_factor * derivatives->viscosity_derivative_wrt_pressure[q]
578+
* (active_cell_data.average_newton_factors ? derivatives->viscosity_derivative_averaging_weights[q] : 1.0);
585579
Assert(std::isfinite(active_cell_data.newton_factor_wrt_pressure_table(cell,q)[i]),
586-
ExcMessage("active_cell_data.newton_factor_wrt_pressure_table is not finite: " + std::to_string(active_cell_data.newton_factor_wrt_pressure_table(cell,q)[i]) +
587-
". Relevant variables are derivatives->viscosity_derivative_wrt_pressure[q] = " + std::to_string(derivatives->viscosity_derivative_wrt_pressure[q]) +
588-
", derivatives->viscosity_derivative_averaging_weights[q] = " + std::to_string(derivatives->viscosity_derivative_averaging_weights[q]) +
589-
", and newton_derivative_scaling_factor = " + std::to_string(newton_derivative_scaling_factor)));
580+
ExcMessage("active_cell_data.newton_factor_wrt_pressure_table is not finite: "
581+
+ std::to_string(active_cell_data.newton_factor_wrt_pressure_table(cell,q)[i])
582+
+ ". Relevant variables are derivatives->viscosity_derivative_wrt_pressure[q] = "
583+
+ std::to_string(derivatives->viscosity_derivative_wrt_pressure[q])
584+
+ (active_cell_data.average_newton_factors ?
585+
", derivatives->viscosity_derivative_averaging_weights[q] = "
586+
+ std::to_string(derivatives->viscosity_derivative_averaging_weights[q])
587+
+ ", and newton_derivative_scaling_factor = "
588+
+ std::to_string(newton_derivative_scaling_factor) :
589+
"")));
590590

591591
for (unsigned int m=0; m<dim; ++m)
592592
for (unsigned int n=0; n<dim; ++n)
@@ -595,24 +595,28 @@ namespace aspect
595595
= effective_strain_rate[m][n];
596596

597597
active_cell_data.newton_factor_wrt_strain_rate_table(cell, q)[m][n][i]
598-
= derivatives->viscosity_derivative_wrt_strain_rate[q][m][n] *
599-
derivatives->viscosity_derivative_averaging_weights[q] *
600-
newton_derivative_scaling_factor * alpha;
598+
= newton_derivative_scaling_factor * alpha * derivatives->viscosity_derivative_wrt_strain_rate[q][m][n]
599+
* (active_cell_data.average_newton_factors ? derivatives->viscosity_derivative_averaging_weights[q] : 1.0);
601600

602601
Assert(std::isfinite(active_cell_data.strain_rate_table(cell, q)[m][n][i]),
603-
ExcMessage("active_cell_data.strain_rate_table has an element which is not finite: " + std::to_string(active_cell_data.strain_rate_table(cell, q)[m][n][i])));
602+
ExcMessage("active_cell_data.strain_rate_table has an element which is not finite: "
603+
+ std::to_string(active_cell_data.strain_rate_table(cell, q)[m][n][i])));
604604
Assert(std::isfinite(active_cell_data.newton_factor_wrt_strain_rate_table(cell, q)[m][n][i]),
605-
ExcMessage("active_cell_data.newton_factor_wrt_strain_rate_table has an element which is not finite: " + std::to_string(active_cell_data.newton_factor_wrt_strain_rate_table(cell, q)[m][n][i])));
605+
ExcMessage("active_cell_data.newton_factor_wrt_strain_rate_table has an element which is not finite: "
606+
+ std::to_string(active_cell_data.newton_factor_wrt_strain_rate_table(cell, q)[m][n][i])));
606607
}
607608

608609
if (active_cell_data.enable_prescribed_dilation)
609610
{
610611
active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]
611612
= derivatives->dilation_derivative_wrt_pressure[q] * newton_derivative_scaling_factor;
612613
Assert(std::isfinite(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]),
613-
ExcMessage("active_cell_data.dilation_derivative_wrt_pressure_table is not finite: " + std::to_string(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i]) +
614-
". Relevant variables are derivatives->dilation_derivative_wrt_pressure[q] = " + std::to_string(derivatives->dilation_derivative_wrt_pressure[q]) +
615-
" and newton_derivative_scaling_factor = " + std::to_string(newton_derivative_scaling_factor)));
614+
ExcMessage("active_cell_data.dilation_derivative_wrt_pressure_table is not finite: "
615+
+ std::to_string(active_cell_data.dilation_derivative_wrt_pressure_table(cell,q)[i])
616+
+ ". Relevant variables are derivatives->dilation_derivative_wrt_pressure[q] = "
617+
+ std::to_string(derivatives->dilation_derivative_wrt_pressure[q])
618+
+ " and newton_derivative_scaling_factor = "
619+
+ std::to_string(newton_derivative_scaling_factor)));
616620

617621
for (unsigned int m=0; m<dim; ++m)
618622
for (unsigned int n=0; n<dim; ++n)

0 commit comments

Comments
 (0)