Revision 0461d2e117ce88704a56dd8bcbf6bf7787991b15 authored by Eh Tan on 08 November 2007, 23:28:46 UTC, committed by Eh Tan on 08 November 2007, 23:28:46 UTC
svn+ssh://svn@geodynamics.org/cig/mc/3D/CitcomS/trunk ........ r8194 | tan2 | 2007-10-30 14:49:58 -0700 (Tue, 30 Oct 2007) | 1 line Compute d(rho)/dr/rho from rho(r) ........ r8195 | tan2 | 2007-10-30 14:50:52 -0700 (Tue, 30 Oct 2007) | 1 line Fixed a bug in dimensionalizing density. Provided the formula of geoid calculation in the comments. Rearranged the order of functions. ........ r8196 | tan2 | 2007-10-30 14:53:50 -0700 (Tue, 30 Oct 2007) | 1 line A post-processing program to project geoid coefficents onto a regular (longitude, latitude) mesh ........ r8197 | tan2 | 2007-10-30 14:54:14 -0700 (Tue, 30 Oct 2007) | 1 line Added the C program project_geoid to the makefile ........ r8199 | tan2 | 2007-10-30 15:29:44 -0700 (Tue, 30 Oct 2007) | 1 line Minor modification ........ r8201 | tan2 | 2007-11-01 16:33:30 -0700 (Thu, 01 Nov 2007) | 1 line Print dv/v=dp/p=1.0 for the 1st Uzawa iteraion ........ r8202 | tan2 | 2007-11-01 16:33:50 -0700 (Thu, 01 Nov 2007) | 1 line Fixed an error in comment ........ r8204 | tan2 | 2007-11-05 17:03:35 -0800 (Mon, 05 Nov 2007) | 1 line Scaled topo with variable gravity. Fixed an error in comment. Rearranged computation. ........ r8205 | tan2 | 2007-11-05 17:03:55 -0800 (Mon, 05 Nov 2007) | 1 line Removed functions related sph. harm in lib/Regional_obsolete.c ........ r8206 | tan2 | 2007-11-05 17:04:20 -0800 (Mon, 05 Nov 2007) | 1 line Shrank the size of sph. harm arrays ........ r8207 | tan2 | 2007-11-05 17:04:43 -0800 (Mon, 05 Nov 2007) | 1 line Init'd some variables about vtk_io, which might be accessed with uninit'd values in output_finalize() ........ r8212 | tan2 | 2007-11-06 15:17:54 -0800 (Tue, 06 Nov 2007) | 1 line Fixed a few memory errors ........ r8213 | tan2 | 2007-11-06 15:18:12 -0800 (Tue, 06 Nov 2007) | 1 line Increase vlowstep to match the default value in pyre ........ r8214 | tan2 | 2007-11-06 15:18:35 -0800 (Tue, 06 Nov 2007) | 1 line Removed unused multigrid parameters ........ r8215 | tan2 | 2007-11-06 15:18:54 -0800 (Tue, 06 Nov 2007) | 1 line Added cgrad solver convergence parameters, increased buoyancy_ratio and lower the # of steps ........ r8226 | tan2 | 2007-11-07 11:51:56 -0800 (Wed, 07 Nov 2007) | 1 line Print a warning when matrix eqn solver not converging ........ r8227 | tan2 | 2007-11-07 11:52:17 -0800 (Wed, 07 Nov 2007) | 1 line Removed comp_el from default output, since it is not required for restart anymore. ........ r8228 | tan2 | 2007-11-07 11:52:39 -0800 (Wed, 07 Nov 2007) | 1 line Decreased the # of processors. This is the only way I can reproduce single-cell convection as in the manual. ........ r8235 | tan2 | 2007-11-08 11:18:26 -0800 (Thu, 08 Nov 2007) | 1 line Dereased the timestep size to reduce artifacts in advection ........ r8236 | tan2 | 2007-11-08 11:18:52 -0800 (Thu, 08 Nov 2007) | 1 line Update NEWS ........ r8237 | tan2 | 2007-11-08 11:19:12 -0800 (Thu, 08 Nov 2007) | 1 line Update the version number ........ r8241 | tan2 | 2007-11-08 13:17:14 -0800 (Thu, 08 Nov 2007) | 1 line Updated file ChangeLog to r8240 ........ r8242 | tan2 | 2007-11-08 13:36:55 -0800 (Thu, 08 Nov 2007) | 1 line Removed binary checkpoint files from makefile, as the file size is too big for distribution. ........ r8243 | tan2 | 2007-11-08 13:38:09 -0800 (Thu, 08 Nov 2007) | 1 line Updated file ChangeLog to r8242 ........ r8244 | tan2 | 2007-11-08 14:31:21 -0800 (Thu, 08 Nov 2007) | 1 line Replaced a system call by std C library remove() and disabled another system call (backup input file). Partially fixed issue130. All remaining system calls are in lib/Output_gzdir.c. ........ r8245 | tan2 | 2007-11-08 14:41:31 -0800 (Thu, 08 Nov 2007) | 1 line Updated file ChangeLog to r8244 ........
1 parent a828fa9
Nodal_mesh.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>
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
/* 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"
void v_from_vector(E)
struct All_variables *E;
{
int m,node;
const int nno = E->lmesh.nno;
for (m=1;m<=E->sphere.caps_per_proc;m++) {
for(node=1;node<=nno;node++) {
E->sphere.cap[m].V[1][node] = E->U[m][E->id[m][node].doff[1]];
E->sphere.cap[m].V[2][node] = E->U[m][E->id[m][node].doff[2]];
E->sphere.cap[m].V[3][node] = E->U[m][E->id[m][node].doff[3]];
if (E->node[m][node] & VBX)
E->sphere.cap[m].V[1][node] = E->sphere.cap[m].VB[1][node];
if (E->node[m][node] & VBY)
E->sphere.cap[m].V[2][node] = E->sphere.cap[m].VB[2][node];
if (E->node[m][node] & VBZ)
E->sphere.cap[m].V[3][node] = E->sphere.cap[m].VB[3][node];
}
}
return;
}
void v_from_vector_pseudo_surf(E)
struct All_variables *E;
{
int m,node;
const int nno = E->lmesh.nno;
double sum_V = 0.0, sum_dV = 0.0, rel_error = 0.0, global_max_error = 0.0;
double tol_error = 1.0e-03;
for (m=1;m<=E->sphere.caps_per_proc;m++) {
for(node=1;node<=nno;node++) {
E->sphere.cap[m].Vprev[1][node] = E->sphere.cap[m].V[1][node];
E->sphere.cap[m].Vprev[2][node] = E->sphere.cap[m].V[2][node];
E->sphere.cap[m].Vprev[3][node] = E->sphere.cap[m].V[3][node];
E->sphere.cap[m].V[1][node] = E->U[m][E->id[m][node].doff[1]];
E->sphere.cap[m].V[2][node] = E->U[m][E->id[m][node].doff[2]];
E->sphere.cap[m].V[3][node] = E->U[m][E->id[m][node].doff[3]];
if (E->node[m][node] & VBX)
E->sphere.cap[m].V[1][node] = E->sphere.cap[m].VB[1][node];
if (E->node[m][node] & VBY)
E->sphere.cap[m].V[2][node] = E->sphere.cap[m].VB[2][node];
if (E->node[m][node] & VBZ)
E->sphere.cap[m].V[3][node] = E->sphere.cap[m].VB[3][node];
sum_dV += (E->sphere.cap[m].V[1][node] - E->sphere.cap[m].Vprev[1][node])*(E->sphere.cap[m].V[1][node] - E->sphere.cap[m].Vprev[1][node])
+ (E->sphere.cap[m].V[2][node] - E->sphere.cap[m].Vprev[2][node])*(E->sphere.cap[m].V[2][node] - E->sphere.cap[m].Vprev[2][node])
+ (E->sphere.cap[m].V[3][node] - E->sphere.cap[m].Vprev[3][node])*(E->sphere.cap[m].V[3][node] - E->sphere.cap[m].Vprev[3][node]);
sum_V += E->sphere.cap[m].V[1][node]*E->sphere.cap[m].V[1][node]
+ E->sphere.cap[m].V[2][node]*E->sphere.cap[m].V[2][node]
+ E->sphere.cap[m].V[3][node]*E->sphere.cap[m].V[3][node];
}
rel_error = sqrt(sum_dV)/sqrt(sum_V);
MPI_Allreduce(&rel_error,&global_max_error,1,MPI_DOUBLE,MPI_MAX,E->parallel.world);
if(global_max_error <= tol_error) E->monitor.stop_topo_loop = 1;
if(E->parallel.me==0)
fprintf(stderr,"global_max_error=%e stop_topo_loop=%d\n",global_max_error,E->monitor.stop_topo_loop);
}
return;
}
void velo_from_element(E,VV,m,el,sphere_key)
struct All_variables *E;
float VV[4][9];
int el,m,sphere_key;
{
int a, node;
double sint, cost, sinf, cosf;
const int ends=enodes[E->mesh.nsd];
const int lev=E->mesh.levmax;
if (sphere_key)
for(a=1;a<=ends;a++) {
node = E->ien[m][el].node[a];
VV[1][a] = E->sphere.cap[m].V[1][node];
VV[2][a] = E->sphere.cap[m].V[2][node];
VV[3][a] = E->sphere.cap[m].V[3][node];
}
else {
for(a=1;a<=ends;a++) {
node = E->ien[m][el].node[a];
sint = E->SinCos[lev][m][0][node];
sinf = E->SinCos[lev][m][1][node];
cost = E->SinCos[lev][m][2][node];
cosf = E->SinCos[lev][m][3][node];
VV[1][a] = E->sphere.cap[m].V[1][node]*cost*cosf
- E->sphere.cap[m].V[2][node]*sinf
+ E->sphere.cap[m].V[3][node]*sint*cosf;
VV[2][a] = E->sphere.cap[m].V[1][node]*cost*sinf
+ E->sphere.cap[m].V[2][node]*cosf
+ E->sphere.cap[m].V[3][node]*sint*sinf;
VV[3][a] = -E->sphere.cap[m].V[1][node]*sint
+ E->sphere.cap[m].V[3][node]*cost;
}
}
return;
}
void velo_from_element_d(E,VV,m,el,sphere_key)
struct All_variables *E;
double VV[4][9];
int el,m,sphere_key;
{
int a, node;
double sint, cost, sinf, cosf;
const int dims=E->mesh.nsd;
const int ends=enodes[E->mesh.nsd];
const int nno=E->lmesh.nno;
const int lev=E->mesh.levmax;
if (sphere_key)
for(a=1;a<=ends;a++) {
node = E->ien[m][el].node[a];
VV[1][a] = E->sphere.cap[m].V[1][node];
VV[2][a] = E->sphere.cap[m].V[2][node];
VV[3][a] = E->sphere.cap[m].V[3][node];
}
else {
for(a=1;a<=ends;a++) {
node = E->ien[m][el].node[a];
sint = E->SinCos[lev][m][0][node];
sinf = E->SinCos[lev][m][1][node];
cost = E->SinCos[lev][m][2][node];
cosf = E->SinCos[lev][m][3][node];
VV[1][a] = E->sphere.cap[m].V[1][node]*cost*cosf
- E->sphere.cap[m].V[2][node]*sinf
+ E->sphere.cap[m].V[3][node]*sint*cosf;
VV[2][a] = E->sphere.cap[m].V[1][node]*cost*sinf
+ E->sphere.cap[m].V[2][node]*cosf
+ E->sphere.cap[m].V[3][node]*sint*sinf;
VV[3][a] = -E->sphere.cap[m].V[1][node]*sint
+ E->sphere.cap[m].V[3][node]*cost;
}
}
return;
}
void p_to_nodes(E,P,PN,lev)
struct All_variables *E;
double **P;
float **PN;
int lev;
{ int e,element,node,j,m;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(node=1;node<=E->lmesh.NNO[lev];node++)
PN[m][node] = 0.0;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(element=1;element<=E->lmesh.NEL[lev];element++)
for(j=1;j<=enodes[E->mesh.nsd];j++) {
node = E->IEN[lev][m][element].node[j];
PN[m][node] += P[m][element] * E->TWW[lev][m][element].node[j] ;
}
(E->exchange_node_f)(E,PN,lev);
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(node=1;node<=E->lmesh.NNO[lev];node++)
PN[m][node] *= E->MASS[lev][m][node];
return;
}
void visc_from_gint_to_nodes(E,VE,VN,lev)
struct All_variables *E;
float **VE,**VN;
int lev;
{
int m,e,i,j,k,n;
const int nsd=E->mesh.nsd;
const int vpts=vpoints[nsd];
const int ends=enodes[nsd];
double temp_visc;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.NNO[lev];i++)
VN[m][i] = 0.0;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(e=1;e<=E->lmesh.NEL[lev];e++) {
temp_visc=0.0;
for(i=1;i<=vpts;i++)
temp_visc += VE[m][(e-1)*vpts + i];
temp_visc = temp_visc/vpts;
for(j=1;j<=ends;j++) {
n = E->IEN[lev][m][e].node[j];
VN[m][n] += E->TWW[lev][m][e].node[j] * temp_visc;
}
}
(E->exchange_node_f)(E,VN,lev);
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(n=1;n<=E->lmesh.NNO[lev];n++)
VN[m][n] *= E->MASS[lev][m][n];
return;
}
void visc_from_nodes_to_gint(E,VN,VE,lev)
struct All_variables *E;
float **VE,**VN;
int lev;
{
int m,e,i,j,k,n;
const int nsd=E->mesh.nsd;
const int vpts=vpoints[nsd];
const int ends=enodes[nsd];
double temp_visc;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(e=1;e<=E->lmesh.NEL[lev];e++)
for(i=1;i<=vpts;i++)
VE[m][(e-1)*vpts+i] = 0.0;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(e=1;e<=E->lmesh.NEL[lev];e++)
for(i=1;i<=vpts;i++) {
temp_visc=0.0;
for(j=1;j<=ends;j++)
temp_visc += E->N.vpt[GNVINDEX(j,i)]*VN[m][E->IEN[lev][m][e].node[j]];
VE[m][(e-1)*vpts+i] = temp_visc;
}
return;
}
void visc_from_gint_to_ele(E,VE,VN,lev)
struct All_variables *E;
float **VE,**VN;
int lev;
{
int m,e,i,j,k,n;
const int nsd=E->mesh.nsd;
const int vpts=vpoints[nsd];
const int ends=enodes[nsd];
double temp_visc;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.NEL[lev];i++)
VN[m][i] = 0.0;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(e=1;e<=E->lmesh.NEL[lev];e++) {
temp_visc=0.0;
for(i=1;i<=vpts;i++)
temp_visc += VE[m][(e-1)*vpts + i];
temp_visc = temp_visc/vpts;
VN[m][e] = temp_visc;
}
return;
}
void visc_from_ele_to_gint(E,VN,VE,lev)
struct All_variables *E;
float **VE,**VN;
int lev;
{
int m,e,i,j,k,n;
const int nsd=E->mesh.nsd;
const int vpts=vpoints[nsd];
const int ends=enodes[nsd];
double temp_visc;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(e=1;e<=E->lmesh.NEL[lev];e++)
for(i=1;i<=vpts;i++)
VE[m][(e-1)*vpts+i] = 0.0;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(e=1;e<=E->lmesh.NEL[lev];e++)
for(i=1;i<=vpts;i++) {
VE[m][(e-1)*vpts+i] = VN[m][e];
}
return;
}
Computing file changes ...
