https://github.com/geodynamics/citcoms
Revision df821c06952128812763ec5b7a944cd178ca5ed3 authored by Thorsten Becker on 08 December 2010, 18:43:58 UTC, committed by Thorsten Becker on 08 December 2010, 18:43:58 UTC
GMT/ggrd compile before).
1 parent cf7495e
Tip revision: df821c06952128812763ec5b7a944cd178ca5ed3 authored by Thorsten Becker on 08 December 2010, 18:43:58 UTC
Added definition of TRUE (1) and FALSE (0) in case undefined (was defined for
Added definition of TRUE (1) and FALSE (0) in case undefined (was defined for
Tip revision: df821c0
Phase_change.c
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
* Clint Conrad, Michael Gurnis, and Eun-seo Choi.
* Copyright (C) 1994-2005, California Institute of Technology.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include <stdio.h>
#include <math.h>
#include <sys/types.h>
#include "global_defs.h"
#include "parsing.h"
#include "phase_change.h"
static void phase_change_apply(struct All_variables *E, double **buoy,
float **B, float **B_b,
float Ra, float clapeyron,
float depth, float transT, float inv_width);
static void calc_phase_change(struct All_variables *E,
float **B, float **B_b,
float Ra, float clapeyron,
float depth, float transT, float inv_width);
static void debug_phase_change(struct All_variables *E, float **B);
void phase_change_allocate(struct All_variables *E)
{
int j;
int nno = E->lmesh.nno;
int nsf = E->lmesh.nsf;
for (j=1;j<=E->sphere.caps_per_proc;j++) {
E->Fas410[j] = (float *) malloc((nno+1)*sizeof(float));
E->Fas410_b[j] = (float *) malloc((nsf+1)*sizeof(float));
E->Fas670[j] = (float *) malloc((nno+1)*sizeof(float));
E->Fas670_b[j] = (float *) malloc((nsf+1)*sizeof(float));
E->Fascmb[j] = (float *) malloc((nno+1)*sizeof(float));
E->Fascmb_b[j] = (float *) malloc((nsf+1)*sizeof(float));
}
return;
}
void phase_change_input(struct All_variables *E)
{
int m = E->parallel.me;
float width;
/* for phase change 410km */
input_float("Ra_410",&(E->control.Ra_410),"0.0",m);
input_float("clapeyron410",&(E->control.clapeyron410),"0.0",m);
input_float("transT410",&(E->control.transT410),"0.0",m);
input_float("width410",&width,"0.0",m);
if (width!=0.0)
E->control.inv_width410 = 1.0/width;
/* for phase change 670km */
input_float("Ra_670",&(E->control.Ra_670),"0.0",m);
input_float("clapeyron670",&(E->control.clapeyron670),"0.0",m);
input_float("transT670",&(E->control.transT670),"0.0",m);
input_float("width670",&width,"0.0",m);
if (width!=0.0)
E->control.inv_width670 = 1.0/width;
/* for phase change CMB */
input_float("Ra_cmb",&(E->control.Ra_cmb),"0.0",m);
input_float("clapeyroncmb",&(E->control.clapeyroncmb),"0.0",m);
input_float("transTcmb",&(E->control.transTcmb),"0.0",m);
input_float("widthcmb",&width,"0.0",m);
if (width!=0.0)
E->control.inv_widthcmb = 1.0/width;
return;
}
void phase_change_apply_410(struct All_variables *E, double **buoy)
{
if (E->control.Ra_410 != 0.0)
phase_change_apply(E, buoy, E->Fas410, E->Fas410_b, E->control.Ra_410,
E->control.clapeyron410, E->viscosity.z410,
E->control.transT410, E->control.inv_width410);
return;
}
void phase_change_apply_670(struct All_variables *E, double **buoy)
{
if (E->control.Ra_670 != 0.0)
phase_change_apply(E, buoy, E->Fas670, E->Fas670_b, E->control.Ra_670,
E->control.clapeyron670, E->viscosity.zlm,
E->control.transT670, E->control.inv_width670);
return;
}
void phase_change_apply_cmb(struct All_variables *E, double **buoy)
{
if (E->control.Ra_cmb != 0.0)
phase_change_apply(E, buoy, E->Fascmb, E->Fascmb_b, E->control.Ra_cmb,
E->control.clapeyroncmb, E->viscosity.zcmb,
E->control.transTcmb, E->control.inv_widthcmb);
return;
}
static void phase_change_apply(struct All_variables *E, double **buoy,
float **B, float **B_b,
float Ra, float clapeyron,
float depth, float transT, float inv_width)
{
int m, i;
calc_phase_change(E, B, B_b, Ra, clapeyron, depth, transT, inv_width);
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.nno;i++)
buoy[m][i] -= Ra * B[m][i];
if (E->control.verbose) {
fprintf(E->fp_out, "Ra=%f, clapeyron=%f, depth=%f, transT=%f, inv_width=%f\n",
Ra, clapeyron, depth, transT, inv_width);
debug_phase_change(E,B);
fflush(E->fp_out);
}
return;
}
static void calc_phase_change(struct All_variables *E,
float **B, float **B_b,
float Ra, float clapeyron,
float depth, float transT, float inv_width)
{
int i,j,k,n,ns,m,nz;
float e_pressure,pt5,one,dz;
pt5 = 0.5;
one = 1.0;
for(m=1;m<=E->sphere.caps_per_proc;m++) {
/* compute phase function B, the concentration of the high pressure
* phase. B is between 0 and 1. */
for(i=1;i<=E->lmesh.nno;i++) {
nz = ((i-1) % E->lmesh.noz) + 1;
dz = (E->sphere.ro-E->sx[m][3][i]) - depth;
/*XXX: dz*rho[nz]*g[nz] is only a approximation for the reduced
* pressure, a more accurate formula is:
* integral(rho(z)*g(z)*dz) from depth_ph to current depth */
e_pressure = dz * E->refstate.rho[nz] * E->refstate.gravity[nz]
- clapeyron * (E->T[m][i] - transT);
B[m][i] = pt5 * (one + tanh(inv_width * e_pressure));
}
/* compute the phase boundary, defined as the depth where B==0.5 */
ns = 0;
for (k=1;k<=E->lmesh.noy;k++)
for (j=1;j<=E->lmesh.nox;j++) {
ns++;
B_b[m][ns]=0.0;
for (i=1;i<E->lmesh.noz;i++) {
n = (k-1)*E->lmesh.noz*E->lmesh.nox + (j-1)*E->lmesh.noz + i;
if (B[m][n]>=pt5 && B[m][n+1]<=pt5)
B_b[m][ns]=(E->sx[m][3][n+1]-E->sx[m][3][n])*(pt5-B[m][n])/(B[m][n+1]-B[m][n])+E->sx[m][3][n];
}
}
}
return;
}
static void debug_phase_change(struct All_variables *E, float **B)
{
int m, j;
fprintf(E->fp_out,"output_phase_change_buoyancy\n");
for(m=1;m<=E->sphere.caps_per_proc;m++) {
fprintf(E->fp_out,"for cap %d\n",E->sphere.capid[m]);
for (j=1;j<=E->lmesh.nno;j++)
fprintf(E->fp_out,"Z = %.6e T = %.6e B[%06d] = %.6e \n",E->sx[m][3][j],E->T[m][j],j,B[m][j]);
}
fflush(E->fp_out);
return;
}
Computing file changes ...