We are hiring ! See our job offers.
##### https://github.com/geodynamics/citcoms
Tip revision: 80d469e
Solver_multigrid.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 "element_definitions.h"
#include "global_defs.h"
#include <math.h>

/*
static double weightIJ[5][5] =
{       { 0.0,    0.0,    0.0,    0.0,    0.0 },
{ 0.0, 0.5625, 0.1875, 0.0625, 0.1875 },
{ 0.0, 0.1875, 0.5625, 0.1875, 0.0625 },
{ 0.0, 0.0625, 0.1875, 0.5625, 0.1875 },
{ 0.0, 0.1875, 0.0625, 0.1875, 0.5625 } };
*/

void set_mg_defaults(E)
struct All_variables *E;
{ void assemble_forces_iterative();
void solve_constrained_flow_iterative();
void mg_allocate_vars();

E->solver_allocate_vars = mg_allocate_vars;
E->build_forcing_term = assemble_forces_iterative;
E->solve_stokes_problem = solve_constrained_flow_iterative;

return;
}

void mg_allocate_vars(E)
struct All_variables *E;
{
return;

}

/* =====================================================
Function to inject data from fine to coarse grid (i.e.
just dropping values at shared grid points.
===================================================== */

void inject_scalar(E,start_lev,AU,AD)
struct All_variables *E;
int start_lev;
float **AU,**AD;  /* data on upper/lower mesh  */

{
int i,m,el,node_coarse,node_fine,sl_minus,eqn,eqn_coarse;
void gather_to_1layer_node ();

const int dims = E->mesh.nsd;
const int ends = enodes[dims];

if(start_lev == E->mesh.levmin)   {
fprintf(E->fp,"Warning, attempting to project below lowest level\n");
return;
}

sl_minus = start_lev-1;

for (m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=E->lmesh.NEL[sl_minus];el++)
for(i=1;i<=ends;i++)       {
node_coarse = E->IEN[sl_minus][m][el].node[i];
node_fine=E->IEN[start_lev][m][E->EL[sl_minus][m][el].sub[i]].node[i];
AD[m][node_coarse] = AU[m][node_fine];
}

return;
}

void inject_vector(E,start_lev,AU,AD)
struct All_variables *E;
int start_lev;
double **AU,**AD;  /* data on upper/lower mesh  */

{
int i,j,m,el,node_coarse,node_fine,sl_minus,eqn_fine,eqn_coarse;

const int dims = E->mesh.nsd;
const int ends = enodes[dims];

void gather_to_1layer_id ();

if(start_lev == E->mesh.levmin)   {
fprintf(E->fp,"Warning, attempting to project below lowest level\n");
return;
}

sl_minus = start_lev-1;

for (m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=E->lmesh.NEL[sl_minus];el++)
for(i=1;i<=ends;i++)       {
node_coarse = E->IEN[sl_minus][m][el].node[i];
node_fine=E->IEN[start_lev][m][E->EL[sl_minus][m][el].sub[i]].node[i];
for (j=1;j<=dims;j++)    {
eqn_fine   = E->ID[start_lev][m][node_fine].doff[j];
eqn_coarse = E->ID[sl_minus][m][node_coarse].doff[j];
AD[m][eqn_coarse] = AU[m][eqn_fine];
}
}

return;
}

/* =====================================================
Function to inject data from coarse to fine grid (i.e.
just dropping values at shared grid points.
===================================================== */

void un_inject_vector(E,start_lev,AD,AU)

struct All_variables *E;
int start_lev;
double **AU,**AD;  /* data on upper/lower mesh  */
{
int i,m;
int el,node,node_plus;
int eqn1,eqn_plus1;
int eqn2,eqn_plus2;
int eqn3,eqn_plus3;

const int dims = E->mesh.nsd;
const int ends = enodes[dims];
const int sl_plus = start_lev+1;
const int neq = E->lmesh.NEQ[sl_plus];
const int nels = E->lmesh.NEL[start_lev];

assert(start_lev != E->mesh.levmax  /* un_injection */);

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=neq;i++)
AU[m][i]=0.0;

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=nels;el++)
for(i=1;i<=ENODES3D;i++)  {
node = E->IEN[start_lev][m][el].node[i];
node_plus=E->IEN[sl_plus][m][E->EL[start_lev][m][el].sub[i]].node[i];

eqn1 = E->ID[start_lev][m][node].doff[1];
eqn2 = E->ID[start_lev][m][node].doff[2];
eqn3 = E->ID[start_lev][m][node].doff[3];
eqn_plus1 = E->ID[sl_plus][m][node_plus].doff[1];
eqn_plus2 = E->ID[sl_plus][m][node_plus].doff[2];
eqn_plus3 = E->ID[sl_plus][m][node_plus].doff[3];
AU[m][eqn_plus1] = AD[m][eqn1];
AU[m][eqn_plus2] = AD[m][eqn2];
AU[m][eqn_plus3] = AD[m][eqn3];
}

