Changeset 26236


Ignore:
Timestamp:
05/04/21 10:09:14 (4 years ago)
Author:
caronlam
Message:

adding new formalism for lovenumbers

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/cores/love_core.cpp

    r25912 r26236  
    99#include "../modules/modules.h"
    1010#include "../solutionsequences/solutionsequences.h"
    11 void love_core(FemModel* femmodel){
    12 
    13         Vector<IssmDouble> *wg    = NULL;
    14         Vector<IssmDouble> *dwdtg = NULL;
    15         IssmDouble          *x    = NULL;
    16         IssmDouble          *y    = NULL;
     11#include "petscblaslapack.h"
     12
     13
     14void love_init(FemModel* femmodel, Matlitho* matlitho){/*{{{*/
     15        // initialize Planet_Mass(r) for efficient computation of gravity, value of surface gravity and inital size of the yi equation system
     16        bool        verbosemod = (int)VerboseModule();
     17
     18        int numlayers = matlitho->numlayers;
     19        IssmDouble* r=matlitho->radius;
     20        IssmDouble* EarthMass=xNewZeroInit<IssmDouble>(numlayers+1);
     21        IssmDouble r1,r2,ro, GG;
     22        femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
     23
     24        for (int i=0;i<numlayers;i++){
     25                  r2 = r[i+1];
     26                  ro = matlitho->density[i];
     27                if (i==0){
     28                        EarthMass[i] = ro*pow(r2,3.0)*4.0*PI/3.0;
     29                }else{
     30                        r1=r[i];
     31                        EarthMass[i] = EarthMass[i-1] + ro*(pow(r2,3.0)-pow(r1,3.0))*4.0*PI/3.0;;
     32                }
     33        }
     34
     35        IssmDouble g0=EarthMass[numlayers-1]*GG/pow(r[numlayers],2.0);
     36        femmodel->parameters->SetParam(g0,LoveG0Enum);
     37        femmodel->parameters->SetParam(EarthMass,numlayers,LoveEarthMassEnum);
     38        xDelete<IssmDouble>(EarthMass);
     39
     40        int nyi=6*(numlayers+1);
     41        int starting_layer=0;
     42        femmodel->parameters->SetParam(nyi,LoveNYiEquationsEnum);
     43        femmodel->parameters->SetParam(starting_layer,LoveStartingLayerEnum);
     44}
     45
     46void GetEarthRheology(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega,FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     47
     48        //returns lame parameters (material rigity) lambda and mu for the right frequency and layer
     49        IssmDouble vi=matlitho->viscosity[layer_index];
     50        IssmDouble mu=matlitho->lame_mu[layer_index];
     51        IssmDouble la=matlitho->lame_lambda[layer_index];
     52        int rheo=matlitho->rheologymodel[layer_index];
     53
     54        if (vi!=0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid
     55                                   //Otherwise return the elastic value
     56                IssmDouble ka=la + 2.0/3.0*mu; //Bulk modulus
     57                if (rheo==2){//EBM
     58                        IssmDouble alpha=matlitho->ebm_alpha[layer_index];
     59                        IssmDouble delta=matlitho->ebm_delta[layer_index];
     60                        IssmDouble taul=matlitho->ebm_taul[layer_index];
     61                        IssmDouble tauh=matlitho->ebm_tauh[layer_index];
     62                        IssmDouble hf1,hf2,U1,U2;
     63                        IssmDouble* a=xNew<IssmDouble>(2);
     64                        IssmDouble* b=xNew<IssmDouble>(1);
     65
     66                        a[0]=1.0;a[1]=1.0+alpha/2.0;
     67                        b[0]=2.0+alpha/2.0;
     68                        //hf1=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*tauh,2.0), 0, 0, 15);
     69                        //hf2=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*taul,2.0), 0, 0, 15);
     70                        //hf1=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*tauh,2.0));
     71                        //hf2=hypergeometric_pFq_({a[0], a[1]}, {b[0]},-pow(omega*taul,2.0));
     72                        U1=(pow(tauh,alpha)-pow(taul,alpha))/alpha-pow(omega,2.0)/(2.0+alpha)*(pow(tauh,2.0+alpha)*hf1-pow(taul,2.0+alpha)*hf2);
     73
     74                        a[0]=1.0;a[1]=.5+alpha/2.0;
     75                        b[0]=1.5+alpha/2.0;
     76                        //hf1=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*tauh,2.0), 0, 0, 15);
     77                        //hf2=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*taul,2.0), 0, 0, 15);
     78                        //hf1=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*tauh,2.0));
     79                        //hf2=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*taul,2.0));
     80                        U2=(pow(tauh,1.0+alpha)*hf1-pow(taul,1.0+alpha)*hf2)/(1.0+alpha);
     81
     82                        mu=mu/(1.0+ alpha*delta/(pow(tauh,alpha)-pow(taul,alpha))*(U1 + omega*U2) -mu/vi/omega);
     83                        la=ka-2.0/3.0*mu;
     84
     85                        xDelete<IssmDouble>(a);
     86                        xDelete<IssmDouble>(b);
     87                }
     88                else if (rheo==1){//Burgers
     89                        IssmDouble vi2=matlitho->burgers_viscosity[layer_index];
     90                        IssmDouble mu2=matlitho->burgers_mu[layer_index];
     91
     92                        mu=mu*omega*(omega+mu2/vi2)/((omega+mu2/vi2)*(omega+mu/vi)+mu/vi2*omega);
     93                        la=ka-2.0/3.0*mu;
     94                }
     95                else{//Maxwell
     96                        la = (la + mu*ka/vi/omega)/(1.0 + mu/vi/omega);
     97                        mu = mu/(1.0+mu/vi/omega);
     98                }
     99         }
     100
     101         *pla=la;
     102         *pmu=mu;
     103
     104} /*}}}*/
     105
     106
     107IssmDouble GetGravity(IssmDouble r2, int layer_index, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     108        //computes gravity at radius r2
     109        IssmDouble* EarthMass;
     110        IssmDouble g, GG;
     111
     112        femmodel->parameters->FindParam(&EarthMass,&matlitho->numlayers,LoveEarthMassEnum);
     113        femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
     114        IssmDouble ro=matlitho->density[layer_index];
     115        IssmDouble M=0;
     116        IssmDouble r1=0;
     117        if (layer_index==0){
     118                M=4.0/3.0*PI*ro*pow(r2,3.0);
     119        }
     120        else{
     121                r1=matlitho->radius[layer_index];
     122                M=EarthMass[layer_index-1]+4.0/3.0*PI*ro*(pow(r2,3.0)-pow(r1,3.0));
     123        }
     124        return  g= GG*M/pow(r2,2.0);
     125}/*}}}*/
     126
     127void fill_yi_prefactor(IssmDouble* yi_prefactor, int* pdeg, IssmDouble* pomega, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     128        //precalculates partial derivative factors for function yi_derivatives
     129 IssmDouble ra=matlitho->radius[matlitho->numlayers];
     130 IssmDouble  g0,r0,mu0, GG;
     131 int nstep,nindex, starting_layer;
     132 femmodel->parameters->FindParam(&g0,LoveG0Enum);
     133 femmodel->parameters->FindParam(&r0,LoveR0Enum);
     134 femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
     135 femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
     136 femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     137 femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
     138
     139 IssmDouble frh,frhg0,fgr0,fgr,fn,rm0,rlm,flm;
     140 IssmDouble xmin,xmax,x,dr;
     141 IssmDouble g,ro, issolid;
     142
     143        if (pomega) { //frequency and degree dependent terms /*{{{*/
     144                IssmDouble la,mu;
     145                IssmDouble* f=xNewZeroInit<IssmDouble>(12);
     146                int deg=*pdeg;
     147                IssmDouble omega=*pomega;       
     148                fn=deg*(deg+1.0);
     149
     150                for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
     151
     152                                ro=matlitho->density[layer_index];
     153                                issolid=matlitho->issolid[layer_index];
     154                                if(issolid==1){
     155                                GetEarthRheology(&la, &mu,layer_index,omega,femmodel,matlitho);   
     156       
     157                                /*_______Expressions*/
     158                                flm=(la+2.0*mu);
     159                                rlm=(3.0*la+2.0*mu)/(la+2.0*mu);
     160                                rm0=mu/mu0;
     161                                frh=ro*ra/mu0;
     162
     163                                f[0]=(-2.0*la/flm);
     164                                f[1]=mu0/flm;
     165                                f[2]=(la*fn/flm);
     166                                f[3]=rm0*rlm;
     167                                f[4]=ro*pow(omega,2.0)*ra/mu0;
     168                                f[5]=(-4.0*mu/flm);
     169                                f[6]=fn*frh;
     170                                f[7]=-(2.0*rm0*rlm)*fn;
     171                                f[8]=1.0/rm0;
     172                                f[9]=-2.0*rm0*rlm;
     173                                f[10]=-la/flm;
     174                                f[11]=2.0*rm0*(la*(2.0*fn-1.0)+2.0*mu*(fn-1.0))/flm;
     175
     176                                xmin=matlitho->radius[layer_index]/ra;
     177                                xmax=(matlitho->radius[layer_index+1])/ra;
     178                                dr = (xmax -xmin)/nstep;
     179                                x=xmin;
     180                                for (int n=0;n<nstep;n++){
     181                                        g=GetGravity(x*ra,layer_index,femmodel,matlitho);
     182
     183                                        nindex=layer_index*nstep*36+n*36;
     184                                        yi_prefactor[nindex+ 0*6+0]= f[0]/x;                      // in dy[0*6+0]
     185                                        yi_prefactor[nindex+ 0*6+1]= f[1];                        // in dy[0*6+1]
     186                                        yi_prefactor[nindex+ 0*6+2]= f[2]/x;                      // in dy[0*6+2]
     187                                        yi_prefactor[nindex+ 1*6+0]= 4.0*(-frh*g+f[3]/x)/x + f[4];// in dy[1*6+0]
     188                                        yi_prefactor[nindex+ 1*6+1]= f[5]/x;                      // in dy[1*6+1]
     189                                        yi_prefactor[nindex+ 1*6+2]= (f[6]*g+f[7]/x)/x;           // in dy[1*6+2]
     190                                        yi_prefactor[nindex+ 2*6+3]= f[8];                        // in dy[2*6+3]
     191                                        yi_prefactor[nindex+ 3*6+0]= (frh*g+f[9]/x)/x;            // in dy[3*6+0]
     192                                        yi_prefactor[nindex+ 3*6+1]= f[10]/x;                     // in dy[3*6+1]
     193                                        yi_prefactor[nindex+ 3*6+2]= f[11]/(x*x) + f[4];          // in dy[3*6+2]
     194                                        x=x+dr;
     195                                }
     196                        }
     197                }
     198
     199                xDelete<IssmDouble>(f);
     200        /*}}}*/
     201        } else if (pdeg) { // degree dependent terms /*{{{*/
     202                int deg=*pdeg;
     203                fn=(deg*(deg+1.0));
     204
     205                for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
     206                        ro=matlitho->density[layer_index];
     207                        issolid=matlitho->issolid[layer_index];
     208       
     209                        /*_______Expressions*/
     210                        fgr=4.0*PI*GG*ro*ra;
     211         
     212                        xmin=matlitho->radius[layer_index]/ra;
     213                        xmax=(matlitho->radius[layer_index+1])/ra;
     214                        dr = (xmax -xmin)/nstep;
     215                        x=xmin;
     216                        for (int n=0;n<nstep;n++){
     217                                 nindex=layer_index*nstep*36+n*36;
     218                                 g=GetGravity(x*ra,layer_index,femmodel,matlitho);
     219       
     220                                 if(issolid==1){
     221                                 yi_prefactor[nindex+ 1*6+3]= fn/x;                  // in dy[1*6+3]
     222                                 yi_prefactor[nindex+ 5*6+2]= -(fgr/g0*fn)/x;        // in dy[5*6+2]
     223                                 yi_prefactor[nindex+ 5*6+4]= fn/(x*x);              // in dy[5*6+4]
     224                                 } else {
     225                                 yi_prefactor[nindex+ 1*6+0]= (-4.0*(fgr/g)+fn/x)/x; // in dy[1*6+0] liquid layer
     226                                 }
     227                                x=x+dr;
     228                        }
     229                }
     230        /*}}}*/
     231        } else { // static terms /*{{{*/
     232                for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
     233                        ro=matlitho->density[layer_index];
     234                        issolid=matlitho->issolid[layer_index];
     235       
     236                        /*_______Expressions*/
     237                        frhg0=ro*g0*ra/mu0;
     238                        fgr=4.0*PI*GG*ro*ra;
     239
     240                        xmin=matlitho->radius[layer_index]/ra;
     241                        xmax=(matlitho->radius[layer_index+1])/ra;
     242                        dr = (xmax -xmin)/nstep;
     243                        x=xmin;
     244                        for (int n=0;n<nstep;n++){
     245                                 g=GetGravity(x*ra,layer_index,femmodel,matlitho);
     246                                 nindex=layer_index*nstep*36+n*36;
     247                                 if(issolid==1){
     248                                         yi_prefactor[nindex+ 1*6+5]= -frhg0;       // in dy[1*6+5]
     249                                         yi_prefactor[nindex+ 2*6+0]= -1.0/x;       // in dy[2*6+0]
     250                                         yi_prefactor[nindex+ 2*6+2]= 1.0/x;        // in dy[2*6+2]
     251                                         yi_prefactor[nindex+ 3*6+3]= -3.0/x;       // in dy[3*6+3]
     252                                         yi_prefactor[nindex+ 3*6+4]= -frhg0/x;     // in dy[3*6+4]
     253                                         yi_prefactor[nindex+ 4*6+0]= fgr/g0;       // in dy[4*6+0]
     254                                         yi_prefactor[nindex+ 4*6+5]= 1.0;          // in dy[4*6+5]
     255                                         yi_prefactor[nindex+ 5*6+5]= -2.0/x;       // in dy[5*6+5]
     256                                 } else {
     257                                         yi_prefactor[nindex+ 0*6+0]= fgr/g;        // in dy[0*6+0] liquid layer
     258                                         yi_prefactor[nindex+ 0*6+1]= 1.0;          // in dy[0*6+1] liquid layer
     259                                         yi_prefactor[nindex+ 1*6+1]= -2.0/x-fgr/g; // in dy[1*6+1] liquid layer
     260                                 }
     261                                x=x+dr;
     262                        }
     263                }
     264        /*}}}*/
     265        }
     266}/*}}}*/
     267
     268IssmDouble* yi_derivatives(IssmDouble* y, int layer_index, int n, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     269        //computes yi derivatives at r=radius[layer_index]+ n/nstep*(radius[layer_index+1]-radius[layer_index])
     270
     271        IssmDouble issolid=matlitho->issolid[layer_index];
     272        int iy,id,ny, nindex, nstep;
     273        IssmDouble* dydx=xNewZeroInit<IssmDouble>(6);
     274        IssmDouble* dy=xNewZeroInit<IssmDouble>(36);
     275        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     276
     277        /*{{{*/ /* For reference:
     278         flm=(la+2.0*mu);
     279         rlm=(3.0*la+2.0*mu)/(la+2.0*mu);
     280         rm0=mu/mu0;
     281         rg0=g/g0;
     282         frh=ro*g*ra/mu0;
     283         fgr=4.0*PI*GG*ro*ra/g0;
     284         fn=(deg*(deg+1.0));
     285         
     286
     287         if(issolid==1){
     288         ny = 6;
     289
     290         dy[0*6+0]= (-2.0*la/flm)/x;
     291         dy[0*6+1]= mu0/flm;
     292         dy[0*6+2]= (la*fn/flm)/x;
     293         dy[0*6+3]= 0.0;
     294         dy[0*6+4]= 0.0;
     295         dy[0*6+5]= 0.0;
     296
     297         dy[1*6+0]=  4.0*(-frh+rm0*rlm/x)/x + ro*pow(omega,2.0)*ra/mu0;
     298         dy[1*6+1]=(-4.0*mu/flm)/x;
     299         dy[1*6+2]= fn*(frh-2.0*rm0*rlm/x)/x;
     300         dy[1*6+3]= fn/x;
     301         dy[1*6+4]= 0.0;
     302         dy[1*6+5]= -frh/rg0;
     303
     304         dy[2*6+0]= -1.0/x;
     305         dy[2*6+1]= 0.0;
     306         dy[2*6+2]= 1.0/x;
     307         dy[2*6+3]= 1/rm0;
     308         dy[2*6+4]= 0.0;
     309         dy[2*6+5]= 0.0;
     310
     311         dy[3*6+0]= (frh-2.0*rm0*rlm/x)/x;
     312         dy[3*6+1]= ( -la/flm)/x;
     313         dy[3*6+2]= (2.0*rm0*(la*(2.0*fn-1.0)+2.0*mu*(fn-1.0))/flm)/(x*x) + ro*pow(omega,2.0)*ra/mu0;
     314         dy[3*6+3]= -3.0/x;
     315         dy[3*6+4]= -(frh/rg0)/x;
     316         dy[3*6+5]= 0.0;
     317
     318         dy[4*6+0]= fgr;
     319         dy[4*6+1]= 0.0;
     320         dy[4*6+2]= 0.0;
     321         dy[4*6+3]= 0.0;
     322         dy[4*6+4]= 0.0;
     323         dy[4*6+5]= 1.0;
     324
     325         dy[5*6+0]= 0.0;
     326         dy[5*6+1]= 0.0;
     327         dy[5*6+2]= -(fgr*fn)/x;
     328         dy[5*6+3]= 0.0;
     329         dy[5*6+4]= fn/(x*x);
     330         dy[5*6+5]= -2.0/x;
     331         
     332         } else {
     333         ny = 2;
     334
     335         dy[0*6+0]= fgr/rg0;
     336         dy[0*6+1]= 1.0;
     337         dy[1*6+0]= (-4.0*(fgr/rg0)+fn/x)/x;
     338         dy[1*6+1]= -2.0/x-fgr/rg0;
     339        }
     340        */ /*}}}*/
     341
     342         if(issolid==1){
     343                 ny = 6;
     344         } else {
     345                 ny = 2;
     346         }
     347         
     348         for (id=0;id<ny;id++){
     349                for (iy=0;iy<ny;iy++){
     350                        nindex=layer_index*nstep*36+n*36;
     351                        dydx[id]=dydx[id]+yi_prefactor[nindex+id*6+iy]*y[iy];
     352                }
     353         }
     354        xDelete<IssmDouble>(dy);
     355        return dydx;
     356}/*}}}*/
     357
     358void propagate_yi_euler(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     359        //yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer?
     360        //euler method
     361        int nstep;
     362        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     363
     364        IssmDouble* dydx=xNewZeroInit<IssmDouble>(6);
     365        IssmDouble dr = (xmax -xmin)/nstep;
     366        IssmDouble x=xmin;
     367    for(int i = 0;i<nstep;i++){
     368        dydx=yi_derivatives(y,layer_index, i,yi_prefactor,femmodel,matlitho);
     369        for (int j=0;j<6;j++){
     370                y[j]+=dydx[j]*dr;
     371        }
     372        x = x + dr;
     373    }
     374        xDelete<IssmDouble>(dydx);
     375}/*}}}*/
     376
     377void propagate_yi_RK2(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     378        //yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer?
     379        //Implements Runge-Kutta 2nd order (midpoint method)
     380        int nstep;
     381        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     382
     383        IssmDouble* k1=xNewZeroInit<IssmDouble>(6);
     384        IssmDouble* k2=xNewZeroInit<IssmDouble>(6);
     385        IssmDouble* k3=xNewZeroInit<IssmDouble>(6);
     386        IssmDouble* k4=xNewZeroInit<IssmDouble>(6);
     387        IssmDouble* y1=xNewZeroInit<IssmDouble>(6);
     388        IssmDouble* y2=xNewZeroInit<IssmDouble>(6);
     389        IssmDouble* y3=xNewZeroInit<IssmDouble>(6);
     390
     391        IssmDouble dr = (xmax -xmin)/nstep;
     392        IssmDouble x=xmin;
     393    for(int i = 0;i<nstep/2;i++){
     394        k1=yi_derivatives(y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
     395        for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
     396        k2=yi_derivatives(y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);               
     397       
     398        for (int j=0;j<6;j++){
     399                  y[j]+=k2[j]*2.0*dr;
     400        }
     401        x = x + 2.0*dr;
     402    }
     403        xDelete<IssmDouble>(k1);
     404        xDelete<IssmDouble>(k2);
     405        xDelete<IssmDouble>(k3);
     406        xDelete<IssmDouble>(k4);
     407        xDelete<IssmDouble>(y1);
     408        xDelete<IssmDouble>(y2);
     409        xDelete<IssmDouble>(y3);
     410}/*}}}*/
     411
     412void propagate_yi_RK4(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     413        //yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer?
     414        //Implements Runge-Kutta 4th order
     415        int nstep;
     416        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     417
     418        IssmDouble* k1=xNewZeroInit<IssmDouble>(6);
     419        IssmDouble* k2=xNewZeroInit<IssmDouble>(6);
     420        IssmDouble* k3=xNewZeroInit<IssmDouble>(6);
     421        IssmDouble* k4=xNewZeroInit<IssmDouble>(6);
     422        IssmDouble* y1=xNewZeroInit<IssmDouble>(6);
     423        IssmDouble* y2=xNewZeroInit<IssmDouble>(6);
     424        IssmDouble* y3=xNewZeroInit<IssmDouble>(6);
     425
     426        IssmDouble dr = (xmax -xmin)/nstep;
     427        IssmDouble x=xmin;
     428    for(int i = 0;i<nstep/2-1;i++){
     429        k1=yi_derivatives(y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
     430        for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
     431        k2=yi_derivatives(y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
     432        for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;}
     433        k3=yi_derivatives(y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
     434        for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;}
     435        k4=yi_derivatives(y3,layer_index, 2*i+2,yi_prefactor,femmodel,matlitho);               
     436       
     437        for (int j=0;j<6;j++){
     438                y[j]+=(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/3.0*dr;         
     439        }
     440        x = x + 2.0*dr;
     441    }
     442
     443        //Last step: we don't know the derivative at xmax, so we will assume the values at xmax-dr
     444        int i=nstep/2;
     445        k1=yi_derivatives(y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
     446        for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
     447        k2=yi_derivatives(y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
     448        for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;}
     449        k3=yi_derivatives(y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
     450        for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;}
     451        k4=yi_derivatives(y3,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);               
     452       
     453        for (int j=0;j<6;j++){
     454                y[j]+=(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/3.0*dr;         
     455        }
     456
     457        x = x + 2.0*dr;
     458
     459        xDelete<IssmDouble>(k1);
     460        xDelete<IssmDouble>(k2);
     461        xDelete<IssmDouble>(k3);
     462        xDelete<IssmDouble>(k4);
     463        xDelete<IssmDouble>(y1);
     464        xDelete<IssmDouble>(y2);
     465        xDelete<IssmDouble>(y3);
     466}/*}}}*/
     467
     468void Innersphere_boundaryconditions(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     469        //fills the boundary conditions at the bottom of layer[layer_index] in yi[0:2][0:5]
     470
     471        IssmDouble r = matlitho->radius[layer_index];
     472        IssmDouble ra=matlitho->radius[matlitho->numlayers];
     473        IssmDouble  g0,r0,mu0, GG;
     474        int nyi;
     475        femmodel->parameters->FindParam(&g0,LoveG0Enum);
     476        femmodel->parameters->FindParam(&r0,LoveR0Enum);
     477        femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
     478        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
     479        femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
     480        IssmDouble ro=matlitho->density[layer_index];
     481        IssmDouble issolid=matlitho->issolid[layer_index];
     482        IssmDouble g=GetGravity(r,layer_index,femmodel,matlitho);
     483        IssmDouble la,mu;
     484
     485        if (layer_index==0){
     486                // radius[0] cannot be 0 for numerical reasons, but below our first interface at radius[0] would in reality be the same material as in the first layer
     487                GetEarthRheology(&la, &mu,layer_index,omega,femmodel,matlitho);   
     488        } else {
     489                GetEarthRheology(&la, &mu,layer_index-1,omega,femmodel,matlitho);   
     490        }   
     491
     492        IssmDouble cst = 4.0*PI*GG*ro;
     493        IssmDouble r2=pow(r,2.0);
     494
     495        yi[0+nyi*0]=1.0*r/ra;
     496        yi[0+nyi*1]=1.0/(r*ra);
     497        yi[0+nyi*2]=0.0;
     498
     499        yi[1+nyi*0]=(2.0*mu*(deg-1.0-3.0/deg) + cst/3.0*ro*r2)/mu0;
     500        yi[1+nyi*1]=(2.0*mu*(deg-1.0)/r2 + cst/3.0*ro)/mu0;
     501        yi[1+nyi*2]=-ro/mu0;
     502
     503        yi[2+nyi*0]=(deg+3.0)/(deg*(deg+1.0))*r/ra;
     504        yi[2+nyi*1]=1.0/(deg*r*ra);
     505        yi[2+nyi*2]=0.0;
     506
     507        yi[3+nyi*0]=2.0*mu*(deg+2.0)/((deg+1.0)*mu0);
     508        yi[3+nyi*1]=2.0*mu*(deg-1.0)/(deg*r2*mu0);
     509        yi[3+nyi*2]=0.0;
     510
     511        yi[4+nyi*0]=0.0;
     512        yi[4+nyi*1]=0.0;
     513        yi[4+nyi*2]=1.0/(g0*ra);
     514
     515        yi[5+nyi*0]=-cst*r/g0;
     516        yi[5+nyi*1]=-cst/(r*g0);
     517        yi[5+nyi*2]=deg/(r*g0);
     518}/*}}}*/
     519
     520void build_yi_system(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho) { /*{{{*/
     521
     522        IssmDouble  g0,r0,mu0,x,ro1, GG;
     523        int nyi,starting_layer, nstep;
     524        femmodel->parameters->FindParam(&g0,LoveG0Enum);
     525        femmodel->parameters->FindParam(&r0,LoveR0Enum);
     526        femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
     527        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
     528        femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
     529        femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
     530        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     531
     532          IssmDouble xmin,xmax,one,ro,g;
     533          IssmDouble ra=matlitho->radius[matlitho->numlayers];
     534
     535        for (int i=0;i<6*(matlitho->numlayers+1);i++){
     536                for(int j=0;j<6*(matlitho->numlayers+1);j++){
     537                        yi[i+6*(matlitho->numlayers+1)*j]=0.0;
     538                }
     539        }
     540
     541        int ny,is,ii,jj;
     542        IssmDouble ystart[6];
     543        for (int k=0;k<6;k++) ystart[k]=0.0;           
     544 
     545        int ici = 0;   // Index of current interface
     546        for (int i = starting_layer; i<matlitho->numlayers;i++){
     547
     548                xmin=matlitho->radius[i]/ra;
     549                xmax=(matlitho->radius[i+1])/ra;
     550
     551                 
     552                if (matlitho->issolid[i]){
     553                                ny = 6;
     554                                is = 0;
     555                                one= 1.0;
     556                } else {       
     557                                ny = 2;
     558                                is = 4;
     559                                one= -1.0;
     560                }
     561       
     562
     563                for (int j = 0;j<ny;j++){
     564                        for (int k=0;k<6;k++){ystart[k]=0.0;}
     565                        ystart[j]= 1.0;
     566                                       
     567                        // Numerical Integration
     568                        //propagate_yi_euler(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
     569                        propagate_yi_RK2(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
     570                        //propagate_yi_RK4(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
     571                        // Boundary Condition matrix - propagation part
     572                        ii = 6*(ici+1)+is;
     573                        jj = 6*(ici+1)+j+is-3;
     574                        for (int kk=0;kk<ny;kk++){
     575                                yi[(ii+kk)+nyi*jj] = ystart[kk]*one;
     576                        }
     577                }
     578                 
     579                 
     580                // Boundary Condition matrix - solid regions
     581                if (matlitho->issolid[i]){
     582                        one = -1.0;
     583                        if (i>0 & !matlitho->issolid[i-1]) one = 1.0;
     584                        for (int j=0;j<6;j++){
     585                                yi[(j+6*ici)+ nyi*(j+6*ici+3)] = one;
     586                        }
     587                } else { // Boundary Condition matrix - liquid regions
     588                        ro1=matlitho->density[i];
     589                        g=GetGravity(matlitho->radius[i], i, femmodel,matlitho);
     590                        ii = 6*ici;
     591                        yi[ii+nyi*(ii+3)] = -1.0;
     592                        yi[ii+nyi*(ii+4+3)] = -g0/g;
     593                        yi[(ii+1)+nyi*(ii+3)]=-ro1*g*ra/mu0;
     594                        yi[(ii+2)+nyi*(ii+1+3)]=-1.0;
     595                        yi[(ii+5)+nyi*(ii+3)]= 4.0*PI*GG*ro1*ra/g0;
     596                        yi[(ii+4)+nyi*(ii+4+3)]=-1.0;
     597                        yi[(ii+5)+nyi*(ii+5+3)]=-1.0;
     598                        g=GetGravity(matlitho->radius[i+1], i,femmodel,matlitho);
     599                        ii = 6*(ici+1);
     600                        yi[ii+nyi*(ii-1)]=-1.0;
     601                        yi[ii+nyi*(ii+1)]=yi[(ii+4)+nyi*(ii+1)]*g0/g; // yi(17,14) solution integration 1 of z5 CMB
     602                        yi[ii+nyi*(ii+2)]=yi[(ii+4)+nyi*(ii+2)]*g0/g; // yi(17,15) solution integration 2 of z5 CMB
     603                          // yi(13,..) y1 CMB
     604                        yi[(ii+1)+nyi*(ii-1)]=-ro1*g*ra/mu0;
     605                        yi[(ii+2)+nyi*(ii)]=-1.0;
     606                        yi[(ii+5)+nyi*(ii-1)]= 4.0*PI*GG*ro1*ra/g0;
     607                }       
     608                ici = ici+1;
     609        }
     610
     611
     612        //-- Internal sphere: integration starts here rather than r=0 for numerical reasons
     613        Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho);
     614
     615        //-- Surface conditions
     616        yi[(nyi-6)+nyi*(nyi-3)]=-1.0;
     617        yi[(nyi-4)+nyi*(nyi-2)]=-1.0;
     618        yi[(nyi-2)+nyi*(nyi-1)]=-1.0;
     619        yi[(nyi-1)+nyi*(nyi-1)]=deg+1.0;
     620
     621        //-- Degree 1 special case
     622        if(deg==1){
     623                for (int i=0;i<nyi;i++){
     624                        yi[(nyi-1)+nyi*i]=0.0;
     625                }
     626                yi[(nyi-1)+nyi*(nyi-1)]=1.0;
     627        }
     628}/*}}}*/
     629
     630void yi_boundary_conditions(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     631
     632        IssmDouble  g0,r0,mu0,ra,rb,rc;
     633        int nyi, forcing_type,icb,cmb;
     634        IssmDouble* EarthMass;
     635        femmodel->parameters->FindParam(&g0,LoveG0Enum);
     636        femmodel->parameters->FindParam(&r0,LoveR0Enum);
     637        femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
     638        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
     639        femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
     640        femmodel->parameters->FindParam(&icb,LoveInnerCoreBoundaryEnum);
     641        femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum);
     642        femmodel->parameters->FindParam(&EarthMass,&matlitho->numlayers,LoveEarthMassEnum);
     643         
     644        // In Case of a Inner Core - Outer Core - Mantle planet and Boundary conditions on these 3 interfaces
     645        ra=matlitho->radius[matlitho->numlayers];       
     646        rb=0;
     647        rc=0;
     648        if (forcing_type<=4){
     649                rc=matlitho->radius[icb];
     650        }
     651        else if (forcing_type<=8){
     652                rb=matlitho->radius[cmb];
     653        }
     654
     655        IssmDouble ro_mean=EarthMass[matlitho->numlayers-1]/(4.0/3.0*PI*pow(ra,3.0));
     656        xDelete<IssmDouble>(EarthMass);
     657
     658        for (int i=0;i<(matlitho->numlayers+1)*6;i++) yi_righthandside[i]=0.0;
     659
     660        switch (forcing_type) {
     661
     662        //-- forcings at the Inner Core Boundary
     663        case 1: //'ICB --Volumetric Potential'
     664                yi_righthandside[6*icb+5]=(deg)/(rc*g0);
     665                yi_righthandside[6*icb+4]=1.0/(ra*g0);
     666                break;
     667        case 2: //'ICB --Pressure'
     668                yi_righthandside[6*icb+1]=-ro_mean/mu0;
     669                break;
     670        case 3://'ICB --Loading'
     671                yi_righthandside[6*icb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rc;
     672                yi_righthandside[6*icb+5]= (2.0*deg+1.0)/(rc*g0);
     673                break;
     674        case 4://'ICB --Tangential Traction'
     675                yi_righthandside[6*icb+3]= ro_mean/mu0;
     676                break;
     677
     678        //--forcings at the Core Mantle Boundary
     679        case 5://'CMB --Volumetric Potential'
     680                yi_righthandside[6*cmb+1]=-ro_mean/mu0*ra/rb;
     681                yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0);
     682                break;
     683        case 6://'CMB --Pressure'
     684                yi_righthandside[6*cmb+1]=-ro_mean/mu0;
     685                break;
     686        case 7://'CMB --Loading'
     687                yi_righthandside[6*cmb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rb;
     688                yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0);
     689                break;
     690        case 8://'CMB --Tangential Traction'
     691                yi_righthandside[6*cmb+3]=-ro_mean/mu0;
     692                break;
     693
     694        //--forcings at the surface
     695        case 9://'SURF--Volumetric Potential'
     696                if (deg>1) yi_righthandside[nyi-1]=(2.0*deg+1.0)/(ra*g0);
     697                break;
     698        case 10://'SURF--Pressure'
     699                yi_righthandside[nyi-5]=-ro_mean/mu0;
     700                break;
     701        case 11://'SURF--Loading'
     702                yi_righthandside[nyi-5]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0);
     703                if (deg>1) yi_righthandside[nyi-1]= (2.0*deg+1.0)/(ra*g0);
     704                break;
     705        case 12://'SURF--Tangential Traction'
     706                yi_righthandside[nyi-3]= ro_mean/mu0;
     707                break;
     708        default:
     709                _error_("love core error: forcing_type not supported yet");
     710        }
     711}/*}}}*/
     712
     713
     714void solve_yi_system(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     715
     716        IssmDouble  g0,r0,mu0,loveratio,underflow_tol;
     717        IssmDouble* frequencies;
     718        int nyi,starting_layer, dummy;
     719        bool allow_layer_deletion;
     720        femmodel->parameters->FindParam(&g0,LoveG0Enum);
     721        femmodel->parameters->FindParam(&r0,LoveR0Enum);
     722        femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
     723        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
     724        femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
     725        femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
     726        femmodel->parameters->FindParam(&underflow_tol,LoveUnderflowTolEnum);
     727        femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum);
     728        IssmDouble ra=matlitho->radius[matlitho->numlayers];
     729        bool exit=false;
     730        int lda,ldb;
     731
     732        for(;!exit;){ //cycles of: attempt to solve the yi system, then delete a layer if necessary
     733                lda=nyi;
     734                ldb=nyi;
     735                IssmDouble*  yilocal=xNew<IssmDouble>(nyi*nyi); // we will need to redeclare these inside here as nyi changes
     736                IssmDouble*  rhslocal=xNew<IssmDouble>(nyi);
     737
     738                //we need to do a local copy of yi,rhs to send them to LAPACK with the appropriate size and to keep the original matrices in case layers are deleted and we need to try again
     739                for (int i=0;i<nyi;i++){
     740                        rhslocal[i]=rhs[i];
     741                        for (int j=0;j<nyi;j++){
     742                                yilocal[i+j*nyi]=yi[i+j*nyi];
     743                        }
     744                }
     745                //-- Resolution
     746                int* ipiv=xNewZeroInit<int>(nyi); //pivot index vector
     747                int info = 0;// error checker
     748                int nrhs=1; // number of right hand size columns
     749
     750                dgesv_(&nyi, &nrhs, yilocal, &lda, ipiv, rhslocal, &ldb, &info);
     751                xDelete<int>(ipiv);
     752
     753                if(VerboseSolution() && info!=0){
     754                        _printf0_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system");
     755                        printf("%s %s %s %s\n", "i","j","yi[i+nyi*j]","rhs[i]");
     756                        for (int i=0;i<nyi;i++){
     757                                for (int j=0;j<nyi;j++){
     758                                        printf("%i %i %e %e\n", i,j,yi[i+nyi*j],rhs[i]);
     759                                }
     760                        }
     761                }
     762
     763                *loveh = rhslocal[nyi-3]*ra*g0;
     764                *lovel = rhslocal[nyi-2]*ra*g0;
     765                *lovek = rhslocal[nyi-1]*ra*g0;
     766
     767                IssmDouble loveh1 = rhslocal[3];
     768                IssmDouble lovel1 = rhslocal[5];
     769                IssmDouble lovek1 = rhslocal[7] - pow(matlitho->radius[starting_layer]/ra,deg)/(g0*ra);
     770
     771                IssmDouble loveh1s = rhslocal[nyi-3];
     772                IssmDouble lovel1s = rhslocal[nyi-2];
     773                IssmDouble lovek1s = rhslocal[nyi-1] - 1.0/(g0*ra);
     774
     775                loveratio = abs(loveh1/loveh1s); //ratio of center to surface love numbers, determines if we should remove layers
     776                if (abs(lovel1/lovel1s) < loveratio) loveratio = abs(lovel1/lovel1s);
     777                if (abs(lovek1/lovek1s) < loveratio) loveratio = abs(lovek1/lovek1s);
     778
     779                if (!allow_layer_deletion || nyi<=12 || omega!=2.0*PI*frequencies[0]){ //We are not allowed to delete layers, or there is only one layer left. We also don't want to delete layers in the middle of a loop on frequencies, as that can lead to a jump that would compromise the inverse laplace transform.
     780                        goto save_results;
     781                }
     782
     783                if (loveratio<=underflow_tol || xIsNan(loveratio) || xIsInf(loveratio)){//We need to delete a layer and try again if the ratio between deepest love number to surface love number is too low (risk of underflow) or garbage
     784                        if(VerboseSolution()){
     785                                printf("\n   Degree: %i, surface/Depth Love number ratio small: %e\n",deg, loveratio);
     786                                printf("    Changing the interface where the integration starts\n");
     787                                printf("    New start interface: r=%gm\n\n", matlitho->radius[starting_layer+1]);
     788                        }
     789                        nyi-=6;
     790                        starting_layer+=1;
     791                        femmodel->parameters->SetParam(nyi,LoveNYiEquationsEnum);
     792                        femmodel->parameters->SetParam(starting_layer,LoveStartingLayerEnum);
     793
     794
     795                        for (int i=0;i<nyi;i++){//shift everything down by 1 layer
     796                                rhs[i]=rhs[i+6];
     797                                for (int j=0;j<nyi;j++){
     798                                        yi[j+i*nyi]=yi[j+6+(i+6)*(nyi+6)];
     799                                }
     800                        }
     801
     802                        Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho); //we move the first interface to the new starting layer. yi[0:2,0:5] will be different
     803                } else { //we are ready to save the outputs and break the main loop
     804
     805                save_results:
     806                        for (int i=0;i<nyi;i++){
     807                                rhs[i]=rhslocal[i];
     808                                for (int j=0;j<nyi;j++){
     809                                        yi[j+i*nyi]=yilocal[j+i*nyi];
     810                                }
     811                        }
     812
     813                        //make sure we can't output numbers from deleted layers
     814                        for (int i=nyi;i<(matlitho->numlayers+1)*6;i++){
     815                                rhs[i]=0.0;
     816                                for (int j=0;j<(matlitho->numlayers+1)*6;j++){
     817                                        yi[j+i*(matlitho->numlayers+1)*6]=0.0;
     818                                }
     819                        }
     820                        for (int i=0;i<nyi;i++){
     821                                for (int j=nyi;j<(matlitho->numlayers+1)*6;j++){
     822                                        yi[j+i*(matlitho->numlayers+1)*6]=0.0;
     823                                }
     824                        }
     825
     826                        exit = true;
     827                }
     828                xDelete<IssmDouble>(yilocal);
     829                xDelete<IssmDouble>(rhslocal);
     830        }
     831        xDelete<IssmDouble>(frequencies);       
     832}/*}}}*/
     833
     834IssmDouble factorial(int n){ /*{{{*/
     835        IssmDouble prod=1;
     836        for (int i=2;i<n+1;i++) prod*=i;
     837        return prod;
     838}/*}}}*/
     839
     840IssmDouble n_C_r(int n, int r){ /*{{{*/
     841        //n choose r
     842        int primes[169] =
     843  {2,    3,    5,    7,   11,   13,   17,   19,   23,   29, 
     844  31,   37,   41,   43,   47,   53,   59,   61,   67,   71, 
     845  73,   79,   83,   89,   97,  101,  103,  107,  109,  113, 
     846 127,  131,  137,  139,  149,  151,  157,  163,  167,  173, 
     847 179,  181,  191,  193,  197,  199,  211,  223,  227,  229, 
     848 233,  239,  241,  251,  257,  263,  269,  271,  277,  281, 
     849 283,  293,  307,  311,  313,  317,  331,  337,  347,  349, 
     850 353,  359,  367,  373,  379,  383,  389,  397,  401,  409, 
     851 419,  421,  431,  433,  439,  443,  449,  457,  461,  463, 
     852 467,  479,  487,  491,  499,  503,  509,  521,  523,  541, 
     853 547,  557,  563,  569,  571,  577,  587,  593,  599,  601, 
     854 607,  613,  617,  619,  631,  641,  643,  647,  653,  659, 
     855 661,  673,  677,  683,  691,  701,  709,  719,  727,  733, 
     856 739,  743,  751,  757,  761,  769,  773,  787,  797,  809, 
     857 811,  821,  823,  827,  829,  839,  853,  857,  859,  863, 
     858 877,  881,  883,  887,  907,  911,  919,  929,  937,  941, 
     859 947,  953,  967,  971,  977,  983,  991,  997, 1009};
     860        int num, den;
     861        num = 1;
     862        den = 1;
     863
     864        for (int i=0;i<r;i++){
     865                num = num*(n-i);
     866                den = den*(i+1);
     867                if (i>0) {
     868                        // Divide out common prime factors
     869                        for (int k=0;k<169;k++){ //169 is the length of the prime vector here
     870                                if ( i % primes[k] == 0) { // modulo
     871                                        num = num/primes[k];
     872                                        den = den/primes[k];
     873                                }
     874                        }
     875                }
     876        }
     877
     878        IssmDouble res;       
     879        return res = num/den;
     880}/*}}}*/
     881
     882
     883IssmDouble* postwidder_coef(int NTit){ /*{{{*/
     884        //Coefficients of the Post-Widder method through Saltzer summation for inverse Laplace transform:
     885        //The Mth iteration estimate will be: f(t)_M = sum_{k=1:2*M}(xi_[M,k] * f(s_k))
     886        //The method is based on equations (2), (6), (7) in:
     887        //Valko PP, Abate J. Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion. Computational Mathematics and Applications. (2004)
     888        //Note that the coefficients xi lack the factor s=k*log(2)/t.
     889        //That is because we are computing the heaviside response of the system rather than its impulse response,
     890        //and Laplace_Transform(Heaviside(t).*f(t)) = f(s)/s. So s cancels out in the sum for f(t)_M.
     891
     892        IssmDouble* xi=xNewZeroInit<IssmDouble>(2*NTit*NTit);
     893        int indxi;
     894        for (int M=1;M<NTit+1;M++){
     895                for (int k=1;k<2*M+1;k++){
     896                indxi=(M-1)*(2*NTit)+k-1;
     897                    for (int j=floor((k+1)/2);j<min(k,M)+1;j++){
     898                            xi[indxi]+=pow(j,M+1.0)/factorial(M)*n_C_r(M,j)*n_C_r(2*j,j)*n_C_r(j,k-j);
     899                    }
     900                    xi[indxi]*=pow(-1.0,k+M)/k;
     901                }
     902        }
     903        return xi;
     904}/*}}}*/
     905
     906void postwidder_transform(IssmDouble* Lovet, IssmDouble* Lovef,int d, int t, int NTit, IssmDouble* xi, FemModel* femmodel){ /*{{{*/
     907        //Computes Lovet for time step t and degree d from the PW coefficients xi and the corresponding 2*NTit frequency samples in Lovef
     908
     909        int nfreq, indxi, indf;
     910        femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
     911        int nt=nfreq/2/NTit;
     912
     913        indf=d*nfreq+t*2*NTit;
     914        IssmDouble LoveM[NTit];
     915
     916        for (int M=1;M<NTit+1;M++){
     917                LoveM[M-1]=0.0;
     918                for (int k=1;k<2*M+1;k++){
     919                        indxi=(M-1)*(2*NTit)+k-1;
     920                        LoveM[M-1]+=xi[indxi]*Lovef[indf+k-1];
     921                }
     922               
     923                //Make sure we are not getting into numerical instability
     924                //Diverging once: ok, we'll give that the benefit of the doubt, it can be an inflexion point in the convergence series
     925                //Diverging twice in a row: we are definitely propagating numerical error: revert to the last stable value and exit
     926                if (M>3){
     927                        if ( abs(LoveM[M-1]-LoveM[M-2]) > abs(LoveM[M-2]-LoveM[M-3]) &&
     928                             abs(LoveM[M-2]-LoveM[M-3]) > abs(LoveM[M-3]-LoveM[M-4]) ){
     929                                Lovet[d*nt+t]=LoveM[M-3];
     930                                return;
     931                        }
     932                }
     933        }
     934        Lovet[d*nt+t]=LoveM[NTit-1];
     935}/*}}}*/
     936
     937void love_freq_to_temporal(IssmDouble* LoveHt,IssmDouble* LoveLt,IssmDouble* LoveKt,IssmDouble* LoveHf,IssmDouble* LoveLf,IssmDouble* LoveKf, FemModel* femmodel){ /*{{{*/
     938        //Transforms all frequency-dependent love numbers into time-dependent love numbers
     939        int nfreq,sh_nmax,sh_nmin,indxi,indf, NTit;
     940        IssmDouble meanh,meanl,meank,dh,dl,dk;
     941
     942        femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
     943        femmodel->parameters->FindParam(&sh_nmax,LoveShNmaxEnum);
     944        femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum);
     945        femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
     946       
     947        int nt=nfreq/2/NTit;
     948        IssmDouble* xi=postwidder_coef(NTit);
     949
     950        for (int d=sh_nmin;d<sh_nmax+1;d++){
     951                for (int t=0;t<nt;t++){
     952                        postwidder_transform(LoveHt,LoveHf,d,t,NTit,xi,femmodel);
     953                        postwidder_transform(LoveLt,LoveLf,d,t,NTit,xi,femmodel);
     954                        postwidder_transform(LoveKt,LoveKf,d,t,NTit,xi,femmodel);
     955                }
     956        }
     957        xDelete<IssmDouble>(xi);
     958}/*}}}*/
     959
     960
     961
     962void love_core(FemModel* femmodel){ /*{{{*/
     963
     964        /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
     965        Matlitho* matlitho=NULL;
     966        for (Object* & object: femmodel->materials->objects){
     967                Material* material=xDynamicCast<Material*>(object);
     968                if(material->ObjectEnum()==MatlithoEnum){
     969                        matlitho=xDynamicCast<Matlitho*>(material);
     970                        break;
     971                }
     972        }
     973        _assert_(matlitho);
    17974
    18975        /*love parameters: */
     976        int         nfreq,nyi,starting_layer,nstep,forcing_type,dummy;
     977        int         sh_nmin,sh_nmax,kernel_index,deleted_layer_offset;
     978        bool        allow_layer_deletion,love_kernels, istemporal;
     979        bool        verbosemod = (int)VerboseModule();
     980        IssmDouble lovek, loveh, lovel, loveratio;
     981        IssmDouble  g0,r0,mu0, underflow_tol, omega;
    19982        IssmDouble *frequencies = NULL;
    20         int         nfreq,dummy;
    21         int         sh_nmin,sh_nmax;
    22         IssmDouble  g0,r0,mu0;
    23         bool        allow_layer_deletion;
    24         bool        love_kernels;
    25         int         forcing_type;
    26         bool        verbosemod = (int)VerboseModule();
    27983
    28984        /*parameters: */
    … …  
    35991
    36992        /*recover love number parameters: */
     993
    37994        femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
    38995        femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); _assert_(nfreq==dummy);
    … …  
    451002        femmodel->parameters->FindParam(&love_kernels,LoveKernelsEnum);
    461003        femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
    47 
    48         /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
    49         Matlitho* matlitho=NULL;
    50         for (Object* & object: femmodel->materials->objects){
    51                 Material* material=xDynamicCast<Material*>(object);
    52                 if(material->ObjectEnum()==MatlithoEnum){
    53                         matlitho=xDynamicCast<Matlitho*>(material);
    54                         break;
    55                 }
    56         }
    57         _assert_(matlitho);
     1004        femmodel->parameters->FindParam(&istemporal,LoveIsTemporalEnum);
     1005        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
    581006
    591007        /*Initialize three love matrices: geoid, vertical and horizontal displacement (real and imaginary parts) */
    60         IssmDouble*  LoveKr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    61         IssmDouble*  LoveHr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    62         IssmDouble*  LoveLr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    63         IssmDouble*  LoveKi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    64         IssmDouble*  LoveHi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    65         IssmDouble*  LoveLi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
     1008        IssmDouble*  LoveKf = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
     1009        IssmDouble*  LoveHf = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
     1010        IssmDouble*  LoveLf = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    661011
    671012        /*Initialize love kernels (real and imaginary parts): */
    … …  
    691014        IssmDouble*  LoveKernelsImag = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
    701015
     1016        love_init(femmodel,matlitho);
     1017        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
     1018        femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
     1019        IssmDouble* yi_prefactor=xNewZeroInit<IssmDouble>(6*6*nstep*matlitho->numlayers);
     1020        IssmDouble* yi_righthandside=xNewZeroInit<IssmDouble>(nyi);
     1021        IssmDouble* yi=xNewZeroInit<IssmDouble>(nyi*nyi);
     1022
     1023        //precalculate yi coefficients that do not depend on degree or frequency
     1024        fill_yi_prefactor(yi_prefactor, NULL, NULL,femmodel, matlitho);
     1025        for(int deg=sh_nmin;deg<sh_nmax+1;deg++){
     1026
     1027                //precalculate yi coefficients that depend on degree but not frequency
     1028                fill_yi_prefactor(yi_prefactor, &deg, NULL,femmodel, matlitho);
     1029
     1030                for (int fr=0;fr<nfreq;fr++){
     1031                        omega=2.0*PI*frequencies[fr]; //angular frequency
     1032
     1033                        //precalculate yi coefficients that depend on degree and frequency
     1034                        fill_yi_prefactor(yi_prefactor, &deg,&omega,femmodel, matlitho);
     1035
     1036                        yi_boundary_conditions(yi_righthandside,deg,femmodel,matlitho);
     1037
     1038                        build_yi_system(yi,deg,omega,yi_prefactor,femmodel,matlitho);
     1039
     1040                        solve_yi_system(&loveh,&lovel,&lovek, deg, omega, yi, yi_righthandside,femmodel, matlitho);
     1041
     1042                        LoveHf[deg*nfreq+fr]=loveh;
     1043                        LoveKf[deg*nfreq+fr]=lovek-1.0;
     1044                        LoveLf[deg*nfreq+fr]=lovel;
     1045
     1046                        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum); //synchronize nyi in case we deleted a layer
     1047                        deleted_layer_offset=(matlitho->numlayers+1)*6-nyi;// =6 per deleted layer
     1048                        kernel_index=deg*nfreq*(matlitho->numlayers+1)*6 +
     1049                                            fr*(matlitho->numlayers+1)*6 +
     1050                                                deleted_layer_offset;
     1051                        for (int i=0;i<nyi;i++){
     1052                                LoveKernelsReal[kernel_index+i]=yi_righthandside[i];
     1053                        }
     1054                }
     1055        }
     1056
     1057
     1058        xDelete<IssmDouble>(yi);
     1059        xDelete<IssmDouble>(yi_righthandside);
     1060
     1061
     1062        //Temporal love numbers
     1063        if (istemporal){
     1064                int NTit;
     1065                femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
     1066                int nt = nfreq/2/NTit;
     1067
     1068                IssmDouble* LoveHt=xNewZeroInit<IssmDouble>((sh_nmax+1)*nt);
     1069                IssmDouble* LoveLt=xNewZeroInit<IssmDouble>((sh_nmax+1)*nt);
     1070                IssmDouble* LoveKt=xNewZeroInit<IssmDouble>((sh_nmax+1)*nt);
     1071
     1072                love_freq_to_temporal(LoveHt,LoveLt,LoveKt,LoveHf,LoveLf,LoveKf,femmodel);
     1073
     1074                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKrEnum,LoveKt,sh_nmax+1,nt,0,0));
     1075                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHrEnum,LoveHt,sh_nmax+1,nt,0,0));
     1076                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLrEnum,LoveLt,sh_nmax+1,nt,0,0));
     1077
     1078                xDelete<IssmDouble>(LoveHt);
     1079                xDelete<IssmDouble>(LoveLt);
     1080                xDelete<IssmDouble>(LoveKt);
     1081        }
     1082
     1083        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKiEnum,LoveKf,sh_nmax+1,nfreq,0,0));
     1084        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHf,sh_nmax+1,nfreq,0,0));
     1085        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLf,sh_nmax+1,nfreq,0,0));
     1086
     1087
    711088        /*call the main module: */
    72         FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHi,LoveLr,LoveLi,LoveKernelsReal,LoveKernelsImag,  //output
    73                         nfreq,frequencies,sh_nmax,sh_nmin,g0,r0,mu0,allow_layer_deletion,forcing_type,verbosemod, //parameter inputs
    74                         matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu,
    75                         matlitho->burgers_viscosity, matlitho->burgers_mu, matlitho->density, matlitho->isburgers, matlitho->issolid //matlitho inputs
    76                         );
     1089        //if (false){
     1090        //FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHi,LoveLr,LoveLi,LoveKernelsReal,LoveKernelsImag,  //output
     1091        //              nfreq,frequencies,sh_nmax,sh_nmin,g0,r0,mu0,allow_layer_deletion,forcing_type,verbosemod, //parameter inputs
     1092        //              matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu,
     1093        //              matlitho->burgers_viscosity, matlitho->burgers_mu, matlitho->density, matlitho->rheologymodel, matlitho->issolid //matlitho inputs
     1094        //              );
     1095        //}
    771096
    781097        /*Add love matrices to results:*/
    79         femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKrEnum,LoveKr,sh_nmax+1,nfreq,0,0));
    80         femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHrEnum,LoveHr,sh_nmax+1,nfreq,0,0));
    81         femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLrEnum,LoveLr,sh_nmax+1,nfreq,0,0));
    82         femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKiEnum,LoveKi,sh_nmax+1,nfreq,0,0));
    83         femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHi,sh_nmax+1,nfreq,0,0));
    84         femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0));
    85         /*Only when love_kernels is set unity*/
     1098        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKrEnum,LoveKr,sh_nmax+1,nfreq,0,0));
     1099        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHrEnum,LoveHr,sh_nmax+1,nfreq,0,0));
     1100        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLrEnum,LoveLr,sh_nmax+1,nfreq,0,0));
     1101        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKiEnum,LoveKi,sh_nmax+1,nfreq,0,0));
     1102        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHi,sh_nmax+1,nfreq,0,0));
     1103        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0));
     1104
     1105        /*Only when love_kernels is on*/
    861106        if (love_kernels==1) {
    871107                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsRealEnum,LoveKernelsReal,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
    88                 femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
    89         }
    90 
    91         /*Free ressources:*/
     1108                //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
     1109        }
     1110
     1111        /*Free resources:*/
    921112        xDelete<IssmDouble>(frequencies);
    93         xDelete<IssmDouble>(LoveKr);
    94         xDelete<IssmDouble>(LoveHr);
    95         xDelete<IssmDouble>(LoveLr);
    96         xDelete<IssmDouble>(LoveKi);
    97         xDelete<IssmDouble>(LoveHi);
    98         xDelete<IssmDouble>(LoveLi);
     1113        xDelete<IssmDouble>(LoveKf);
     1114        xDelete<IssmDouble>(LoveHf);
     1115        xDelete<IssmDouble>(LoveLf);
     1116        //xDelete<IssmDouble>(LoveKi);
     1117        //xDelete<IssmDouble>(LoveHi);
     1118        //xDelete<IssmDouble>(LoveLi);
    991119        xDelete<IssmDouble>(LoveKernelsReal);
    1001120        xDelete<IssmDouble>(LoveKernelsImag);
    101 }
     1121} /*}}}*/
     1122
Note: See TracChangeset for help on using the changeset viewer.