Changeset 26249

Timestamp:

05/06/21 20:44:39 (4 years ago)

Author:

caronlam

Message:

BUG: in love_core: dydx was not reinitialized at each call of yi_derivatives; la, mu not properly handled in all cases in GetEarthRheology

File:

• issm/trunk-jpl/src/c/cores/love_core.cpp (modified) (14 diffs)

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

@@ -164 +164 @@
 
         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=la0 + 2.0/3.0*mu0; //Bulk modulus
                 if (rheo==2){//EBM
@@ -203 +202 @@
 
                         mu=mu0*omega*(omega+mu2/vi2)/((omega+mu2/vi2)*(omega+mu0/vi)+mu0/vi2*omega);
-                        la=ka-2.0/3.0*mu0;
+                        la=ka-2.0/3.0*mu;
                 }
                 else{//Maxwell
@@ -210 +209 @@
                 }
         }
-
+        else{//Otherwise return the elastic value
+        la=la0;
+        mu=mu0;
+        }
         *pla=la;
         *pmu=mu;
@@ -444 +446 @@
                            dy[1*6+1]= -2.0/x-fgr/rg0;
                            }
-                           */ /*}}}*/
+        */ /*}}}*/
 
         if(issolid==1){
@@ -453 +455 @@
 
         for (id=0;id<ny;id++){
+                dydx[id]=0.0;
                 for (iy=0;iy<ny;iy++){
                         nindex=layer_index*nstep*36+n*36;
-                        dydx[id]=dydx[id]+yi_prefactor[nindex+id*6+iy]*y[iy];
+                        dydx[id]+=yi_prefactor[nindex+id*6+iy]*y[iy];
                 }
         }
         return;
 }/*}}}*/
-template <typename doubletype> void        propagate_yi_euler(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+template <typename doubletype> void        propagate_yi_euler(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?
         //euler method
@@ -466 +469 @@
         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); 
 
-        IssmDouble* dydx=xNewZeroInit<IssmDouble>(6);
+        doubletype* dydx=xNewZeroInit<doubletype>(6);
         IssmDouble dr = (xmax -xmin)/nstep;
         IssmDouble x=xmin;
@@ -476 +479 @@
                 x = x + dr;
         }
-        xDelete<IssmDouble>(dydx);
+        xDelete<doubletype>(dydx);
 }/*}}}*/
 template <typename doubletype> void        propagate_yi_RK2(doubletype* y, IssmDouble xmin, IssmDouble xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
@@ -681 +684 @@
                 if (matlitho->issolid[i]){
                         one = -1.0;
-                        if (i>0 && !matlitho->issolid[i-1]) one = 1.0;
+                        if (i>0) if (!matlitho->issolid[i-1]) one = 1.0;
                         for (int j=0;j<6;j++){
                                 yi[(j+6*ici)+ nyi*(j+6*ici+3)] = one;
@@ -891 +894 @@
                 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
                         if(VerboseSolution()){
-                                printf("\n   Degree: %i, surface/Depth Love number ratio small: %e\n",deg, loveratio);
-                                printf("    Changing the interface where the integration starts\n");
-                                printf("    New start interface: r=%gm\n\n", matlitho->radius[starting_layer+1]);
+                                _printf_("\n   Degree: " << deg <<", surface/Depth Love number ratio small: "<<loveratio<<"\n");
+                                _printf_("    Changing the interface where the integration starts\n");
+                                _printf_("    New start interface: r="<< matlitho->radius[starting_layer+1] <<"m\n\n");
                         }
                         nyi-=6;
@@ -1005 +1008 @@
         bool        allow_layer_deletion,love_kernels, istemporal;
         bool        verbosemod = (int)VerboseModule();
-        IssmDouble  g0,r0,mu0, underflow_tol;
+        IssmDouble  mu0, underflow_tol;
         IssmDouble *frequencies = NULL;
         bool        save_results;
@@ -1040 +1043 @@
         femmodel->parameters->FindParam(&sh_nmax,LoveShNmaxEnum);
         femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum);
-        femmodel->parameters->FindParam(&g0,LoveG0Enum);
-        femmodel->parameters->FindParam(&r0,LoveR0Enum);
         femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
         femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
@@ -1126 +1127 @@
         femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveLiEnum,LoveLf,sh_nmax+1,nfreq,0,0));
 