return;
}

/* =======================================================================================
Interpolation from coarse grid to fine. See the appology attached to project() if you get
stressed out by node based assumptions. If it makes you feel any better, I don't like
it much either.
======================================================================================= */

void interp_vector(E,start_lev,AD,AU)

struct All_variables *E;
int start_lev;
double **AD,**AU;  /* data on upper/lower mesh  */
{
void un_inject_vector();
void fill_in_gaps();
void from_rtf_to_xyz();
void from_xyz_to_rtf();
void scatter_to_nlayer_id();

int i,j,k,m;
float x1,x2;
float n1,n2;
int noxz,node0,node1,node2;
int eqn0,eqn1,eqn2;

const int level = start_lev + 1;
const int dims =E->mesh.nsd;
const int ends= enodes[dims];

const int nox = E->lmesh.NOX[level];
const int noz = E->lmesh.NOZ[level];
const int noy = E->lmesh.NOY[level];
const int high_eqn = E->lmesh.NEQ[level];

if (start_lev==E->mesh.levmax) return;

from_rtf_to_xyz(E,start_lev,AD,AU);    /* transform in xyz coordinates */
un_inject_vector(E,start_lev,AU,E->temp); /*  information from lower level */
fill_in_gaps(E,E->temp,level);
from_xyz_to_rtf(E,level,E->temp,AU);      /* get back to rtf coordinates */

return;

}

/*  ==============================================
function to project viscosity down to all the
levels in the problem. (no gaps for vbcs)
==============================================  */

void project_viscosity(E)
struct All_variables *E;

{
int lv,i,j,k,m,sl_minus;
void inject_scalar();
void project_scalar();
void project_scalar_e();
void inject_scalar_e();
void visc_from_gint_to_nodes();
void visc_from_nodes_to_gint();
void visc_from_gint_to_ele();
void visc_from_ele_to_gint();

const int nsd=E->mesh.nsd;
const int vpts=vpoints[nsd];

float *viscU[NCS],*viscD[NCS];

lv = E->mesh.levmax;

for(m=1;m<=E->sphere.caps_per_proc;m++)  {
viscU[m]=(float *)malloc((1+E->lmesh.NNO[lv])*sizeof(float));
viscD[m]=(float *)malloc((1+vpts*E->lmesh.NEL[lv-1])*sizeof(float));
}

for(lv=E->mesh.levmax;lv>E->mesh.levmin;lv--)     {

sl_minus = lv -1;

if (E->viscosity.smooth_cycles==1)  {
/*      visc_from_gint_to_nodes(E,E->EVI[lv],viscU,lv);
project_scalar(E,lv,viscU,viscD);
visc_from_nodes_to_gint(E,viscD,E->EVI[sl_minus],sl_minus);  */

project_scalar(E,lv,E->VI[lv],E->VI[sl_minus]);
visc_from_nodes_to_gint(E,E->VI[sl_minus],E->EVI[sl_minus],sl_minus);
}
else if (E->viscosity.smooth_cycles==2)   {
visc_from_gint_to_ele(E,E->EVI[lv],viscU,lv);
inject_scalar_e(E,lv,viscU,E->EVI[sl_minus]);
}
else if (E->viscosity.smooth_cycles==3)   {
visc_from_gint_to_ele(E,E->EVI[lv],viscU,lv);
project_scalar_e(E,lv,viscU,viscD);
visc_from_ele_to_gint(E,viscD,E->EVI[sl_minus],sl_minus);
}
else if (E->viscosity.smooth_cycles==0)  {
/*      visc_from_gint_to_nodes(E,E->EVI[lv],viscU,lv);
inject_scalar(E,lv,viscU,viscD);
visc_from_nodes_to_gint(E,viscD,E->EVI[sl_minus],sl_minus); */

inject_scalar(E,lv,E->VI[lv],E->VI[sl_minus]);
visc_from_nodes_to_gint(E,E->VI[sl_minus],E->EVI[sl_minus],sl_minus);
}

/*        for(m=1;m<=E->sphere.caps_per_proc;m++) {
for (i=1;i<=E->lmesh.NEL[lv-1];i++)
fprintf (E->fp_out,"%d %g\n",i,viscD[m][i]);
for (i=1;i<=E->lmesh.NEL[lv];i++)
fprintf (E->fp_out,"%d %g\n",i,viscU[m][i]);
}
*/
}

for(m=1;m<=E->sphere.caps_per_proc;m++)  {
free((void *)viscU[m]);
free((void *)viscD[m]);
}

return;
}

