Skip to content

Commit

Permalink
commenting some some stuff out and changing a bit...
Browse files Browse the repository at this point in the history
  • Loading branch information
@MFraters
MFraters committed Feb 19, 2018
1 parent 949d99c commit c956037
Show file tree
Hide file tree
Showing 4 changed files with 12 additions and 11 deletions.
2 changes: 1 addition & 1 deletion build_9.0/nonlinear_channel_flow_tractions_Newton_Stokes.prm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
set Dimension = 2
set CFL number = 1.0
set Maximum time step = 1
set End time = 1
set End time = 0
set Start time = 0
set Adiabatic surface temperature = 0
set Surface pressure = 0
Expand Down
11 changes: 6 additions & 5 deletions source/material_model/drucker_prager_compositions.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,7 @@ namespace aspect
//const double cohesion = 1e8;
//const double phi = 30 * numbers::PI/180;
//strictly speaking the derivative is this: 0.5 * ((1/stress_exponent)-1) * std::pow(2,2) * out.viscosities[i] * (1/(edot_ii*edot_ii)) * deviator(in.strain_rate[i])
composition_viscosities_derivatives[c] = -0.5*(composition_viscosities[c]/(strain_rate_norm_square))*deviator_strain_rate;//*(cohesion*std::cos(phi)+in.pressure[i]*std::sin(phi));//alpha * 1 * composition_viscosities[c] * (1/(edot_ii * edot_ii + eref * eref)) * in.strain_rate[i];
composition_viscosities_derivatives[c] = -(composition_viscosities[c]/(strain_rate_norm_square))*deviator_strain_rate;//*(cohesion*std::cos(phi)+in.pressure[i]*std::sin(phi));//alpha * 1 * composition_viscosities[c] * (1/(edot_ii * edot_ii + eref * eref)) * in.strain_rate[i];
/**
* EQ 2: Raids euqation of powerlaw
*/
Expand Down Expand Up @@ -317,7 +317,7 @@ namespace aspect
//const double deviator_strain_rate_norm = std::sqrt(deviator_strain_rate*deviator_strain_rate);

//std::cout << "strain_rate = " << strain_rate << ", strain-rated diff" << std::max(std::fabs(strain_rate[strain_rate_indices]), min_strain_rate[c]*min_strain_rate[c]);
const SymmetricTensor<2,dim> strain_rate_difference_plus = deviator_strain_rate + 0.5 * std::max(edot_ii, min_strain_rate[c]*min_strain_rate[c]) * finite_difference_accuracy
const SymmetricTensor<2,dim> strain_rate_difference_plus = deviator_strain_rate + std::max(edot_ii, min_strain_rate[c]*min_strain_rate[c]) * finite_difference_accuracy
* Utilities::nth_basis_for_symmetric_tensors<dim>(component);
const double second_invariant_strain_rate_difference_plus = compute_second_invariant(strain_rate_difference_plus, min_strain_rate[c]);
const double eta_component_plus = compute_viscosity(second_invariant_strain_rate_difference_plus,pressure,c,prefactor[c],alpha,eref,min_visc[c],max_visc[c]);
Expand All @@ -334,16 +334,17 @@ namespace aspect
const double eta_component_min = compute_viscosity(second_invariant_strain_rate_difference_min,pressure,c,prefactor[c],alpha,eref,min_visc[c],max_visc[c]);
//std::cout << ", second_invariant_strain_rate_difference = " << second_invariant_strain_rate_difference << ", edot_ii = " << edot_ii;

