https://github.com/geodynamics/citcoms
Revision 3a23d009ea2727d8d8c920647550f8493ae19e25 authored by Eh Tan on 17 January 2012, 02:41:15 UTC, committed by Eh Tan on 17 January 2012, 02:41:15 UTC
1 parent 15a09f1
Tip revision: 3a23d009ea2727d8d8c920647550f8493ae19e25 authored by Eh Tan on 17 January 2012, 02:41:15 UTC
When python 2.7 is used, ask users to download pythia manually. See issue622 and issue606.
When python 2.7 is used, ask users to download pythia manually. See issue622 and issue606.
Tip revision: 3a23d00
Sphere_harmonics.c
/* Functions relating to the building and use of mesh locations ... */
#include <math.h>
#include <sys/types.h>
#include "element_definitions.h"
#include "global_defs.h"
#include <stdlib.h>
static void compute_sphereh_table(struct All_variables *);
/* ======================================================================
====================================================================== */
void set_sphere_harmonics(E)
struct All_variables *E;
{
int m,node,ll,mm,i,j;
i=0;
for (ll=0;ll<=E->output.llmax;ll++)
for (mm=0;mm<=ll;mm++) {
E->sphere.hindex[ll][mm] = i;
i++;
}
E->sphere.hindice = i;
/* spherical harmonic coeff (0=cos, 1=sin)
for surface topo, cmb topo and geoid */
for (i=0;i<=1;i++) {
E->sphere.harm_geoid[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
E->sphere.harm_geoid_from_bncy[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
E->sphere.harm_geoid_from_bncy_botm[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
E->sphere.harm_geoid_from_tpgt[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
E->sphere.harm_geoid_from_tpgb[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
E->sphere.harm_tpgt[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
E->sphere.harm_tpgb[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
}
compute_sphereh_table(E);
return;
}
/* =====================================
Generalized Legendre polynomials
=====================================*/
double modified_plgndr_a(int l, int m, double t)
{
int i,ll;
double x,fact1,fact2,fact,pll,pmm,pmmp1,somx2,plgndr;
const double three=3.0;
const double two=2.0;
const double one=1.0;
x = cos(t);
pmm=one;
if(m>0) {
somx2=sqrt((one-x)*(one+x));
fact1= three;
fact2= two;
for (i=1;i<=m;i++) {
fact=sqrt(fact1/fact2);
pmm = -pmm*fact*somx2;
fact1+= two;
fact2+= two;
}
}
if (l==m)
plgndr = pmm;
else {
pmmp1 = x*sqrt(two*m+three)*pmm;
if(l==m+1)
plgndr = pmmp1;
else {
for (ll=m+2;ll<=l;ll++) {
fact1= sqrt((4.0*ll*ll-one)*(double)(ll-m)/(double)(ll+m));
fact2= sqrt((2.0*ll+one)*(ll-m)*(ll+m-one)*(ll-m-one)
/(double)((two*ll-three)*(ll+m)));
pll = ( x*fact1*pmmp1-fact2*pmm)/(ll-m);
pmm = pmmp1;
pmmp1 = pll;
}
plgndr = pll;
}
}
plgndr /= sqrt(4.0*M_PI);
if (m!=0) plgndr *= sqrt(two);
return plgndr;
}
/* =========================================================
expand the field TG into spherical harmonics
========================================================= */
void sphere_expansion(E,TG,sphc,sphs)
struct All_variables *E;
float **TG,*sphc,*sphs;
{
int el,nint,d,p,i,m,j,es,mm,ll,rand();
void sum_across_surf_sph1();
for (i=0;i<E->sphere.hindice;i++) {
sphc[i] = 0.0;
sphs[i] = 0.0;
}
for (m=1;m<=E->sphere.caps_per_proc;m++)
for (es=1;es<=E->lmesh.snel;es++) {
for (ll=0;ll<=E->output.llmax;ll++)
for (mm=0; mm<=ll; mm++) {
p = E->sphere.hindex[ll][mm];
for(nint=1;nint<=onedvpoints[E->mesh.nsd];nint++) {
for(d=1;d<=onedvpoints[E->mesh.nsd];d++) {
j = E->sien[m][es].node[d];
sphc[p] += TG[m][E->sien[m][es].node[d]]
* E->sphere.tablesplm[m][j][p]
* E->sphere.tablescosf[m][j][mm]
* E->M.vpt[GMVINDEX(d,nint)]
* E->surf_det[m][nint][es];
sphs[p] += TG[m][E->sien[m][es].node[d]]
* E->sphere.tablesplm[m][j][p]
* E->sphere.tablessinf[m][j][mm]
* E->M.vpt[GMVINDEX(d,nint)]
* E->surf_det[m][nint][es];
}
}
} /* end for ll and mm */
}
sum_across_surf_sph1(E,sphc,sphs);
return;
}
void debug_sphere_expansion(struct All_variables *E)
{
/* expand temperature field (which should be a sph. harm. load)
* and output the expansion coeff. to stderr
*/
int m, i, j, k, p, node;
int ll, mm;
float *TT[NCS], *sph_harm[2];
for(m=1;m<=E->sphere.caps_per_proc;m++)
TT[m] = (float *) malloc ((E->lmesh.nsf+1)*sizeof(float));
/* sin coeff */
sph_harm[0] = (float*)malloc(E->sphere.hindice*sizeof(float));
/* cos coeff */
sph_harm[1] = (float*)malloc(E->sphere.hindice*sizeof(float));
for(k=1;k<=E->lmesh.noz;k++) {
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.noy;i++)
for(j=1;j<=E->lmesh.nox;j++) {
node= k + (j-1)*E->lmesh.noz + (i-1)*E->lmesh.nox*E->lmesh.noz;
p = j + (i-1)*E->lmesh.nox;
TT[m][p] = E->T[m][node];
}
/* expand TT into spherical harmonics */
sphere_expansion(E, TT, sph_harm[0], sph_harm[1]);
/* only the first nprocz CPU needs output */
if(E->parallel.me < E->parallel.nprocz) {
for (ll=0;ll<=E->output.llmax;ll++)
for (mm=0; mm<=ll; mm++) {
p = E->sphere.hindex[ll][mm];
fprintf(stderr, "T expanded layer=%d ll=%d mm=%d -- %12g %12g\n",
k+E->lmesh.nzs-1, ll, mm,
sph_harm[0][p], sph_harm[1][p]);
}
}
}
return;
}
/* ==================================================*/
/* ==================================================*/
static void compute_sphereh_table(E)
struct All_variables *E;
{
double modified_plgndr_a();
int m,node,ll,mm,i,j,p;
double t,f,mmf;
for(m=1;m<=E->sphere.caps_per_proc;m++) {
E->sphere.tablesplm[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
E->sphere.tablescosf[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
E->sphere.tablessinf[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
for (i=1;i<=E->lmesh.nsf;i++) {
E->sphere.tablesplm[m][i]= (double *)malloc((E->sphere.hindice)*sizeof(double));
E->sphere.tablescosf[m][i]= (double *)malloc((E->output.llmax+1)*sizeof(double));
E->sphere.tablessinf[m][i]= (double *)malloc((E->output.llmax+1)*sizeof(double));
}
}
for(m=1;m<=E->sphere.caps_per_proc;m++) {
for (j=1;j<=E->lmesh.nsf;j++) {
node = j*E->lmesh.noz;
f=E->sx[m][2][node];
t=E->sx[m][1][node];
for (mm=0;mm<=E->output.llmax;mm++) {
mmf = (double)(mm)*f;
E->sphere.tablescosf[m][j][mm] = cos( mmf );
E->sphere.tablessinf[m][j][mm] = sin( mmf );
}
for (ll=0;ll<=E->output.llmax;ll++)
for (mm=0;mm<=ll;mm++) {
p = E->sphere.hindex[ll][mm];
E->sphere.tablesplm[m][j][p] = modified_plgndr_a(ll,mm,t) ;
}
}
}
return;
}
Computing file changes ...