/* ==================================================== */
void inject_scalar_e(E,start_lev,AU,AD)

struct All_variables *E;
int start_lev;
float **AU,**AD;  /* data on upper/lower mesh  */
{
int i,j,m;
int el,node,e;
float average,w;
void gather_to_1layer_ele ();

const int sl_minus = start_lev-1;
const int nels_minus=E->lmesh.NEL[start_lev-1];
const int dims=E->mesh.nsd;
const int ends=enodes[E->mesh.nsd];
const int vpts=vpoints[E->mesh.nsd];
const int n_minus=nels_minus*vpts;

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=n_minus;i++)
AD[m][i] = 0.0;

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=nels_minus;el++)
for(i=1;i<=ENODES3D;i++)                {
e = E->EL[sl_minus][m][el].sub[i];
AD[m][(el-1)*vpts+i] = AU[m][e];
}

return;
}

/* ==================================================== */
void project_scalar_e(E,start_lev,AU,AD)

struct All_variables *E;
int start_lev;
float **AU,**AD;  /* data on upper/lower mesh  */
{
int i,j,m;
int el,node,e;
float average,w;
void gather_to_1layer_ele ();

const int sl_minus = start_lev-1;
const int nels_minus=E->lmesh.NEL[start_lev-1];
const int  dims=E->mesh.nsd;
const int ends=enodes[E->mesh.nsd];
const double weight=(double) 1.0/ends;
const int vpts=vpoints[E->mesh.nsd];
const int n_minus=nels_minus*vpts;

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=n_minus;i++)
AD[m][i] = 0.0;

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=nels_minus;el++)    {
average=0.0;
for(i=1;i<=ENODES3D;i++) {
e = E->EL[sl_minus][m][el].sub[i];
average += AU[m][e];
}

AD[m][el] = average*weight;
}
return;
}

/* ==================================================== */
void project_scalar(E,start_lev,AU,AD)

struct All_variables *E;
int start_lev;
float **AU,**AD;  /* data on upper/lower mesh  */
{
int i,j,m;
int el,node,node1;
float average,w;
void gather_to_1layer_node ();

const int sl_minus = start_lev-1;
const int nno_minus=E->lmesh.NNO[start_lev-1];
const int nels_minus=E->lmesh.NEL[start_lev-1];
const int  dims=E->mesh.nsd;
const int ends=enodes[E->mesh.nsd];
const double weight=(double) 1.0/ends;

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nno_minus;i++)
AD[m][i] = 0.0;

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=nels_minus;el++)
for(i=1;i<=ENODES3D;i++) {
average=0.0;
node1 = E->EL[sl_minus][m][el].sub[i];
for(j=1;j<=ENODES3D;j++)                     {
node=E->IEN[start_lev][m][node1].node[j];
average += AU[m][node];
}

w=weight*average;

node= E->IEN[sl_minus][m][el].node[i];

AD[m][node] += w * E->TWW[sl_minus][m][el].node[i];
}

(E->exchange_node_f)(E,AD,sl_minus);

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nno_minus;i++)  {
AD[m][i] *= E->MASS[sl_minus][m][i];
}

return;
}

/* this is prefered scheme with averages */

void project_vector(E,start_lev,AU,AD,ic)

