Skip to content

Commit

Permalink
Add transmutation parameters to input file (#120)
Browse files Browse the repository at this point in the history
  • Loading branch information
@kulakovri
kulakovri authored Dec 22, 2023
1 parent 6417022 commit f2f312d
Show file tree
Hide file tree
Showing 7 changed files with 241 additions and 242 deletions.
4 changes: 4 additions & 0 deletions MDLIB/InputOutput.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1324,6 +1324,10 @@ Input ReadInputFile( char *fileName ) {
materials.axx[k] = ReadMatProps( fin, "axx", k, 1.0 ); // plastic directional factor xx
materials.azz[k] = ReadMatProps( fin, "azz", k, 1.0 ); // plastic directional factor zz
materials.ayy[k] = ReadMatProps( fin, "ayy", k, 1.0 ); // plastic directional factor yy
// temperature driven transition between phases
materials.transmutation[k] = (int)ReadMatProps( fin, "transmutation", k, 0 ); // switcher: 0 - no transition; -1 transition happens when below transition_temperature; 1 transition happens when above transition_temperature;
materials.transmutation_phase[k] = (int)ReadMatProps( fin, "transmutation_phase", k, 1 ); // phase to set
materials.transmutation_temperature[k] = ReadMatProps( fin, "transmutation_temperature", k, 1400.0 ); // temperature boundary
// Check if any flow law is active
int sum = abs(materials.cstv[k]) + abs(materials.pwlv[k]) + abs(materials.linv[k]) + abs(materials.gbsv[k]) + abs(materials.expv[k]);
if ( sum == 0 ) {
Expand Down
2 changes: 2 additions & 0 deletions MDLIB/Main_DOODZ.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -399,6 +399,8 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {

//------------------------------------------------------------------------------------------------------------------------------//

TransmutateMarkers(&particles, &input.materials, input.scaling.T);

clock_t t_omp_step = (double)omp_get_wtime();

printf(GREEN "Number of particles = %d\n" RESET, particles.Nb_part );
Expand Down
22 changes: 22 additions & 0 deletions MDLIB/ParticleRoutines.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3201,6 +3201,28 @@ void CountPartCell ( markers* particles, grid *mesh, params model, surface topo,

}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/

void TransmutateMarkers(markers *particles, mat_prop *materials, double scaling_T) {
#pragma omp parallel for shared ( particles )
int transmutated_particles = 0;
for (int k=0; k<particles->Nb_part-1; k++) {
const int phase = particles->phase[k];
const double temperature = (particles->T[k] * scaling_T) - zeroC;
const double transmutation = materials->transmutation[phase];

if ((transmutation == -1 && materials->transmutation_temperature[phase] >= temperature)
|| (transmutation == 1 && materials->transmutation_temperature[phase] <= temperature)) {
particles->phase[k] = materials->transmutation_phase[phase];
particles->dual[k] = materials->transmutation_phase[phase];
transmutated_particles++;
}
}
printf("Transmutated particles = %03d\n", transmutated_particles);
}

/*--------------------------------------------------------------------------------------------------------------------*/
/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
/*--------------------------------------------------------------------------------------------------------------------*/
2 changes: 2 additions & 0 deletions MDLIB/include/mdoodz.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,8 @@ typedef struct {
double aniso_angle[20], ani_fac_max[20], aniso_factor[20]; //ani_fac_v[20], ani_fac_e[20], ani_fac_p[20]
double axx[20], azz[20], ayy[20];
int ani_fstrain[20];
int transmutation[20], transmutation_phase[20];
double transmutation_temperature[20];
} mat_prop;

char *GetSetupFileName(int nargs, char *args[]);
Expand Down
4 changes: 3 additions & 1 deletion MDLIB/mdoodz-private.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -577,4 +577,6 @@ void SetBCs(SetBCs_ff setBCs, MdoodzInput *instance, grid *mesh);

void ValidateSetup(MdoodzSetup *setup, MdoodzInput *instance);

void SetBCs_MD6( grid*, params*, scale, markers*, mat_prop*, surface* );
void SetBCs_MD6( grid*, params*, scale, markers*, mat_prop*, surface* );

void TransmutateMarkers(markers *particles, mat_prop *materials, double scaling_T);
129 changes: 51 additions & 78 deletions SETS/NeckingReview.c
View file
Original file line number Diff line number Diff line change
@@ -1,75 +1,49 @@
#include "math.h"
#include "mdoodz.h"
#include "stdbool.h"
#include "stdlib.h"
#include "stdio.h"

int SetDualPhase(MdoodzInput *input, Coordinates coordinate, int phase) {

int dual_phase = phase;
double Lx = input->model.xmax - input->model.xmin;
double Lz = input->model.zmax - input->model.zmin;
double Ax, Az;

// Set checkerboard for phase 0
Az = sin( 48.0*2.0*M_PI*coordinate.z / Lz );
if ( Az>0.0 && dual_phase==0 ) {
dual_phase += input->model.Nb_phases;
}

// Set checkerboard for phase 1
Az = sin( 48.0*2.0*M_PI*coordinate.z / Lz );
if ( Az>0.0 && dual_phase==1 ) {
dual_phase += input->model.Nb_phases;
}

// Set checkerboard for phase 2
Az = sin( 48.0*2.0*M_PI*coordinate.z / Lz );
if ( Az>0.0 && dual_phase==2 ) {
dual_phase += input->model.Nb_phases;
}

return dual_phase;
}

#include "stdlib.h"

double SetSurfaceZCoord(MdoodzInput *instance, double x_coord) {
const double TopoLevel = -0.0e3 / instance->scaling.L;
const double basin_width = 30.0e3 / instance->scaling.L;
const double h_pert = instance->model.user3 / instance->scaling.L;
const double x = x_coord - instance->model.user4 / instance->scaling.L;
return TopoLevel + h_pert * exp( - (x*x) / (2.0*basin_width*basin_width) );
const double TopoLevel = -0.0e3 / instance->scaling.L;
return TopoLevel;
}

double SetNoise(MdoodzInput *instance, Coordinates coordinates, int phase) {
const double x = coordinates.x - instance->model.user4 / instance->scaling.L;
const double z = coordinates.z;
const double x = coordinates.x - instance->model.user4 / instance->scaling.L;
const double z = coordinates.z;
const double basin_width = 30.0e3 / instance->scaling.L;
const double noise = ((double) rand() / (double) RAND_MAX) - 0.5;
const double filter_x = exp( - (x*x)/ (2.0*basin_width*basin_width) );
const double filter_z = exp( - (z*z)/ (2.0*basin_width*basin_width*4.0) );
return noise * filter_x * filter_z;
const double noise = ((double) rand() / (double) RAND_MAX) - 0.5;
const double filter_x = exp(-(x * x) / (2.0 * basin_width * basin_width));
const double filter_z = exp(-(z * z) / (2.0 * basin_width * basin_width * 4.0));
return noise * filter_x * filter_z;
}

static Ellipse GetMicaEllipse(double centreX, double centreZ) {
return (Ellipse) {
.radiusZ = 20e3,
.radiusX = 10e3,
.centreX = centreX,
.centreZ = centreZ,
.angle = 0,
};
}

int SetPhase(MdoodzInput *instance, Coordinates coordinates) {
const double basin_width = 30.0e3 / instance->scaling.L;
const double h_pert = instance->model.user3 / instance->scaling.L;
const double lithosphereThickness = instance->model.user1 / instance->scaling.L;
const double crustThickness = instance->model.user2 / instance->scaling.L;
const double perturbationAmplitude = instance->model.user3 / instance->scaling.L;
const double x = coordinates.x - instance->model.user4 / instance->scaling.L;
const double mohoLevel = -crustThickness + 0*h_pert * exp( - (x*x) / (2.0*basin_width*basin_width) );
const bool isBelowLithosphere = coordinates.z < -lithosphereThickness;
const bool isAboveMoho = coordinates.z > mohoLevel;
const bool isLowerCrust = coordinates.z < (mohoLevel + 0.5*crustThickness);

if (isAboveMoho) {
if (isLowerCrust) {
// remove perturbation
// moho temperature should be around 500 degrees
// add weak inclusion
const double lithosphereThickness = instance->model.user1 / instance->scaling.L;
const double crustThickness = instance->model.user2 / instance->scaling.L;
const double mohoLevel = -crustThickness;
const bool isBelowLithosphere = coordinates.z < -lithosphereThickness;
const bool isAboveMoho = coordinates.z > mohoLevel;
const bool isMicaEllipse = IsEllipseCoordinates(coordinates, GetMicaEllipse(0e3, -45e3), instance->scaling.L);

if (isMicaEllipse) {
return 1;
}
else {
} else if (isAboveMoho) {
return 0;
}
} else if (isBelowLithosphere) {
return 3;
} else {
Expand All @@ -81,13 +55,13 @@ double SetTemperature(MdoodzInput *instance, Coordinates coordinates) {
const double lithosphereThickness = instance->model.user1 / instance->scaling.L;
const double surfaceTemperature = 273.15 / instance->scaling.T;
const double mantleTemperature = (1330.0 + 273.15) / instance->scaling.T;
const double Tamp = 50.0 / instance->scaling.T;
const double x = coordinates.x - instance->model.user4 / instance->scaling.L;
const double z = coordinates.z;
const double basin_width = 30.0e3 / instance->scaling.L;
const double filter_x = Tamp * exp( - (x*x)/ (2.0*basin_width*basin_width) );

const double Tamp = 50.0 / instance->scaling.T;
const double x = coordinates.x - instance->model.user4 / instance->scaling.L;
const double z = coordinates.z;
const double basin_width = 30.0e3 / instance->scaling.L;
const double filter_x = Tamp * exp(-(x * x) / (2.0 * basin_width * basin_width));

const double particleTemperature = ((mantleTemperature - surfaceTemperature) / lithosphereThickness) * (-coordinates.z) + surfaceTemperature;
if (particleTemperature > mantleTemperature) {
return mantleTemperature + filter_x;
Expand All @@ -110,8 +84,8 @@ char SetBCPType(MdoodzInput *instance, POSITION position) {
}

SetBC SetBCT(MdoodzInput *instance, POSITION position, double particleTemperature) {
SetBC bc;
double surface_temperature = zeroC / instance->scaling.T;
SetBC bc;
double surface_temperature = zeroC / instance->scaling.T;
double mantle_temperature = (1330. + zeroC) / instance->scaling.T;
if (position == S) {
bc.type = constant_temperature;
Expand All @@ -120,7 +94,7 @@ SetBC SetBCT(MdoodzInput *instance, POSITION position, double particleTemperatur
if (position == free_surface || position == N) {
bc.type = constant_temperature;
bc.value = surface_temperature;
}
}
if (position == W || position == E) {
bc.type = constant_heatflux;
bc.value = 0.;
Expand All @@ -141,25 +115,23 @@ void AddCrazyConductivity(MdoodzInput *input) {
int main(int nargs, char *args[]) {
// Input file name
char *input_file;
if ( nargs < 2 ) {
asprintf(&input_file, "NeckingReview.txt"); // Default
}
else {
if (nargs < 2) {
asprintf(&input_file, "NeckingReview.txt");// Default
} else {
printf("dodo %s\n", args[1]);
asprintf(&input_file, "%s", args[1]); // Custom
asprintf(&input_file, "%s", args[1]);// Custom
}
printf("Running MDoodz7.0 using %s\n", input_file);
srand(69); // Force random generator seed for reproducibility
srand(69);// Force random generator seed for reproducibility
MdoodzSetup setup = {
.BuildInitialTopography = &(BuildInitialTopography_ff){
.SetSurfaceZCoord = SetSurfaceZCoord,
},
.SetParticles = &(SetParticles_ff){
.SetPhase = SetPhase,
.SetTemperature = SetTemperature,
.SetGrainSize = SetGrainSize,
.SetNoise = SetNoise,
.SetDualPhase = SetDualPhase,
.SetPhase = SetPhase,
.SetTemperature = SetTemperature,
.SetGrainSize = SetGrainSize,
.SetNoise = SetNoise,

},
.SetBCs = &(SetBCs_ff){
Expand All @@ -174,3 +146,4 @@ int main(int nargs, char *args[]) {
RunMDOODZ(input_file, &setup);
free(input_file);
}

Loading

0 comments on commit f2f312d

Please sign in to comment.