+        /*Only when love_kernels is on*/
+        if (love_kernels==1) {
+                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernels,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
+        }
+
+        /*Free resources:*/
+        xDelete<IssmDouble>(frequencies);
+        xDelete<doubletype>(LoveKf);
+        xDelete<doubletype>(LoveHf);
+        xDelete<doubletype>(LoveLf);
+        xDelete<doubletype>(LoveKernels);
+
+        /* Legacy for fortran core, to be removed after complete validation */
+
+        //IssmDouble g0,r0;
+        //femmodel->parameters->FindParam(&g0,LoveG0Enum);
+        //femmodel->parameters->FindParam(&r0,LoveR0Enum);
+        //IssmDouble* LoveKr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+        //IssmDouble* LoveHr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+        //IssmDouble* LoveLr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+        //IssmDouble* LoveKi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+        //IssmDouble* LoveHi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+        //IssmDouble* LoveLi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+
+        /*Initialize love kernels (real and imaginary parts): */
+        //IssmDouble* LoveKernelsr= xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
+        //IssmDouble* LoveKernelsi= xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
 
         /*call the main module: */
         //if (false){
-        //FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHi,LoveLr,LoveLi,LoveKernelsReal,LoveKernelsImag,  //output
+        //FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHr,LoveLr,LoveLi,LoveKernelsr,LoveKernelsi,  //output
         //              nfreq,frequencies,sh_nmax,sh_nmin,g0,r0,mu0,allow_layer_deletion,forcing_type,verbosemod, //parameter inputs
         //              matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu,
         //              matlitho->burgers_viscosity, matlitho->burgers_mu, matlitho->density, matlitho->rheologymodel, matlitho->issolid //matlitho inputs
         //              );
+
+
+        /*Only when love_kernels is on*/
+        //if (love_kernels==1) {
+                //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernelsr,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
+                //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
         //}
+
 
         /*Add love matrices to results:*/
@@ -1144 +1179 @@
         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0));
 
-        /*Only when love_kernels is on*/
-        if (love_kernels==1) {
-                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernels,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
-                //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
-        }
-
-        /*Free resources:*/
-        xDelete<IssmDouble>(frequencies);
-        xDelete<doubletype>(LoveKf);
-        xDelete<doubletype>(LoveHf);
-        xDelete<doubletype>(LoveLf);
+
+
+        //xDelete<IssmDouble>(LoveKr);
+        //xDelete<IssmDouble>(LoveHr);
+        //xDelete<IssmDouble>(LoveLr);
+        //xDelete<IssmDouble>(LoveKernelsr);
+        //xDelete<IssmDouble>(LoveKi);
+        //xDelete<IssmDouble>(LoveHi);
         //xDelete<IssmDouble>(LoveLi);
-        xDelete<doubletype>(LoveKernels);
+        //xDelete<IssmDouble>(LoveKernelsi);
+
 } /*}}}*/
 
@@ -1173 +1206 @@
 template void        propagate_yi_RK2<IssmDouble>(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
 template void        propagate_yi_RK2<IssmComplex>(IssmComplex* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
+template void        propagate_yi_euler<IssmDouble>(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
+template void        propagate_yi_euler<IssmComplex>(IssmComplex* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
 template void        Innersphere_boundaryconditions<IssmDouble>(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
 template void        Innersphere_boundaryconditions<IssmComplex>(IssmComplex* yi, int layer_index, int deg, IssmComplex omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);