struct All_variables *E;
int start_lev,ic;
double **AU,**AD;  /* data on upper/lower mesh  */
{
int i,j,m;
int el,node1,node,e1;
int eqn1,eqn_minus1;
int eqn2,eqn_minus2;
int eqn3,eqn_minus3;
double average1,average2,average3,w,weight;
float CPU_time(),time;

void from_rtf_to_xyz();
void from_xyz_to_rtf();
void gather_to_1layer_id ();

const int sl_minus = start_lev-1;
const int neq_minus=E->lmesh.NEQ[start_lev-1];
const int nno_minus=E->lmesh.NNO[start_lev-1];
const int nels_minus=E->lmesh.NEL[start_lev-1];
const int  dims=E->mesh.nsd;
const int ends=enodes[E->mesh.nsd];

if (ic==1)
weight = 1.0;
else
weight=(double) 1.0/ends;

if (start_lev==E->mesh.levmin) return;

/* convert into xyz coordinates */
from_rtf_to_xyz(E,start_lev,AU,E->temp);

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=0;i<=neq_minus;i++)
E->temp1[m][i] = 0.0;

/* smooth in xyz coordinates */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=nels_minus;el++)
for(i=1;i<=ENODES3D;i++) {
node= E->IEN[sl_minus][m][el].node[i];
average1=average2=average3=0.0;
e1 = E->EL[sl_minus][m][el].sub[i];
for(j=1;j<=ENODES3D;j++) {
node1=E->IEN[start_lev][m][e1].node[j];
average1 += E->temp[m][E->ID[start_lev][m][node1].doff[1]];
average2 += E->temp[m][E->ID[start_lev][m][node1].doff[2]];
average3 += E->temp[m][E->ID[start_lev][m][node1].doff[3]];
}
w = weight*E->TWW[sl_minus][m][el].node[i];

E->temp1[m][E->ID[sl_minus][m][node].doff[1]] += w * average1;
E->temp1[m][E->ID[sl_minus][m][node].doff[2]] += w * average2;
E->temp1[m][E->ID[sl_minus][m][node].doff[3]] += w * average3;
}

(E->solver.exchange_id_d)(E, E->temp1, sl_minus);

for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=nno_minus;i++)  {
E->temp1[m][E->ID[sl_minus][m][i].doff[1]] *= E->MASS[sl_minus][m][i];
E->temp1[m][E->ID[sl_minus][m][i].doff[2]] *= E->MASS[sl_minus][m][i];
E->temp1[m][E->ID[sl_minus][m][i].doff[3]] *= E->MASS[sl_minus][m][i];
}

/* back into rtf coordinates */
from_xyz_to_rtf(E,sl_minus,E->temp1,AD);

return;
}

/* ================================================= */
void from_xyz_to_rtf(E,level,xyz,rtf)
struct All_variables *E;
int level;
double **rtf,**xyz;
{

int i,j,m,eqn1,eqn2,eqn3;
double cost,cosf,sint,sinf;

for (m=1;m<=E->sphere.caps_per_proc;m++)
for (i=1;i<=E->lmesh.NNO[level];i++)  {
eqn1 = E->ID[level][m][i].doff[1];
eqn2 = E->ID[level][m][i].doff[2];
eqn3 = E->ID[level][m][i].doff[3];
sint = E->SinCos[level][m][0][i];
sinf = E->SinCos[level][m][1][i];
cost = E->SinCos[level][m][2][i];
cosf = E->SinCos[level][m][3][i];
rtf[m][eqn1] = xyz[m][eqn1]*cost*cosf
+ xyz[m][eqn2]*cost*sinf
- xyz[m][eqn3]*sint;
rtf[m][eqn2] = -xyz[m][eqn1]*sinf
+ xyz[m][eqn2]*cosf;
rtf[m][eqn3] = xyz[m][eqn1]*sint*cosf
+ xyz[m][eqn2]*sint*sinf
+ xyz[m][eqn3]*cost;
}

return;
}

