Index: /issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp (revision 26241) +++ /issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp (revision 26242) @@ -42,4 +42,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.love.inner_core_boundary",LoveInnerCoreBoundaryEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.love.core_mantle_boundary",LoveCoreMantleBoundaryEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.love.complex_computation",LoveComplexComputationEnum)); }/*}}}*/ void LoveAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26241) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26242) @@ -66,6 +66,6 @@ /*Initialize solid earth motion and sea level: */ - InputUpdateFromConstantx(inputs,elements,0.,BedEastEnum); - InputUpdateFromConstantx(inputs,elements,0.,BedNorthEnum); + iomodel->ConstantToInput(inputs,elements,0.,BedEastEnum,P1Enum); + iomodel->ConstantToInput(inputs,elements,0.,BedNorthEnum,P1Enum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26241) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26242) @@ -7169,18 +7169,19 @@ IssmDouble sealevel[3]={0}; + int i,e,l,nbar; if(loads->sealevelloads){ - for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; - for (int e=0;esealevelloads[e]+loads->loads[e]); } - for(int l=0;lnbar[l]; - for (int e=0;enbar[l]; + for (e=0;esubloads[l][e]); } if(l==SLGEOM_OCEAN){ - for (int e=0;esubsealevelloads[e]); } @@ -7190,12 +7191,12 @@ } else{ //this is the initial convolution where only loads are provided - for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; - for (int e=0;eloads[e]); } - for(int l=0;lnbar[l]; - for (int e=0;enbar[l]; + for (e=0;esubloads[l][e]); } @@ -7206,10 +7207,10 @@ /*store values computed on vertices, but don't repeat the computation if another element already computed it!:*/ if(percpu){ - for(int i=0;ilids[this->vertices[i]->lid]==this->lid)sealevelout[this->vertices[i]->lid]=sealevel[i]; } } else{ - for(int i=0;i inline GenericExternalResult::GenericExternalResult(int in_id, int in_enum_type,IssmComplex* in_values, int in_M,int in_N,int in_step,IssmDouble in_time){/*{{{*/ + + id = in_id; + M = in_M; + N = in_N; + + EnumToStringx(&this->result_name,in_enum_type); + + step = in_step; + time = in_time; + + /*Copy result in values*/ + if(M*N){ + value=xNew(M*N); + xMemCpy(value,in_values,M*N); + } + else value=NULL; +} +/*}}}*/ +template <> inline void GenericExternalResult::WriteData(FILE* fid,bool io_gather){ /*{{{*/ + + int my_rank; + int type; + int rows,cols; + char *name = NULL; + IssmPDouble passiveDouble; + IssmDouble* reals=NULL; + IssmDouble* imags=NULL; + + /*recover my_rank:*/ + my_rank=IssmComm::GetRank(); + + if(io_gather){ + /*we are gathering the data on cpu 0, don't write on other cpus: */ + if(my_rank) return; + } + + /*First write enum: */ + int length=(strlen(this->result_name)+1)*sizeof(char); + fwrite(&length,sizeof(int),1,fid); + fwrite(this->result_name,length,1,fid); + + /*Now write time and step: */ + passiveDouble=reCast(time); + fwrite(&passiveDouble,sizeof(IssmPDouble),1,fid); + fwrite(&step,sizeof(int),1,fid); + + /*writing a IssmComplex array, type is 3:*/ + type=5; + fwrite(&type,sizeof(int),1,fid); + rows=this->M; + fwrite(&rows,sizeof(int),1,fid); + cols=this->N; + fwrite(&cols,sizeof(int),1,fid); + + /*write complex array into two real arrays:*/ + reals=xNew(cols*rows); + imags=xNew(cols*rows); + for(int i=0;i inline GenericExternalResult*>::GenericExternalResult(int in_id, int in_enum_type,Vector* in_values, int in_M,int in_N,int in_step,IssmDouble in_time){/*{{{*/ Index: /issm/trunk-jpl/src/c/cores/love_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/love_core.cpp (revision 26241) +++ /issm/trunk-jpl/src/c/cores/love_core.cpp (revision 26242) @@ -9,48 +9,161 @@ #include "../modules/modules.h" #include "../solutionsequences/solutionsequences.h" -#include - -void love_init(FemModel* femmodel, Matlitho* matlitho){/*{{{*/ - // initialize Planet_Mass(r) for efficient computation of gravity, value of surface gravity and inital size of the yi equation system - bool verbosemod = (int)VerboseModule(); - - int numlayers = matlitho->numlayers; - IssmDouble* r=matlitho->radius; - IssmDouble* EarthMass=xNewZeroInit(numlayers+1); - IssmDouble r1,r2,ro, GG; - femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); - - for (int i=0;idensity[i]; - if (i==0){ - EarthMass[i] = ro*pow(r2,3.0)*4.0*PI/3.0; - }else{ - r1=r[i]; - EarthMass[i] = EarthMass[i-1] + ro*(pow(r2,3.0)-pow(r1,3.0))*4.0*PI/3.0;; - } - } - - IssmDouble g0=EarthMass[numlayers-1]*GG/pow(r[numlayers],2.0); - femmodel->parameters->SetParam(g0,LoveG0Enum); - femmodel->parameters->SetParam(EarthMass,numlayers,LoveEarthMassEnum); - xDelete(EarthMass); - - int nyi=6*(numlayers+1); - int starting_layer=0; - femmodel->parameters->SetParam(nyi,LoveNYiEquationsEnum); - femmodel->parameters->SetParam(starting_layer,LoveStartingLayerEnum); -}/*}}}*/ -void GetEarthRheology(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega,FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ +#include "petscblaslapack.h" + +/*local definitions:*/ +class LoveVariables{ /*{{{*/ + + public: + IssmDouble g0; + IssmDouble r0; + IssmDouble* EarthMass; + int nyi; + int starting_layer; + + LoveVariables(){ /*{{{*/ + g0=0; + r0=0; + EarthMass=NULL; + nyi=0; + starting_layer=0; + } /*}}}*/ + LoveVariables(IssmDouble* EarthMassin,IssmDouble g0in,IssmDouble r0in,int nyiin,int starting_layerin){ /*{{{*/ + EarthMass=EarthMassin; + g0=g0in; + r0=r0in; + nyi=nyiin; + starting_layer=starting_layerin; + } /*}}}*/ + ~LoveVariables(){}; +}; /*}}}*/ + +/*self contained support routines used by cores below:*/ +template doubletype angular_frequency(IssmDouble frequency); //pure declaration +template <> IssmDouble angular_frequency(IssmDouble frequency){ /*{{{*/ + return 2.0*PI*frequency; +} /*}}}*/ +template <> IssmComplex angular_frequency(IssmDouble frequency){ /*{{{*/ + IssmComplex value=reCast(complex(0,1))*2.0*PI*reCast(frequency); + return value; +} /*}}}*/ +template void allgesv(int* pnyi, int* pnrhs, doubletype* yilocal, int* plda, int* ipiv, doubletype* rhslocal, int* pldb, int* pinfo); +template <> void allgesv(int* pnyi, int* pnrhs, IssmDouble* yilocal, int* plda, int* ipiv, IssmDouble* rhslocal, int* pldb, int* pinfo){ /*{{{*/ + dgesv_(pnyi, pnrhs, yilocal, plda, ipiv, rhslocal, pldb, pinfo); +} /*}}}*/ +template <> void allgesv(int* pnyi, int* pnrhs, IssmComplex* yilocal, int* plda, int* ipiv, IssmComplex* rhslocal, int* pldb, int* pinfo){ /*{{{*/ + _error_("zgesv_ not linked correctly yet! "); + //zgesv_(pnyi, pnrhs, yilocal, plda, ipiv, rhslocal, pldb, pinfo); +} /*}}}*/ +template doubletype factorial(int n){ /*{{{*/ + doubletype prod=1; + for (int i=2;i doubletype n_C_r(int n, int r){ /*{{{*/ + //n choose r + int primes[169] = + {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, + 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, + 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, + 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, + 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, + 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, + 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, + 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, + 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, + 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, + 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, + 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, + 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, + 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, + 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, + 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, + 947, 953, 967, 971, 977, 983, 991, 997, 1009}; + int num, den; + num = 1; + den = 1; + + for (int i=0;i0) { + // Divide out common prime factors + for (int k=0;k<169;k++){ //169 is the length of the prime vector here + if ( i % primes[k] == 0) { // modulo + num = num/primes[k]; + den = den/primes[k]; + } + } + } + } + + doubletype res; + return res = num/den; +}/*}}}*/ +template doubletype* postwidder_coef(int NTit){ /*{{{*/ + //Coefficients of the Post-Widder method through Saltzer summation for inverse Laplace transform: + //The Mth iteration estimate will be: f(t)_M = sum_{k=1:2*M}(xi_[M,k] * f(s_k)) + //The method is based on equations (2), (6), (7) in: + //Valko PP, Abate J. Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion. Computational Mathematics and Applications. (2004) + //Note that the coefficients xi lack the factor s=k*log(2)/t. + //That is because we are computing the heaviside response of the system rather than its impulse response, + //and Laplace_Transform(Heaviside(t).*f(t)) = f(s)/s. So s cancels out in the sum for f(t)_M. + + doubletype* xi=xNewZeroInit(2*NTit*NTit); + int indxi; + for (int M=1;M(M)*n_C_r(M,j)*n_C_r(2*j,j)*n_C_r(j,k-j); + } + xi[indxi]*=pow(-1.0,k+M)/k; + } + } + return xi; +}/*}}}*/ +template void postwidder_transform(doubletype* Lovet, doubletype* Lovef,int d, int t, int NTit, doubletype* xi, FemModel* femmodel){ /*{{{*/ + //Computes Lovet for time step t and degree d from the PW coefficients xi and the corresponding 2*NTit frequency samples in Lovef + + int nfreq, indxi, indf; + femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum); + int nt=nfreq/2/NTit; + + indf=d*nfreq+t*2*NTit; + doubletype LoveM[NTit]; + + for (int M=1;M3){ + if ( abs(LoveM[M-1]-LoveM[M-2]) > abs(LoveM[M-2]-LoveM[M-3]) && + abs(LoveM[M-2]-LoveM[M-3]) > abs(LoveM[M-3]-LoveM[M-4]) ){ + Lovet[d*nt+t]=LoveM[M-3]; + return; + } + } + } + Lovet[d*nt+t]=LoveM[NTit-1]; +}/*}}}*/ +template void GetEarthRheology(doubletype* pla, doubletype* pmu, int layer_index, doubletype omega, Matlitho* matlitho){ /*{{{*/ //returns lame parameters (material rigity) lambda and mu for the right frequency and layer - IssmDouble vi=matlitho->viscosity[layer_index]; - IssmDouble mu=matlitho->lame_mu[layer_index]; - IssmDouble la=matlitho->lame_lambda[layer_index]; - int rheo=matlitho->rheologymodel[layer_index]; + doubletype mu,la; + + IssmDouble vi=matlitho->viscosity[layer_index]; + IssmDouble mu0=matlitho->lame_mu[layer_index]; + IssmDouble la0=matlitho->lame_lambda[layer_index]; + int rheo=matlitho->rheologymodel[layer_index]; if (vi!=0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid - //Otherwise return the elastic value - IssmDouble ka=la + 2.0/3.0*mu; //Bulk modulus + //Otherwise return the elastic value + IssmDouble ka=la0 + 2.0/3.0*mu0; //Bulk modulus if (rheo==2){//EBM IssmDouble alpha=matlitho->ebm_alpha[layer_index]; @@ -58,5 +171,6 @@ IssmDouble taul=matlitho->ebm_taul[layer_index]; IssmDouble tauh=matlitho->ebm_tauh[layer_index]; - IssmDouble hf1,hf2,U1,U2; + IssmDouble hf1,hf2; + doubletype U1,U2; IssmDouble* a=xNew(2); IssmDouble* b=xNew(1); @@ -79,5 +193,5 @@ mu=mu/(1.0+ alpha*delta/(pow(tauh,alpha)-pow(taul,alpha))*(U1 + omega*U2) -mu/vi/omega); - la=ka-2.0/3.0*mu; + la=ka-2.0/3.0*mu; xDelete(a); @@ -88,23 +202,23 @@ IssmDouble mu2=matlitho->burgers_mu[layer_index]; - mu=mu*omega*(omega+mu2/vi2)/((omega+mu2/vi2)*(omega+mu/vi)+mu/vi2*omega); - la=ka-2.0/3.0*mu; - } + mu=mu0*omega*(omega+mu2/vi2)/((omega+mu2/vi2)*(omega+mu0/vi)+mu0/vi2*omega); + la=ka-2.0/3.0*mu0; + } else{//Maxwell - la = (la + mu*ka/vi/omega)/(1.0 + mu/vi/omega); - mu = mu/(1.0+mu/vi/omega); - } - } - - *pla=la; - *pmu=mu; + la = (la0 + mu0*ka/vi/omega)/(1.0 + mu0/vi/omega); + mu = mu0/(1.0+mu0/vi/omega); + } + } + + *pla=la; + *pmu=mu; } /*}}}*/ -IssmDouble GetGravity(IssmDouble r2, int layer_index, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ +IssmDouble GetGravity(IssmDouble r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars){ /*{{{*/ //computes gravity at radius r2 IssmDouble* EarthMass; IssmDouble g, GG; - femmodel->parameters->FindParam(&EarthMass,&matlitho->numlayers,LoveEarthMassEnum); + EarthMass=vars->EarthMass; femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); IssmDouble ro=matlitho->density[layer_index]; @@ -120,34 +234,36 @@ return g= GG*M/pow(r2,2.0); }/*}}}*/ -void fill_yi_prefactor(IssmDouble* yi_prefactor, int* pdeg, IssmDouble* pomega, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ +template void fill_yi_prefactor(doubletype* yi_prefactor, int* pdeg, doubletype* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ //precalculates partial derivative factors for function yi_derivatives - IssmDouble ra=matlitho->radius[matlitho->numlayers]; - IssmDouble g0,r0,mu0, GG; - int nstep,nindex, starting_layer; - femmodel->parameters->FindParam(&g0,LoveG0Enum); - femmodel->parameters->FindParam(&r0,LoveR0Enum); - femmodel->parameters->FindParam(&mu0,LoveMu0Enum); - femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); - femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); - femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum); - - IssmDouble frh,frhg0,fgr0,fgr,fn,rm0,rlm,flm; - IssmDouble xmin,xmax,x,dr; - IssmDouble g,ro, issolid; + IssmDouble ra=matlitho->radius[matlitho->numlayers]; + IssmDouble g0,r0,mu0, GG; + int nstep,nindex, starting_layer; + + femmodel->parameters->FindParam(&mu0,LoveMu0Enum); + femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); + femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + + g0=vars->g0; + r0=vars->r0; + starting_layer=vars->starting_layer; + + doubletype frh,frhg0,fgr0,fgr,fn,rm0,rlm,flm; + IssmDouble xmin,xmax,x,dr; + IssmDouble g,ro, issolid; if (pomega) { //frequency and degree dependent terms /*{{{*/ - IssmDouble la,mu; - IssmDouble* f=xNewZeroInit(12); + doubletype la,mu; + doubletype f[12]; int deg=*pdeg; - IssmDouble omega=*pomega; + doubletype omega=*pomega; fn=deg*(deg+1.0); for (int layer_index=starting_layer;layer_indexnumlayers;layer_index++){ - ro=matlitho->density[layer_index]; - issolid=matlitho->issolid[layer_index]; - if(issolid==1){ - GetEarthRheology(&la, &mu,layer_index,omega,femmodel,matlitho); - + ro=matlitho->density[layer_index]; + issolid=matlitho->issolid[layer_index]; + if(issolid==1){ + GetEarthRheology(&la, &mu,layer_index,omega,matlitho); + /*_______Expressions*/ flm=(la+2.0*mu); @@ -172,7 +288,7 @@ xmax=(matlitho->radius[layer_index+1])/ra; dr = (xmax -xmin)/nstep; - x=xmin; + x=xmin; for (int n=0;n(f); - /*}}}*/ + /*}}}*/ } else if (pdeg) { // degree dependent terms /*{{{*/ - int deg=*pdeg; - fn=(deg*(deg+1.0)); + int deg=*pdeg; + fn=(deg*(deg+1.0)); for (int layer_index=starting_layer;layer_indexnumlayers;layer_index++){ ro=matlitho->density[layer_index]; issolid=matlitho->issolid[layer_index]; - + /*_______Expressions*/ fgr=4.0*PI*GG*ro*ra; - + xmin=matlitho->radius[layer_index]/ra; xmax=(matlitho->radius[layer_index+1])/ra; dr = (xmax -xmin)/nstep; - x=xmin; + x=xmin; for (int n=0;nnumlayers;layer_index++){ ro=matlitho->density[layer_index]; issolid=matlitho->issolid[layer_index]; - + /*_______Expressions*/ frhg0=ro*g0*ra/mu0; @@ -236,118 +350,115 @@ xmax=(matlitho->radius[layer_index+1])/ra; dr = (xmax -xmin)/nstep; - x=xmin; - for (int n=0;n void yi_derivatives(doubletype* dydx, doubletype* y, int layer_index, int n, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ //computes yi derivatives at r=radius[layer_index]+ n/nstep*(radius[layer_index+1]-radius[layer_index]) IssmDouble issolid=matlitho->issolid[layer_index]; int iy,id,ny, nindex, nstep; - IssmDouble* dydx=xNewZeroInit(6); - IssmDouble* dy=xNewZeroInit(36); femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); /*{{{*/ /* For reference: - flm=(la+2.0*mu); - rlm=(3.0*la+2.0*mu)/(la+2.0*mu); - rm0=mu/mu0; - rg0=g/g0; - frh=ro*g*ra/mu0; - fgr=4.0*PI*GG*ro*ra/g0; - fn=(deg*(deg+1.0)); - - - if(issolid==1){ - ny = 6; - - dy[0*6+0]= (-2.0*la/flm)/x; - dy[0*6+1]= mu0/flm; - dy[0*6+2]= (la*fn/flm)/x; - dy[0*6+3]= 0.0; - dy[0*6+4]= 0.0; - dy[0*6+5]= 0.0; - - dy[1*6+0]= 4.0*(-frh+rm0*rlm/x)/x + ro*pow(omega,2.0)*ra/mu0; - dy[1*6+1]=(-4.0*mu/flm)/x; - dy[1*6+2]= fn*(frh-2.0*rm0*rlm/x)/x; - dy[1*6+3]= fn/x; - dy[1*6+4]= 0.0; - dy[1*6+5]= -frh/rg0; - - dy[2*6+0]= -1.0/x; - dy[2*6+1]= 0.0; - dy[2*6+2]= 1.0/x; - dy[2*6+3]= 1/rm0; - dy[2*6+4]= 0.0; - dy[2*6+5]= 0.0; - - dy[3*6+0]= (frh-2.0*rm0*rlm/x)/x; - dy[3*6+1]= ( -la/flm)/x; - 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; - dy[3*6+3]= -3.0/x; - dy[3*6+4]= -(frh/rg0)/x; - dy[3*6+5]= 0.0; - - dy[4*6+0]= fgr; - dy[4*6+1]= 0.0; - dy[4*6+2]= 0.0; - dy[4*6+3]= 0.0; - dy[4*6+4]= 0.0; - dy[4*6+5]= 1.0; - - dy[5*6+0]= 0.0; - dy[5*6+1]= 0.0; - dy[5*6+2]= -(fgr*fn)/x; - dy[5*6+3]= 0.0; - dy[5*6+4]= fn/(x*x); - dy[5*6+5]= -2.0/x; - - } else { - ny = 2; - - dy[0*6+0]= fgr/rg0; - dy[0*6+1]= 1.0; - dy[1*6+0]= (-4.0*(fgr/rg0)+fn/x)/x; - dy[1*6+1]= -2.0/x-fgr/rg0; - } - */ /*}}}*/ - - if(issolid==1){ - ny = 6; - } else { - ny = 2; - } - - for (id=0;id(dy); - return dydx; -}/*}}}*/ -void propagate_yi_euler(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ + } + } + return; +}/*}}}*/ +template void propagate_yi_euler(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ //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? //euler method @@ -355,17 +466,17 @@ femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); - IssmDouble* dydx=xNewZeroInit(6); + IssmDouble* dydx=xNewZeroInit(6); IssmDouble dr = (xmax -xmin)/nstep; - IssmDouble x=xmin; - for(int i = 0;i(dydx,y,layer_index, i,yi_prefactor,femmodel,matlitho); + for (int j=0;j<6;j++){ + y[j]+=dydx[j]*dr; + } + x = x + dr; + } xDelete(dydx); }/*}}}*/ -void propagate_yi_RK2(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ +template void propagate_yi_RK2(doubletype* y, IssmDouble xmin, IssmDouble xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ //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? //Implements Runge-Kutta 2nd order (midpoint method) @@ -373,24 +484,73 @@ femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); - IssmDouble* k1=xNewZeroInit(6); - IssmDouble* k2=xNewZeroInit(6); - IssmDouble* k3=xNewZeroInit(6); - IssmDouble* k4=xNewZeroInit(6); - IssmDouble* y1=xNewZeroInit(6); - IssmDouble* y2=xNewZeroInit(6); - IssmDouble* y3=xNewZeroInit(6); + doubletype k1[6]={0}; + doubletype k2[6]={0}; + doubletype k3[6]={0}; + doubletype k4[6]={0}; + doubletype y1[6]={0}; + doubletype y2[6]={0}; + doubletype y3[6]={0}; IssmDouble dr = (xmax -xmin)/nstep; - IssmDouble x=xmin; - for(int i = 0;i(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho); + for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;} + yi_derivatives(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + + for (int j=0;j<6;j++){ + y[j]+=k2[j]*2.0*dr; + } + x = x + 2.0*dr; + } +}/*}}}*/ +template void propagate_yi_RK4(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ + //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? + //Implements Runge-Kutta 4th order + int nstep; + femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + + IssmDouble* k1=xNewZeroInit(6); + IssmDouble* k2=xNewZeroInit(6); + IssmDouble* k3=xNewZeroInit(6); + IssmDouble* k4=xNewZeroInit(6); + IssmDouble* y1=xNewZeroInit(6); + IssmDouble* y2=xNewZeroInit(6); + IssmDouble* y3=xNewZeroInit(6); + + IssmDouble dr = (xmax -xmin)/nstep; + IssmDouble x=xmin; + for(int i = 0;i(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho); + for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;} + yi_derivatives(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;} + yi_derivatives(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;} + yi_derivatives(k4,y3,layer_index, 2*i+2,yi_prefactor,femmodel,matlitho); + + for (int j=0;j<6;j++){ + y[j]+=(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/3.0*dr; + } + x = x + 2.0*dr; + } + + //Last step: we don't know the derivative at xmax, so we will assume the values at xmax-dr + int i=nstep/2; + yi_derivatives(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho); for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;} - k2=yi_derivatives(y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); - - for (int j=0;j<6;j++){ - y[j]+=k2[j]*2.0*dr; - } - x = x + 2.0*dr; - } + yi_derivatives(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;} + yi_derivatives(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;} + yi_derivatives(k4,y3,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho); + + for (int j=0;j<6;j++){ + y[j]+=(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/3.0*dr; + } + + x = x + 2.0*dr; + xDelete(k1); xDelete(k2); @@ -401,60 +561,5 @@ xDelete(y3); }/*}}}*/ -void propagate_yi_RK4(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ - //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? - //Implements Runge-Kutta 4th order - int nstep; - femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); - - IssmDouble* k1=xNewZeroInit(6); - IssmDouble* k2=xNewZeroInit(6); - IssmDouble* k3=xNewZeroInit(6); - IssmDouble* k4=xNewZeroInit(6); - IssmDouble* y1=xNewZeroInit(6); - IssmDouble* y2=xNewZeroInit(6); - IssmDouble* y3=xNewZeroInit(6); - - IssmDouble dr = (xmax -xmin)/nstep; - IssmDouble x=xmin; - for(int i = 0;i(k1); - xDelete(k2); - xDelete(k3); - xDelete(k4); - xDelete(y1); - xDelete(y2); - xDelete(y3); -}/*}}}*/ -void Innersphere_boundaryconditions(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ +template void Innersphere_boundaryconditions(doubletype* yi, int layer_index, int deg, doubletype omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ //fills the boundary conditions at the bottom of layer[layer_index] in yi[0:2][0:5] @@ -463,19 +568,22 @@ IssmDouble g0,r0,mu0, GG; int nyi; - femmodel->parameters->FindParam(&g0,LoveG0Enum); - femmodel->parameters->FindParam(&r0,LoveR0Enum); - femmodel->parameters->FindParam(&mu0,LoveMu0Enum); - femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum); + + femmodel->parameters->FindParam(&mu0,LoveMu0Enum); femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); + + g0=vars->g0; + r0=vars->r0; + nyi=vars->nyi; + IssmDouble ro=matlitho->density[layer_index]; IssmDouble issolid=matlitho->issolid[layer_index]; - IssmDouble g=GetGravity(r,layer_index,femmodel,matlitho); - IssmDouble la,mu; + IssmDouble g=GetGravity(r,layer_index,femmodel,matlitho,vars); + doubletype la,mu; if (layer_index==0){ // 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 - GetEarthRheology(&la, &mu,layer_index,omega,femmodel,matlitho); + GetEarthRheology(&la, &mu,layer_index,omega,matlitho); } else { - GetEarthRheology(&la, &mu,layer_index-1,omega,femmodel,matlitho); + GetEarthRheology(&la, &mu,layer_index-1,omega,matlitho); } @@ -506,19 +614,22 @@ yi[5+nyi*1]=-cst/(r*g0); yi[5+nyi*2]=deg/(r*g0); -}/*}}}*/ -void build_yi_system(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho) { /*{{{*/ + +}/*}}}*/ +template void build_yi_system(doubletype* yi, int deg, doubletype omega, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars) { /*{{{*/ IssmDouble g0,r0,mu0,x,ro1, GG; int nyi,starting_layer, nstep; - femmodel->parameters->FindParam(&g0,LoveG0Enum); - femmodel->parameters->FindParam(&r0,LoveR0Enum); + + g0=vars->g0; + r0=vars->r0; + nyi=vars->nyi; + starting_layer=vars->starting_layer; + femmodel->parameters->FindParam(&mu0,LoveMu0Enum); - femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum); femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); - femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum); femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); - IssmDouble xmin,xmax,one,ro,g; - IssmDouble ra=matlitho->radius[matlitho->numlayers]; + IssmDouble xmin,xmax,one,ro,g; + IssmDouble ra=matlitho->radius[matlitho->numlayers]; for (int i=0;i<6*(matlitho->numlayers+1);i++){ @@ -529,7 +640,7 @@ int ny,is,ii,jj; - IssmDouble ystart[6]; + doubletype ystart[6]; for (int k=0;k<6;k++) ystart[k]=0.0; - + int ici = 0; // Index of current interface for (int i = starting_layer; inumlayers;i++){ @@ -538,23 +649,23 @@ xmax=(matlitho->radius[i+1])/ra; - + if (matlitho->issolid[i]){ - ny = 6; - is = 0; - one= 1.0; + ny = 6; + is = 0; + one= 1.0; } else { - ny = 2; - is = 4; - one= -1.0; - } - + ny = 2; + is = 4; + one= -1.0; + } + for (int j = 0;j(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho); //propagate_yi_RK4(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho); // Boundary Condition matrix - propagation part @@ -562,9 +673,9 @@ jj = 6*(ici+1)+j+is-3; for (int kk=0;kkissolid[i]){ @@ -576,5 +687,5 @@ } else { // Boundary Condition matrix - liquid regions ro1=matlitho->density[i]; - g=GetGravity(matlitho->radius[i], i, femmodel,matlitho); + g=GetGravity(matlitho->radius[i], i, femmodel,matlitho,vars); ii = 6*ici; yi[ii+nyi*(ii+3)] = -1.0; @@ -585,10 +696,10 @@ yi[(ii+4)+nyi*(ii+4+3)]=-1.0; yi[(ii+5)+nyi*(ii+5+3)]=-1.0; - g=GetGravity(matlitho->radius[i+1], i,femmodel,matlitho); + g=GetGravity(matlitho->radius[i+1], i,femmodel,matlitho,vars); ii = 6*(ici+1); yi[ii+nyi*(ii-1)]=-1.0; yi[ii+nyi*(ii+1)]=yi[(ii+4)+nyi*(ii+1)]*g0/g; // yi(17,14) solution integration 1 of z5 CMB yi[ii+nyi*(ii+2)]=yi[(ii+4)+nyi*(ii+2)]*g0/g; // yi(17,15) solution integration 2 of z5 CMB - // yi(13,..) y1 CMB + // yi(13,..) y1 CMB yi[(ii+1)+nyi*(ii-1)]=-ro1*g*ra/mu0; yi[(ii+2)+nyi*(ii)]=-1.0; @@ -599,8 +710,8 @@ - //-- Internal sphere: integration starts here rather than r=0 for numerical reasons - Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho); - - //-- Surface conditions + //-- Internal sphere: integration starts here rather than r=0 for numerical reasons + Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho,vars); + + //-- Surface conditions yi[(nyi-6)+nyi*(nyi-3)]=-1.0; yi[(nyi-4)+nyi*(nyi-2)]=-1.0; @@ -608,5 +719,5 @@ yi[(nyi-1)+nyi*(nyi-1)]=deg+1.0; - //-- Degree 1 special case + //-- Degree 1 special case if(deg==1){ for (int i=0;i void yi_boundary_conditions(doubletype* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars){ /*{{{*/ IssmDouble g0,r0,mu0,ra,rb,rc; - int nyi, forcing_type,icb,cmb; + int nyi, forcing_type,icb,cmb,starting_layer; IssmDouble* EarthMass; - femmodel->parameters->FindParam(&g0,LoveG0Enum); - femmodel->parameters->FindParam(&r0,LoveR0Enum); + + g0=vars->g0; + r0=vars->r0; + nyi=vars->nyi; + starting_layer=vars->starting_layer; + EarthMass=vars->EarthMass; + femmodel->parameters->FindParam(&mu0,LoveMu0Enum); - femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum); femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum); femmodel->parameters->FindParam(&icb,LoveInnerCoreBoundaryEnum); femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum); - femmodel->parameters->FindParam(&EarthMass,&matlitho->numlayers,LoveEarthMassEnum); - - // In Case of a Inner Core - Outer Core - Mantle planet and Boundary conditions on these 3 interfaces + + // In Case of a Inner Core - Outer Core - Mantle planet and Boundary conditions on these 3 interfaces ra=matlitho->radius[matlitho->numlayers]; rb=0; @@ -642,5 +757,4 @@ IssmDouble ro_mean=EarthMass[matlitho->numlayers-1]/(4.0/3.0*PI*pow(ra,3.0)); - xDelete(EarthMass); for (int i=0;i<(matlitho->numlayers+1)*6;i++) yi_righthandside[i]=0.0; @@ -648,55 +762,55 @@ switch (forcing_type) { - //-- forcings at the Inner Core Boundary - case 1: //'ICB --Volumetric Potential' - yi_righthandside[6*icb+5]=(deg)/(rc*g0); - yi_righthandside[6*icb+4]=1.0/(ra*g0); - break; - case 2: //'ICB --Pressure' - yi_righthandside[6*icb+1]=-ro_mean/mu0; - break; - case 3://'ICB --Loading' - yi_righthandside[6*icb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rc; - yi_righthandside[6*icb+5]= (2.0*deg+1.0)/(rc*g0); - break; - case 4://'ICB --Tangential Traction' - yi_righthandside[6*icb+3]= ro_mean/mu0; - break; - - //--forcings at the Core Mantle Boundary - case 5://'CMB --Volumetric Potential' - yi_righthandside[6*cmb+1]=-ro_mean/mu0*ra/rb; - yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0); - break; - case 6://'CMB --Pressure' - yi_righthandside[6*cmb+1]=-ro_mean/mu0; - break; - case 7://'CMB --Loading' - yi_righthandside[6*cmb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rb; - yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0); - break; - case 8://'CMB --Tangential Traction' - yi_righthandside[6*cmb+3]=-ro_mean/mu0; - break; - - //--forcings at the surface - case 9://'SURF--Volumetric Potential' - if (deg>1) yi_righthandside[nyi-1]=(2.0*deg+1.0)/(ra*g0); - break; - case 10://'SURF--Pressure' - yi_righthandside[nyi-5]=-ro_mean/mu0; - break; - case 11://'SURF--Loading' - yi_righthandside[nyi-5]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0); - if (deg>1) yi_righthandside[nyi-1]= (2.0*deg+1.0)/(ra*g0); - break; - case 12://'SURF--Tangential Traction' - yi_righthandside[nyi-3]= ro_mean/mu0; - break; - default: - _error_("love core error: forcing_type not supported yet"); - } -}/*}}}*/ -void solve_yi_system(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/ + //-- forcings at the Inner Core Boundary + case 1: //'ICB --Volumetric Potential' + yi_righthandside[6*icb+5]=(deg)/(rc*g0); + yi_righthandside[6*icb+4]=1.0/(ra*g0); + break; + case 2: //'ICB --Pressure' + yi_righthandside[6*icb+1]=-ro_mean/mu0; + break; + case 3://'ICB --Loading' + yi_righthandside[6*icb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rc; + yi_righthandside[6*icb+5]= (2.0*deg+1.0)/(rc*g0); + break; + case 4://'ICB --Tangential Traction' + yi_righthandside[6*icb+3]= ro_mean/mu0; + break; + + //--forcings at the Core Mantle Boundary + case 5://'CMB --Volumetric Potential' + yi_righthandside[6*cmb+1]=-ro_mean/mu0*ra/rb; + yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0); + break; + case 6://'CMB --Pressure' + yi_righthandside[6*cmb+1]=-ro_mean/mu0; + break; + case 7://'CMB --Loading' + yi_righthandside[6*cmb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rb; + yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0); + break; + case 8://'CMB --Tangential Traction' + yi_righthandside[6*cmb+3]=-ro_mean/mu0; + break; + + //--forcings at the surface + case 9://'SURF--Volumetric Potential' + if (deg>1) yi_righthandside[nyi-1]=(2.0*deg+1.0)/(ra*g0); + break; + case 10://'SURF--Pressure' + yi_righthandside[nyi-5]=-ro_mean/mu0; + break; + case 11://'SURF--Loading' + yi_righthandside[nyi-5]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0); + if (deg>1) yi_righthandside[nyi-1]= (2.0*deg+1.0)/(ra*g0); + break; + case 12://'SURF--Tangential Traction' + yi_righthandside[nyi-3]= ro_mean/mu0; + break; + default: + _error_("love core error: forcing_type not supported yet"); + } +}/*}}}*/ +template void solve_yi_system(doubletype* loveh, doubletype* lovel, doubletype* lovek, int deg, doubletype omega, doubletype* yi, doubletype* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/ IssmDouble g0,r0,mu0,loveratio,underflow_tol; @@ -704,13 +818,17 @@ int nyi,starting_layer, dummy; bool allow_layer_deletion; - femmodel->parameters->FindParam(&g0,LoveG0Enum); - femmodel->parameters->FindParam(&r0,LoveR0Enum); + IssmDouble* EarthMass=NULL; + + g0=vars->g0; + r0=vars->r0; + nyi=vars->nyi; + starting_layer=vars->starting_layer; + EarthMass=vars->EarthMass; + femmodel->parameters->FindParam(&mu0,LoveMu0Enum); - femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum); - femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum); femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum); femmodel->parameters->FindParam(&underflow_tol,LoveUnderflowTolEnum); femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); - IssmDouble ra=matlitho->radius[matlitho->numlayers]; + IssmDouble ra=matlitho->radius[matlitho->numlayers]; bool exit=false; int lda,ldb; @@ -719,6 +837,6 @@ lda=nyi; ldb=nyi; - IssmDouble* yilocal=xNew(nyi*nyi); // we will need to redeclare these inside here as nyi changes - IssmDouble* rhslocal=xNew(nyi); + doubletype* yilocal=xNew(nyi*nyi); // we will need to redeclare these inside here as nyi changes + doubletype* rhslocal=xNew(nyi); //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 @@ -734,8 +852,9 @@ int nrhs=1; // number of right hand size columns - dgesv_(&nyi, &nrhs, yilocal, &lda, ipiv, rhslocal, &ldb, &info); + allgesv(&nyi, &nrhs, yilocal, &lda, ipiv, rhslocal, &ldb, &info); + xDelete(ipiv); - - if(VerboseSolution() && info!=0){ + + if(VerboseModule() && info!=0){ _printf0_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system"); printf("%s %s %s %s\n", "i","j","yi[i+nyi*j]","rhs[i]"); @@ -751,11 +870,11 @@ *lovek = rhslocal[nyi-1]*ra*g0; - IssmDouble loveh1 = rhslocal[3]; - IssmDouble lovel1 = rhslocal[5]; - IssmDouble lovek1 = rhslocal[7] - pow(matlitho->radius[starting_layer]/ra,deg)/(g0*ra); - - IssmDouble loveh1s = rhslocal[nyi-3]; - IssmDouble lovel1s = rhslocal[nyi-2]; - IssmDouble lovek1s = rhslocal[nyi-1] - 1.0/(g0*ra); + doubletype loveh1 = rhslocal[3]; + doubletype lovel1 = rhslocal[5]; + doubletype lovek1 = rhslocal[7] - pow(matlitho->radius[starting_layer]/ra,deg)/(g0*ra); + + doubletype loveh1s = rhslocal[nyi-3]; + doubletype lovel1s = rhslocal[nyi-2]; + doubletype lovek1s = rhslocal[nyi-1] - 1.0/(g0*ra); loveratio = abs(loveh1/loveh1s); //ratio of center to surface love numbers, determines if we should remove layers @@ -763,6 +882,9 @@ if (abs(lovek1/lovek1s) < loveratio) loveratio = abs(lovek1/lovek1s); - 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. + if (!allow_layer_deletion || nyi<=12 || omega!=angular_frequency(frequencies[0])){ goto save_results; + /*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.*/ } @@ -775,7 +897,6 @@ nyi-=6; starting_layer+=1; - femmodel->parameters->SetParam(nyi,LoveNYiEquationsEnum); - femmodel->parameters->SetParam(starting_layer,LoveStartingLayerEnum); - + vars->nyi=nyi; + vars->starting_layer=starting_layer; for (int i=0;i(yi, starting_layer, deg, omega, femmodel, matlitho,vars); //we move the first interface to the new starting layer. yi[0:2,0:5] will be different } else { //we are ready to save the outputs and break the main loop - save_results: +save_results: for (int i=0;i(yilocal); - xDelete(rhslocal); - } + } + xDelete(yilocal); + xDelete(rhslocal); + } xDelete(frequencies); -}/*}}}*/ -IssmDouble factorial(int n){ /*{{{*/ - IssmDouble prod=1; - for (int i=2;i0) { - // Divide out common prime factors - for (int k=0;k<169;k++){ //169 is the length of the prime vector here - if ( i % primes[k] == 0) { // modulo - num = num/primes[k]; - den = den/primes[k]; - } - } - } - } - - IssmDouble res; - return res = num/den; -}/*}}}*/ -IssmDouble* postwidder_coef(int NTit){ /*{{{*/ - //Coefficients of the Post-Widder method through Saltzer summation for inverse Laplace transform: - //The Mth iteration estimate will be: f(t)_M = sum_{k=1:2*M}(xi_[M,k] * f(s_k)) - //The method is based on equations (2), (6), (7) in: - //Valko PP, Abate J. Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion. Computational Mathematics and Applications. (2004) - //Note that the coefficients xi lack the factor s=k*log(2)/t. - //That is because we are computing the heaviside response of the system rather than its impulse response, - //and Laplace_Transform(Heaviside(t).*f(t)) = f(s)/s. So s cancels out in the sum for f(t)_M. - - IssmDouble* xi=xNewZeroInit(2*NTit*NTit); - int indxi; - for (int M=1;Mparameters->FindParam(&nfreq,LoveNfreqEnum); - int nt=nfreq/2/NTit; - - indf=d*nfreq+t*2*NTit; - IssmDouble LoveM[NTit]; - - for (int M=1;M3){ - if ( abs(LoveM[M-1]-LoveM[M-2]) > abs(LoveM[M-2]-LoveM[M-3]) && - abs(LoveM[M-2]-LoveM[M-3]) > abs(LoveM[M-3]-LoveM[M-4]) ){ - Lovet[d*nt+t]=LoveM[M-3]; - return; - } - } - } - Lovet[d*nt+t]=LoveM[NTit-1]; -}/*}}}*/ -void love_freq_to_temporal(IssmDouble* LoveHt,IssmDouble* LoveLt,IssmDouble* LoveKt,IssmDouble* LoveHf,IssmDouble* LoveLf,IssmDouble* LoveKf, FemModel* femmodel){ /*{{{*/ + +}/*}}}*/ +template void love_freq_to_temporal(doubletype* LoveHt,doubletype* LoveLt,doubletype* LoveKt,doubletype* LoveHf,doubletype* LoveLf,doubletype* LoveKf, FemModel* femmodel){ /*{{{*/ //Transforms all frequency-dependent love numbers into time-dependent love numbers int nfreq,sh_nmax,sh_nmin,indxi,indf, NTit; - IssmDouble meanh,meanl,meank,dh,dl,dk; femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum); @@ -924,21 +947,82 @@ femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum); femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum); - + int nt=nfreq/2/NTit; - IssmDouble* xi=postwidder_coef(NTit); + doubletype* xi=postwidder_coef(NTit); for (int d=sh_nmin;d(xi); -}/*}}}*/ -void love_core(FemModel* femmodel){ /*{{{*/ + postwidder_transform(LoveHt,LoveHf,d,t,NTit,xi,femmodel); + postwidder_transform(LoveLt,LoveLf,d,t,NTit,xi,femmodel); + postwidder_transform(LoveKt,LoveKf,d,t,NTit,xi,femmodel); + } + } + xDelete(xi); +}/*}}}*/ + +/*templated cores:*/ +LoveVariables* love_init(FemModel* femmodel, Matlitho* matlitho){/*{{{*/ + + /*initialize Planet_Mass(r) for efficient computation of gravity, value of surface gravity and inital size of the yi equation system*/ + + bool verbosemod = (int)VerboseModule(); + int numlayers = matlitho->numlayers; + IssmDouble* r = matlitho->radius; + IssmDouble r1,r2,ro, GG; + + /*outputs:*/ + IssmDouble* EarthMass=NULL; + IssmDouble g0,r0; + int nyi,starting_layer; + + femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum); + EarthMass=xNewZeroInit(numlayers+1); + + for (int i=0;idensity[i]; + if (i==0){ + EarthMass[i] = ro*pow(r2,3.0)*4.0*PI/3.0; + }else{ + r1=r[i]; + EarthMass[i] = EarthMass[i-1] + ro*(pow(r2,3.0)-pow(r1,3.0))*4.0*PI/3.0;; + } + } + + g0=EarthMass[numlayers-1]*GG/pow(r[numlayers],2.0); + r0=r[numlayers]; + nyi=6*(numlayers+1); + starting_layer=0; + + return new LoveVariables(EarthMass,g0,r0,nyi,starting_layer); + +} /*}}}*/ +template void love_core_template(FemModel* femmodel){ /*{{{*/ + + Matlitho* matlitho=NULL; + int nfreq,nyi,starting_layer,nstep,forcing_type,dummy; + int sh_nmin,sh_nmax,kernel_index,deleted_layer_offset; + bool allow_layer_deletion,love_kernels, istemporal; + bool verbosemod = (int)VerboseModule(); + IssmDouble g0,r0,mu0, underflow_tol; + IssmDouble *frequencies = NULL; + bool save_results; + bool complex_computation; + + doubletype omega; + doubletype lovek, loveh, lovel, loveratio; + + doubletype* LoveKf = NULL; + doubletype* LoveHf = NULL; + doubletype* LoveLf = NULL; + doubletype* LoveKernels= NULL; + LoveVariables* vars=NULL; + doubletype* yi_prefactor=NULL; + doubletype* yi_righthandside=NULL; + doubletype* yi=NULL; + + if(VerboseSolution()) _printf0_(" computing LOVE numbers\n"); /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/ - Matlitho* matlitho=NULL; for (Object* & object: femmodel->materials->objects){ Material* material=xDynamicCast(object); @@ -950,23 +1034,6 @@ _assert_(matlitho); - /*love parameters: */ - int nfreq,nyi,starting_layer,nstep,forcing_type,dummy; - int sh_nmin,sh_nmax,kernel_index,deleted_layer_offset; - bool allow_layer_deletion,love_kernels, istemporal; - bool verbosemod = (int)VerboseModule(); - IssmDouble lovek, loveh, lovel, loveratio; - IssmDouble g0,r0,mu0, underflow_tol, omega; - IssmDouble *frequencies = NULL; - - /*parameters: */ - bool save_results; - - if(VerboseSolution()) _printf0_(" computing LOVE numbers\n"); - - /*Recover some parameters: */ + /*recover parameters: */ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - - /*recover love number parameters: */ - femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum); femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); _assert_(nfreq==dummy); @@ -981,39 +1048,39 @@ femmodel->parameters->FindParam(&istemporal,LoveIsTemporalEnum); femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); + femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum); + /*Initialize three love matrices: geoid, vertical and horizontal displacement (real and imaginary parts) */ - IssmDouble* LoveKf = xNewZeroInit(nfreq*(sh_nmax+1)); - IssmDouble* LoveHf = xNewZeroInit(nfreq*(sh_nmax+1)); - IssmDouble* LoveLf = xNewZeroInit(nfreq*(sh_nmax+1)); + LoveKf = xNewZeroInit(nfreq*(sh_nmax+1)); + LoveHf = xNewZeroInit(nfreq*(sh_nmax+1)); + LoveLf = xNewZeroInit(nfreq*(sh_nmax+1)); /*Initialize love kernels (real and imaginary parts): */ - IssmDouble* LoveKernelsReal = xNewZeroInit(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6); - IssmDouble* LoveKernelsImag = xNewZeroInit(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6); - - love_init(femmodel,matlitho); - femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum); - femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum); - IssmDouble* yi_prefactor=xNewZeroInit(6*6*nstep*matlitho->numlayers); - IssmDouble* yi_righthandside=xNewZeroInit(nyi); - IssmDouble* yi=xNewZeroInit(nyi*nyi); + LoveKernels= xNewZeroInit(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6); + + vars=love_init(femmodel,matlitho); + + yi_prefactor=xNewZeroInit(6*6*nstep*matlitho->numlayers); + yi_righthandside=xNewZeroInit(vars->nyi); + yi=xNewZeroInit(vars->nyi*vars->nyi); //precalculate yi coefficients that do not depend on degree or frequency - fill_yi_prefactor(yi_prefactor, NULL, NULL,femmodel, matlitho); + fill_yi_prefactor(yi_prefactor, NULL, NULL,femmodel, matlitho,vars); for(int deg=sh_nmin;deg(yi_prefactor, °, NULL,femmodel, matlitho,vars); for (int fr=0;fr(frequencies[fr]); //precalculate yi coefficients that depend on degree and frequency - fill_yi_prefactor(yi_prefactor, °,&omega,femmodel, matlitho); - - yi_boundary_conditions(yi_righthandside,deg,femmodel,matlitho); - - build_yi_system(yi,deg,omega,yi_prefactor,femmodel,matlitho); - - solve_yi_system(&loveh,&lovel,&lovek, deg, omega, yi, yi_righthandside,femmodel, matlitho); + fill_yi_prefactor(yi_prefactor, °,&omega,femmodel, matlitho,vars); + + yi_boundary_conditions(yi_righthandside,deg,femmodel,matlitho,vars); + + build_yi_system(yi,deg,omega,yi_prefactor,femmodel,matlitho,vars); + + solve_yi_system(&loveh,&lovel,&lovek, deg, omega, yi, yi_righthandside,femmodel, matlitho,vars); LoveHf[deg*nfreq+fr]=loveh; @@ -1021,44 +1088,41 @@ LoveLf[deg*nfreq+fr]=lovel; - femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum); //synchronize nyi in case we deleted a layer - deleted_layer_offset=(matlitho->numlayers+1)*6-nyi;// =6 per deleted layer - kernel_index=deg*nfreq*(matlitho->numlayers+1)*6 + - fr*(matlitho->numlayers+1)*6 + - deleted_layer_offset; - for (int i=0;i(yi); - xDelete(yi_righthandside); + deleted_layer_offset=(matlitho->numlayers+1)*6-vars->nyi;// =6 per deleted layer + kernel_index=deg*nfreq*(matlitho->numlayers+1)*6 + fr*(matlitho->numlayers+1)*6 + deleted_layer_offset; + for (int i=0;inyi;i++){ + LoveKernels[kernel_index+i]=yi_righthandside[i]; + } + } + } + + + xDelete(yi); + xDelete(yi_righthandside); //Temporal love numbers - if (istemporal){ + if (istemporal && !complex_computation){ int NTit; femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum); int nt = nfreq/2/NTit; - IssmDouble* LoveHt=xNewZeroInit((sh_nmax+1)*nt); - IssmDouble* LoveLt=xNewZeroInit((sh_nmax+1)*nt); - IssmDouble* LoveKt=xNewZeroInit((sh_nmax+1)*nt); - - love_freq_to_temporal(LoveHt,LoveLt,LoveKt,LoveHf,LoveLf,LoveKf,femmodel); - - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKrEnum,LoveKt,sh_nmax+1,nt,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHrEnum,LoveHt,sh_nmax+1,nt,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLrEnum,LoveLt,sh_nmax+1,nt,0,0)); - - xDelete(LoveHt); - xDelete(LoveLt); - xDelete(LoveKt); - } - - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKiEnum,LoveKf,sh_nmax+1,nfreq,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHiEnum,LoveHf,sh_nmax+1,nfreq,0,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLiEnum,LoveLf,sh_nmax+1,nfreq,0,0)); + doubletype* LoveHt=xNewZeroInit((sh_nmax+1)*nt); + doubletype* LoveLt=xNewZeroInit((sh_nmax+1)*nt); + doubletype* LoveKt=xNewZeroInit((sh_nmax+1)*nt); + + love_freq_to_temporal(LoveHt,LoveLt,LoveKt,LoveHf,LoveLf,LoveKf,femmodel); + + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKrEnum,LoveKt,sh_nmax+1,nt,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHrEnum,LoveHt,sh_nmax+1,nt,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLrEnum,LoveLt,sh_nmax+1,nt,0,0)); + + xDelete(LoveHt); + xDelete(LoveLt); + xDelete(LoveKt); + } + + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKiEnum,LoveKf,sh_nmax+1,nfreq,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveHiEnum,LoveHf,sh_nmax+1,nfreq,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveLiEnum,LoveLf,sh_nmax+1,nfreq,0,0)); @@ -1082,5 +1146,5 @@ /*Only when love_kernels is on*/ if (love_kernels==1) { - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKernelsRealEnum,LoveKernelsReal,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernels,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0)); //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0)); } @@ -1088,11 +1152,45 @@ /*Free resources:*/ xDelete(frequencies); - xDelete(LoveKf); - xDelete(LoveHf); - xDelete(LoveLf); - //xDelete(LoveKi); - //xDelete(LoveHi); + xDelete(LoveKf); + xDelete(LoveHf); + xDelete(LoveLf); //xDelete(LoveLi); - xDelete(LoveKernelsReal); - xDelete(LoveKernelsImag); + xDelete(LoveKernels); } /*}}}*/ + +/*cores and template instantiations:*/ +/*template instantiations :{{{*/ +template void love_core_template(FemModel* femmodel); +template void love_core_template(FemModel* femmodel); +template void fill_yi_prefactor(IssmDouble* yi_prefactor, int* pdeg, IssmDouble* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void fill_yi_prefactor(IssmComplex* yi_prefactor, int* pdeg, IssmComplex* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void GetEarthRheology(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega, Matlitho* matlitho); +template void GetEarthRheology(IssmComplex* pla, IssmComplex* pmu, int layer_index, IssmComplex omega, Matlitho* matlitho); +template void yi_boundary_conditions(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars); +template void yi_boundary_conditions(IssmComplex* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars); +template void yi_derivatives(IssmDouble* dydx, IssmDouble* y, int layer_index, int n, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho); +template void yi_derivatives(IssmComplex* dydx, IssmComplex* y, int layer_index, int n, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho); +template void propagate_yi_RK2(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho); +template void propagate_yi_RK2(IssmComplex* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho); +template void Innersphere_boundaryconditions(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void Innersphere_boundaryconditions(IssmComplex* yi, int layer_index, int deg, IssmComplex omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void build_yi_system(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars); +template void build_yi_system(IssmComplex* yi, int deg, IssmComplex omega, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars); +template void solve_yi_system(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template void solve_yi_system(IssmComplex* loveh, IssmComplex* lovel, IssmComplex* lovek, int deg, IssmComplex omega, IssmComplex* yi, IssmComplex* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars); +template IssmDouble factorial(int n); +template IssmDouble* postwidder_coef(int NTit); +template IssmDouble n_C_r(int n, int r); +template void postwidder_transform(IssmDouble* Lovet, IssmDouble* Lovef,int d, int t, int NTit, IssmDouble* xi, FemModel* femmodel); + +/*}}}*/ +void love_core(FemModel* femmodel){ /*{{{*/ + + bool complex_computation; + + femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum); + + if(complex_computation) love_core_template(femmodel); + else love_core_template(femmodel); + +} /*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26241) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26242) @@ -256,13 +256,24 @@ syn keyword cConstant LockFileNameEnum syn keyword cConstant LoveAllowLayerDeletionEnum +syn keyword cConstant LoveCoreMantleBoundaryEnum +syn keyword cConstant LoveEarthMassEnum syn keyword cConstant LoveForcingTypeEnum syn keyword cConstant LoveFrequenciesEnum +syn keyword cConstant LoveIsTemporalEnum syn keyword cConstant LoveG0Enum +syn keyword cConstant LoveGravitationalConstantEnum +syn keyword cConstant LoveInnerCoreBoundaryEnum +syn keyword cConstant LoveComplexComputationEnum +syn keyword cConstant LoveIntStepsPerLayerEnum syn keyword cConstant LoveKernelsEnum syn keyword cConstant LoveMu0Enum syn keyword cConstant LoveNfreqEnum +syn keyword cConstant LoveNTemporalIterationsEnum +syn keyword cConstant LoveNYiEquationsEnum syn keyword cConstant LoveR0Enum syn keyword cConstant LoveShNmaxEnum syn keyword cConstant LoveShNminEnum +syn keyword cConstant LoveStartingLayerEnum +syn keyword cConstant LoveUnderflowTolEnum syn keyword cConstant MassFluxSegmentsEnum syn keyword cConstant MassFluxSegmentsPresentEnum @@ -783,4 +794,13 @@ syn keyword cConstant SealevelGrotm2Enum syn keyword cConstant SealevelGrotm3Enum +syn keyword cConstant SealevelGUrotm1Enum +syn keyword cConstant SealevelGUrotm2Enum +syn keyword cConstant SealevelGUrotm3Enum +syn keyword cConstant SealevelGNrotm1Enum +syn keyword cConstant SealevelGNrotm2Enum +syn keyword cConstant SealevelGNrotm3Enum +syn keyword cConstant SealevelGErotm1Enum +syn keyword cConstant SealevelGErotm2Enum +syn keyword cConstant SealevelGErotm3Enum syn keyword cConstant SealevelRSLBarystaticEnum syn keyword cConstant SealevelRSLRateEnum @@ -1528,4 +1548,5 @@ syn keyword cType GenericParam syn keyword cType GiaDeflectionCoreArgs +syn keyword cType GrdLoads syn keyword cType Hook syn keyword cType Input Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26241) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26242) @@ -258,4 +258,5 @@ LoveGravitationalConstantEnum, LoveInnerCoreBoundaryEnum, + LoveComplexComputationEnum, LoveIntStepsPerLayerEnum, LoveKernelsEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26241) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26242) @@ -258,13 +258,24 @@ case LockFileNameEnum : return "LockFileName"; case LoveAllowLayerDeletionEnum : return "LoveAllowLayerDeletion"; + case LoveCoreMantleBoundaryEnum : return "LoveCoreMantleBoundary"; + case LoveEarthMassEnum : return "LoveEarthMass"; case LoveForcingTypeEnum : return "LoveForcingType"; case LoveFrequenciesEnum : return "LoveFrequencies"; + case LoveIsTemporalEnum : return "LoveIsTemporal"; case LoveG0Enum : return "LoveG0"; + case LoveGravitationalConstantEnum : return "LoveGravitationalConstant"; + case LoveInnerCoreBoundaryEnum : return "LoveInnerCoreBoundary"; + case LoveComplexComputationEnum : return "LoveComplexComputation"; + case LoveIntStepsPerLayerEnum : return "LoveIntStepsPerLayer"; case LoveKernelsEnum : return "LoveKernels"; case LoveMu0Enum : return "LoveMu0"; case LoveNfreqEnum : return "LoveNfreq"; + case LoveNTemporalIterationsEnum : return "LoveNTemporalIterations"; + case LoveNYiEquationsEnum : return "LoveNYiEquations"; case LoveR0Enum : return "LoveR0"; case LoveShNmaxEnum : return "LoveShNmax"; case LoveShNminEnum : return "LoveShNmin"; + case LoveStartingLayerEnum : return "LoveStartingLayer"; + case LoveUnderflowTolEnum : return "LoveUnderflowTol"; case MassFluxSegmentsEnum : return "MassFluxSegments"; case MassFluxSegmentsPresentEnum : return "MassFluxSegmentsPresent"; @@ -785,4 +796,13 @@ case SealevelGrotm2Enum : return "SealevelGrotm2"; case SealevelGrotm3Enum : return "SealevelGrotm3"; + case SealevelGUrotm1Enum : return "SealevelGUrotm1"; + case SealevelGUrotm2Enum : return "SealevelGUrotm2"; + case SealevelGUrotm3Enum : return "SealevelGUrotm3"; + case SealevelGNrotm1Enum : return "SealevelGNrotm1"; + case SealevelGNrotm2Enum : return "SealevelGNrotm2"; + case SealevelGNrotm3Enum : return "SealevelGNrotm3"; + case SealevelGErotm1Enum : return "SealevelGErotm1"; + case SealevelGErotm2Enum : return "SealevelGErotm2"; + case SealevelGErotm3Enum : return "SealevelGErotm3"; case SealevelRSLBarystaticEnum : return "SealevelRSLBarystatic"; case SealevelRSLRateEnum : return "SealevelRSLRate"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26241) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26242) @@ -264,13 +264,24 @@ if(stage==3){ if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum; + else if (strcmp(name,"LoveCoreMantleBoundary")==0) return LoveCoreMantleBoundaryEnum; + else if (strcmp(name,"LoveEarthMass")==0) return LoveEarthMassEnum; else if (strcmp(name,"LoveForcingType")==0) return LoveForcingTypeEnum; else if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum; + else if (strcmp(name,"LoveIsTemporal")==0) return LoveIsTemporalEnum; else if (strcmp(name,"LoveG0")==0) return LoveG0Enum; + else if (strcmp(name,"LoveGravitationalConstant")==0) return LoveGravitationalConstantEnum; + else if (strcmp(name,"LoveInnerCoreBoundary")==0) return LoveInnerCoreBoundaryEnum; + else if (strcmp(name,"LoveComplexComputation")==0) return LoveComplexComputationEnum; + else if (strcmp(name,"LoveIntStepsPerLayer")==0) return LoveIntStepsPerLayerEnum; else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum; else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum; else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; + else if (strcmp(name,"LoveNTemporalIterations")==0) return LoveNTemporalIterationsEnum; + else if (strcmp(name,"LoveNYiEquations")==0) return LoveNYiEquationsEnum; else if (strcmp(name,"LoveR0")==0) return LoveR0Enum; else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum; + else if (strcmp(name,"LoveStartingLayer")==0) return LoveStartingLayerEnum; + else if (strcmp(name,"LoveUnderflowTol")==0) return LoveUnderflowTolEnum; else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum; else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum; @@ -372,5 +383,8 @@ else if (strcmp(name,"SealevelchangeGRigid")==0) return SealevelchangeGRigidEnum; else if (strcmp(name,"SealevelchangeGElastic")==0) return SealevelchangeGElasticEnum; - else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum; else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum; else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; @@ -383,8 +397,5 @@ else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum; else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum; + else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum; else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum; else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum; @@ -495,5 +506,8 @@ else if (strcmp(name,"TimesteppingType")==0) return TimesteppingTypeEnum; else if (strcmp(name,"ToMITgcmComm")==0) return ToMITgcmCommEnum; - else if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum; else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum; else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum; @@ -506,8 +520,5 @@ else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; + else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; @@ -618,5 +629,8 @@ else if (strcmp(name,"DistanceToCalvingfront")==0) return DistanceToCalvingfrontEnum; else if (strcmp(name,"DistanceToGroundingline")==0) return DistanceToGroundinglineEnum; - else if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum; else if (strcmp(name,"Domain2Dvertical")==0) return Domain2DverticalEnum; else if (strcmp(name,"Domain3D")==0) return Domain3DEnum; @@ -629,8 +643,5 @@ else if (strcmp(name,"EffectivePressureTransient")==0) return EffectivePressureTransientEnum; else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; + else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; else if (strcmp(name,"EplHead")==0) return EplHeadEnum; else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum; @@ -741,5 +752,8 @@ else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; else if (strcmp(name,"NewDamage")==0) return NewDamageEnum; - else if (strcmp(name,"Node")==0) return NodeEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"Node")==0) return NodeEnum; else if (strcmp(name,"OmegaAbsGradient")==0) return OmegaAbsGradientEnum; else if (strcmp(name,"OceantransportSpcbottompressure")==0) return OceantransportSpcbottompressureEnum; @@ -752,8 +766,5 @@ else if (strcmp(name,"Radar")==0) return RadarEnum; else if (strcmp(name,"RadarAttenuationMacGregor")==0) return RadarAttenuationMacGregorEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum; + else if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum; else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum; else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum; @@ -803,4 +814,13 @@ else if (strcmp(name,"SealevelGrotm2")==0) return SealevelGrotm2Enum; else if (strcmp(name,"SealevelGrotm3")==0) return SealevelGrotm3Enum; + else if (strcmp(name,"SealevelGUrotm1")==0) return SealevelGUrotm1Enum; + else if (strcmp(name,"SealevelGUrotm2")==0) return SealevelGUrotm2Enum; + else if (strcmp(name,"SealevelGUrotm3")==0) return SealevelGUrotm3Enum; + else if (strcmp(name,"SealevelGNrotm1")==0) return SealevelGNrotm1Enum; + else if (strcmp(name,"SealevelGNrotm2")==0) return SealevelGNrotm2Enum; + else if (strcmp(name,"SealevelGNrotm3")==0) return SealevelGNrotm3Enum; + else if (strcmp(name,"SealevelGErotm1")==0) return SealevelGErotm1Enum; + else if (strcmp(name,"SealevelGErotm2")==0) return SealevelGErotm2Enum; + else if (strcmp(name,"SealevelGErotm3")==0) return SealevelGErotm3Enum; else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum; else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum; @@ -855,5 +875,8 @@ else if (strcmp(name,"SmbD")==0) return SmbDEnum; else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum; - else if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum; else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum; else if (strcmp(name,"SmbDailydsradiation")==0) return SmbDailydsradiationEnum; @@ -875,8 +898,5 @@ else if (strcmp(name,"SmbEC")==0) return SmbECEnum; else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; + else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; else if (strcmp(name,"SmbEla")==0) return SmbElaEnum; else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum; @@ -978,5 +998,8 @@ else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum; else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum; - else if (strcmp(name,"Temperature")==0) return TemperatureEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Temperature")==0) return TemperatureEnum; else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum; else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; @@ -998,8 +1021,5 @@ else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; else if (strcmp(name,"VxObs")==0) return VxObsEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"VxShear")==0) return VxShearEnum; + else if (strcmp(name,"VxShear")==0) return VxShearEnum; else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum; else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; @@ -1101,5 +1121,8 @@ else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum; - else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum; else if (strcmp(name,"Outputdefinition80")==0) return Outputdefinition80Enum; else if (strcmp(name,"Outputdefinition81")==0) return Outputdefinition81Enum; @@ -1121,8 +1144,5 @@ else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum; else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum; + else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum; else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum; else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum; @@ -1224,5 +1244,8 @@ else if (strcmp(name,"FileParam")==0) return FileParamEnum; else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum; - else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum; else if (strcmp(name,"FloatingMeltRate")==0) return FloatingMeltRateEnum; @@ -1244,8 +1267,5 @@ else if (strcmp(name,"Gradient1")==0) return Gradient1Enum; else if (strcmp(name,"Gradient2")==0) return Gradient2Enum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"Gradient3")==0) return Gradient3Enum; + else if (strcmp(name,"Gradient3")==0) return Gradient3Enum; else if (strcmp(name,"Gradient4")==0) return Gradient4Enum; else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum; @@ -1347,5 +1367,8 @@ else if (strcmp(name,"MinVy")==0) return MinVyEnum; else if (strcmp(name,"MinVz")==0) return MinVzEnum; - else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum; else if (strcmp(name,"Moulin")==0) return MoulinEnum; else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; @@ -1367,8 +1390,5 @@ else if (strcmp(name,"OceantransportSolution")==0) return OceantransportSolutionEnum; else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum; + else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum; else if (strcmp(name,"Open")==0) return OpenEnum; else if (strcmp(name,"Option")==0) return OptionEnum; @@ -1470,5 +1490,8 @@ else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum; else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum; - else if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum; + else stage=13; + } + if(stage==13){ + if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum; else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum; else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum; @@ -1490,8 +1513,5 @@ else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum; else if (strcmp(name,"Vertices")==0) return VerticesEnum; - else stage=13; - } - if(stage==13){ - if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; + else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum; else if (strcmp(name,"Water")==0) return WaterEnum; else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum; Index: /issm/trunk-jpl/src/c/shared/Numerics/types.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Numerics/types.h (revision 26241) +++ /issm/trunk-jpl/src/c/shared/Numerics/types.h (revision 26242) @@ -13,4 +13,5 @@ #include +#include /*here are our abstracted types: inspired on petsc */ @@ -24,17 +25,23 @@ /*ADOLC typedefs*/ #include "adolc/adolc.h" -typedef adouble IssmDouble; /*for active variables*/ -typedef double IssmPDouble; /*for passive variables*/ +typedef adouble IssmDouble; /*for active variables*/ +typedef acomplex IssmComplex; /*for active variables*/ +typedef double IssmPDouble; /*for passive variables*/ +typedef std::complex IssmPComplex; /*for passive variables*/ #elif defined(_HAVE_CODIPACK_) && !defined(_WRAPPERS_) /*CoDiPack typedefs*/ #include -typedef codi::RealReverse IssmDouble; -typedef double IssmPDouble; +typedef codi::RealReverse IssmDouble; +//typedef codi::ComplexReverse IssmComplex; //not defined so far +typedef double IssmPDouble; +typedef std::complex IssmPComplex; /*Non-AD typedefs*/ #else -typedef double IssmDouble; -typedef IssmDouble IssmPDouble; +typedef double IssmDouble; +typedef std::complex IssmComplex; +typedef IssmDouble IssmPDouble; +typedef IssmComplex IssmPComplex; #endif Index: /issm/trunk-jpl/src/m/classes/fourierlove.m =================================================================== --- /issm/trunk-jpl/src/m/classes/fourierlove.m (revision 26241) +++ /issm/trunk-jpl/src/m/classes/fourierlove.m (revision 26242) @@ -19,9 +19,11 @@ istemporal = 0; n_temporal_iterations = 0; - time = 0; + time = 0; love_kernels = 0; forcing_type = 0; - inner_core_boundary = 0; + inner_core_boundary = 0; core_mantle_boundary = 0; + complex_computation = 0; + end methods (Static) @@ -63,4 +65,5 @@ self.inner_core_boundary=1; self.core_mantle_boundary=2; + self.complex_computation=0; end % }}} function disp(self) % {{{ @@ -102,4 +105,5 @@ md = checkfield(md,'fieldname','love.love_kernels','values',[0 1]); md = checkfield(md,'fieldname','love.forcing_type','NaN',1,'Inf',1,'numel',1,'>',0, '<=', 12); + md = checkfield(md,'fieldname','love.complex_computation','NaN',1,'Inf',1,'numel',1,'values',[0 1]); md = checkfield(md,'fieldname','love.istemporal','values',[0 1]); @@ -140,4 +144,5 @@ WriteData(fid,prefix,'object',self,'fieldname','istemporal','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','n_temporal_iterations','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','complex_computation','format','Boolean'); %note: no need to marshall the time vector, we have frequencies WriteData(fid,prefix,'object',self,'fieldname','love_kernels','format','Boolean'); Index: /issm/trunk-jpl/src/m/solve/loadresultfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/solve/loadresultfromdisk.m (revision 26241) +++ /issm/trunk-jpl/src/m/solve/loadresultfromdisk.m (revision 26242) @@ -49,4 +49,9 @@ N=fread(fid,1,'int'); field=fread(fid,[N M],'int')'; + elseif type==5, + N=fread(fid,1,'int'); + fieldr=fread(fid,[N M],'double')'; + fieldi=fread(fid,[N M],'double')'; + field=complex(fieldr,fieldi); else error(['cannot read data of type ' num2str(type) ]); Index: /issm/trunk-jpl/test/NightlyRun/test2084.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2084.m (revision 26241) +++ /issm/trunk-jpl/test/NightlyRun/test2084.m (revision 26242) @@ -26,5 +26,5 @@ md.materials.lame_lambda=md.materials.lame_mu*0+5e17; md.materials.issolid=[1 0 1 1 1 1]'; -md.materials.rheologymodel=2*ones(md.materials.numlayers,1); +md.materials.rheologymodel=0*ones(md.materials.numlayers,1); md.materials.burgers_mu=md.materials.lame_mu/3; md.materials.burgers_viscosity=md.materials.viscosity/10;