std::cout << std::setprecision(12) << "eta_component_plus = " << eta_component_plus << ", eta_component = " << eta_component << ", composition_viscosities[c] = " << composition_viscosities[c] << ", eta_component_min " << eta_component_min
<< ", eta_component_plus-composition_viscosities[c] = " << eta_component_plus-composition_viscosities[c] <<", eta_component_plus-eta_component_min = " << eta_component_plus-eta_component_min << ", eta_component-composition_viscosities[c]" << eta_component-composition_viscosities[c] << std::endl;
//std::cout << std::setprecision(12) << "eta_component_plus = " << eta_component_plus << ", eta_component = " << eta_component << ", composition_viscosities[c] = " << composition_viscosities[c] << ", eta_component_min " << eta_component_min
// << ", eta_component_plus-composition_viscosities[c] = " << eta_component_plus-composition_viscosities[c] <<", eta_component_plus-eta_component_min = " << eta_component_plus-eta_component_min << ", eta_component-composition_viscosities[c]" << eta_component-composition_viscosities[c] << std::endl;
// compute the difference between the viscosity with and without the strain-rate difference.
double viscosity_derivative = eta_component_plus - eta_component_min;
double viscosity_derivative = eta_component_plus - eta_component;
if (viscosity_derivative != 0)
{
// when the difference is non-zero, divide by the difference.
viscosity_derivative /= std::max(edot_ii, min_strain_rate[c]*min_strain_rate[c]) * finite_difference_accuracy;
}
composition_viscosities_derivatives[c][strain_rate_indices] = viscosity_derivative;
//std::cout << "analytical = " << composition_viscosities_derivatives[c][strain_rate_indices] << ", fd = " << viscosity_derivative << std::endl;
//std::cout << ", viscosity_derivative = " << viscosity_derivative << std::endl;
//}
}
Expand Down
4 changes: 2 additions & 2 deletions source/simulator/assemblers/newton_stokes.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -313,10 +313,10 @@ namespace aspect
:
1;

if((velocity_block_stabilization & Newton::Parameters::Stabilization::PD) != Newton::Parameters::Stabilization::none)
/*if((velocity_block_stabilization & Newton::Parameters::Stabilization::PD) != Newton::Parameters::Stabilization::none)
std::cout << "alpha is used resuting in: " << alpha << std::endl;
else
std::cout << "alpha is not really computed: " << alpha << std::endl;
std::cout << "alpha is not really computed: " << alpha << std::endl;*/

// symmetrize when the stabilization is symmetric or SPD
if ((velocity_block_stabilization & Newton::Parameters::Stabilization::symmetric) != Newton::Parameters::Stabilization::none)
Expand Down
6 changes: 3 additions & 3 deletions source/utilities.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2415,16 +2415,16 @@ namespace aspect
if ((strain_rate.norm() == 0) || (dviscosities_dstrain_rate.norm() == 0))
return 1;

const double norm_a_b = std::sqrt((deviator(strain_rate)*deviator(strain_rate))*(dviscosities_dstrain_rate*dviscosities_dstrain_rate));
const double norm_a_b = std::sqrt((strain_rate*strain_rate)*(dviscosities_dstrain_rate*dviscosities_dstrain_rate));;//std::sqrt((deviator(strain_rate)*deviator(strain_rate))*(dviscosities_dstrain_rate*dviscosities_dstrain_rate));
const double contract_b_a = (dviscosities_dstrain_rate*strain_rate);
const double one_minus_part = 1 - (contract_b_a / norm_a_b);
const double denom = one_minus_part * one_minus_part * norm_a_b;

std::cout << std::endl;
/*std::cout << std::endl;
std::cout << "strain_rate = " << strain_rate << ", deviator(strain_rate) = " << deviator(strain_rate) << ", dviscosities_dstrain_rate = " << dviscosities_dstrain_rate << ", norm_a_b = " << norm_a_b << ", contract_b_a = " << contract_b_a << std::endl;
std::cout << "eta = " << eta << ", 2 * eta = " << 2*eta << ", denom = " << denom << std::endl;
if(denom > 0)
std::cout << ", ((2.0 * eta) / denom) = " << ((2.0 * eta) / denom) << std::endl;
std::cout << ", ((2.0 * eta) / denom) = " << ((2.0 * eta) / denom) << std::endl;*/
// the case denom == 0 (smallest eigenvalue is zero), should return one,
// and it does here, because C_safety * 2.0 * eta is always larger then zero.
if (denom <= SPD_safety_factor * 2.0 * eta)
Expand Down

0 comments on commit c956037

Please sign in to comment.