Skip to content

corrected the error with CPO initial_volume #25

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: add_drexpp
Choose a base branch
from
Open

corrected the error with CPO initial_volume #25

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 11 additions & 3 deletions source/particle/property/crystal_preferred_orientation.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@ namespace aspect
CrystalPreferredOrientation<dim>::compute_random_rotation_matrix(Tensor<2,3> &rotation_matrix) const
{

const double dt = this -> get_time();
Copy link
Owner

@MFraters MFraters Apr 12, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

dt is usally used for a time difference, you can just call the variable time.

// This function is based on an article in Graphic Gems III, written by James Arvo, Cornell University (p 116-120).
// The original code can be found on http://www.realtimerendering.com/resources/GraphicsGems/gemsiii/rand_rotation.c
// and is licenced according to this website with the following licence:
Expand All @@ -95,7 +96,7 @@ namespace aspect

double theta = 2.0 * numbers::PI * one; // Rotation about the pole (Z)
double phi = 2.0 * numbers::PI * two; // For direction of pole deflection.
double z = 2.0* three; //For magnitude of pole deflection.
double z = dt != 0.0 ? 0.1 * three : 2.0 * three; //For magnitude of pole deflection.

// Compute a vector V used for distributing points over the sphere
// via the reflection I - V Transpose(V). This formulation of V
Expand Down Expand Up @@ -191,7 +192,7 @@ namespace aspect
for (unsigned int grain_i = 0; grain_i < n_grains ; ++grain_i)
{
// set volume fraction
if (n_grains < n_grains/10.)
if (grain_i < n_grains/10.)
volume_fractions_grains[mineral_i][grain_i] = initial_volume_fraction;
else
{
Expand Down Expand Up @@ -635,7 +636,8 @@ namespace aspect
const double homologous_T = 1770+273.15; // in Kelvin I assume. TODO: make this a function of temperature and pressure? Both are avaible as variables with those names here.
const double aggregate_recrystalization_increment = const_C*exp(const_g*temperature/homologous_T)*strain; // TODO: compute actual value with equation 5
// compute paleowatt/pizometer grainsize
const double recrystalized_grain_size = 0.5; // TODO: actually compute the grainsize with the paleowatt/pizometer

const double recrystalized_grain_size = 0.01; // TODO: actually compute the grainsize with the paleowatt/pizometer
const double half_recrystalized_grain_size = 0.5 * recrystalized_grain_size;
const double recrystalized_grain_volume = (4./3.)*numbers::PI*half_recrystalized_grain_size*half_recrystalized_grain_size*half_recrystalized_grain_size;

Expand Down Expand Up @@ -712,6 +714,12 @@ namespace aspect
const double stress_invariant = std::sqrt(std::max(-second_invariant(deviatoric_stress), 0.));

const double diffusion_creep_strain_rate = stress_invariant/(2.0*diffusion_viscosity);
//const double eff_vis =std::max(second_invariant(deviatoric_strain_rate), 0.);
//const long double pre_exponent_term =(3*1*1.92*std::pow(1/10.0,10.))/(0.1*stress_invariant*((eff_vis-diffusion_creep_strain_rate)*3));
//const long double exponent_term = -1*(400*1000)/(8.314*temperature);
//const long double recrystalized_grain_size = pre_exponent_term*std::pow(std::exp(exponent_term),0.25); // TODO: actually compute the grainsize with the paleowatt/pizometer
//const long double half_recrystalized_grain_size = 0.5 * recrystalized_grain_size;
//const long double recrystalized_grain_volume = (4./3.)*numbers::PI*half_recrystalized_grain_size*half_recrystalized_grain_size*half_recrystalized_grain_size;

return compute_derivatives_drexpp(cpo_index,
data,
Expand Down