/* ================================================= */
void from_rtf_to_xyz(E,level,rtf,xyz)
struct All_variables *E;
int level;
double **rtf,**xyz;
{

int i,j,m,eqn1,eqn2,eqn3;
double cost,cosf,sint,sinf;

for (m=1;m<=E->sphere.caps_per_proc;m++)
for (i=1;i<=E->lmesh.NNO[level];i++)  {
eqn1 = E->ID[level][m][i].doff[1];
eqn2 = E->ID[level][m][i].doff[2];
eqn3 = E->ID[level][m][i].doff[3];
sint = E->SinCos[level][m][0][i];
sinf = E->SinCos[level][m][1][i];
cost = E->SinCos[level][m][2][i];
cosf = E->SinCos[level][m][3][i];
xyz[m][eqn1] = rtf[m][eqn1]*cost*cosf
- rtf[m][eqn2]*sinf
+ rtf[m][eqn3]*sint*cosf;
xyz[m][eqn2] = rtf[m][eqn1]*cost*sinf
+ rtf[m][eqn2]*cosf
+ rtf[m][eqn3]*sint*sinf;
xyz[m][eqn3] = -rtf[m][eqn1]*sint
+ rtf[m][eqn3]*cost;

}

return;
}

/* ========================================================== */
void fill_in_gaps(E,temp,level)
struct All_variables *E;
int level;
double **temp;
{

int i,j,k,m;
float x1,x2;
float n1,n2;
int rnoz,noxz,node0,node1,node2;
int eqn0,eqn1,eqn2;

const int dims =E->mesh.nsd;
const int ends= enodes[dims];

const int nox = E->lmesh.NOX[level];
const int noz = E->lmesh.NOZ[level];
const int noy = E->lmesh.NOY[level];
const int sl_minus = level-1;

for(m=1;m<=E->sphere.caps_per_proc;m++)       {
n1 = n2 =0.5;
noxz = nox*noz;
for(k=1;k<=noy;k+=2)          /* Fill in gaps in x direction */
for(j=1;j<=noz;j+=2)
for(i=2;i<nox;i+=2)  {
node0 = j + (i-1)*noz + (k-1)*noxz; /* this node */
node1 = node0 - noz;
node2 = node0 + noz;

/* now for each direction */

eqn0=E->ID[level][m][node0].doff[1];
eqn1=E->ID[level][m][node1].doff[1];
eqn2=E->ID[level][m][node2].doff[1];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];

eqn0=E->ID[level][m][node0].doff[2];
eqn1=E->ID[level][m][node1].doff[2];
eqn2=E->ID[level][m][node2].doff[2];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];

eqn0=E->ID[level][m][node0].doff[3];
eqn1=E->ID[level][m][node1].doff[3];
eqn2=E->ID[level][m][node2].doff[3];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
}

n1 = n2 =0.5;
for(i=1;i<=nox;i++)   /* Fill in gaps in y direction */
for(j=1;j<=noz;j+=2)
for(k=2;k<noy;k+=2)   {
node0 = j + (i-1)*noz + (k-1)*noxz; /* this node */
node1 = node0 - noxz;
node2 = node0 + noxz;

eqn0=E->ID[level][m][node0].doff[1];
eqn1=E->ID[level][m][node1].doff[1];
eqn2=E->ID[level][m][node2].doff[1];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];

eqn0=E->ID[level][m][node0].doff[2];
eqn1=E->ID[level][m][node1].doff[2];
eqn2=E->ID[level][m][node2].doff[2];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];

eqn0=E->ID[level][m][node0].doff[3];
eqn1=E->ID[level][m][node1].doff[3];
eqn2=E->ID[level][m][node2].doff[3];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
}

for(j=2;j<noz;j+=2)	  {
x1 = E->sphere.R[level][j] - E->sphere.R[level][j-1];
x2 = E->sphere.R[level][j+1] - E->sphere.R[level][j];
n1 = x2/(x1+x2);
n2 = 1.0-n1;
for(k=1;k<=noy;k++)          /* Fill in gaps in z direction */
for(i=1;i<=nox;i++)  {
node0 = j + (i-1)*noz + (k-1)*noxz; /* this node */
node1 = node0 - 1;
node2 = node0 + 1;

eqn0=E->ID[level][m][node0].doff[1];
eqn1=E->ID[level][m][node1].doff[1];
eqn2=E->ID[level][m][node2].doff[1];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];

eqn0=E->ID[level][m][node0].doff[2];
eqn1=E->ID[level][m][node1].doff[2];
eqn2=E->ID[level][m][node2].doff[2];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];

eqn0=E->ID[level][m][node0].doff[3];
eqn1=E->ID[level][m][node1].doff[3];
eqn2=E->ID[level][m][node2].doff[3];
temp[m][eqn0] = n1*temp[m][eqn1]+n2*temp[m][eqn2];
}
}
}         /* end for m */

return;
}
``````