Changeset 20047

Timestamp:

02/02/16 15:28:39 (9 years ago)

Author:

seroussi

Message:

NEW: finished mantle class

File:

• issm/trunk-jpl/src/c/classes/Elements/Element.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -376 +376 @@
 
         /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
         delete gauss;
         return divergence;
@@ -1574 +1575 @@
         IssmDouble  mantleconductivity,nusselt,dtbg,plumeradius,topplumedepth,bottomplumedepth,plumex,plumey;
         IssmDouble  crustthickness,uppercrustthickness,uppercrustheat,lowercrustheat;
-        //IssmDouble* base     = xNew<IssmDouble>(numvertices);
-        //IssmDouble* bed      = xNew<IssmDouble>(numvertices);
+        IssmDouble  crustheat,plumeheat,dt,middleplumedepth,a,e,eprime,A0,lambda,Alambda,dAlambda;
+        IssmDouble  x,y,z,c;
         IssmDouble* values   = xNew<IssmDouble>(numvertices);
+        IssmDouble *xyz_list = NULL;
 
         parameters->FindParam(&mantleconductivity,BasalforcingsMantleconductivityEnum);
@@ -1591 +1593 @@
         parameters->FindParam(&lowercrustheat,BasalforcingsLowercrustheatEnum);
 
-        //this->GetInputListOnVertices(base,BaseEnum);
-        //this->GetInputListOnVertices(bed,BedEnum);
+        this->GetVerticesCoordinates(&xyz_list);
+        c=plumeradius;
+        a=(bottomplumedepth-topplumedepth)/2.;
+        e=pow(a*a-c*c,1./2.)/a;
+        A0=(1-e*e)/(e*e*e)*(1/2*log((1+e)/(1-e))-e);
+        middleplumedepth=-topplumedepth-crustthickness-a;
         for(int i=0;i<numvertices;i++){
-                values[i]=uppercrustheat*uppercrustthickness+ lowercrustheat*(crustthickness-uppercrustthickness);
+                y=xyz_list[i*3+0]-plumex;
+                z=xyz_list[i*3+1]-plumey;
+                x=(a+topplumedepth+crustthickness);
+                lambda=(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))/2;
+                dAlambda=(-8*a*pow(c,2)*x*(-2*pow(a,2)+2*pow(c,2)+sqrt(2)*sqrt((a-c)*(a+c))*sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))*(pow(a,4)*(pow(y,2)+pow(z,2))+pow(c,4)*(pow(y,2)+pow(z,2))+pow(pow(x,2)+pow(y,2)+pow(z,2),2)*(pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))+pow(c,2)*(pow(x,4)-pow(x,2)*(pow(y,2)+pow(z,2))-(pow(y,2)+pow(z,2))*(2*pow(y,2)+2*pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))+pow(a,2)*(-pow(x,4)+pow(x,2)*(pow(y,2)+pow(z,2))+(pow(y,2)+pow(z,2))*(-2*pow(c,2)+2*(pow(y,2)+pow(z,2))+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))))/(sqrt((a-c)*(a+c))*sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))*pow(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))),3.5)*pow(-(sqrt(2)*sqrt((a-c)*(a+c)))+sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))),2)*(sqrt(2)*sqrt((a-c)*(a+c))+sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))));
+                eprime=pow((a*a-plumeradius*plumeradius)/(a*a+lambda),1./2.);
+                Alambda=(1.-e*e)/(e*e*e)*(1./2.*log((1.+eprime)/(1.-eprime))-eprime);
+                dt=dtbg-(nusselt-1.)/(1.+A0*(nusselt-1.))*(Alambda*dtbg+middleplumedepth*dtbg*dAlambda);
+                plumeheat=mantleconductivity*dt;
+                crustheat=uppercrustheat*uppercrustthickness+lowercrustheat*(crustthickness-uppercrustthickness);
+                values[i]=crustheat+plumeheat;
         }
 
         this->AddInput(BasalforcingsGeothermalfluxEnum,values,P1Enum);
-        //xDelete<IssmDouble>(base);
-        //xDelete<IssmDouble>(bed);
+        xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(values);
 
@@ -2669 +2684 @@
         /*Clean up and return*/
         xDelete<IssmDouble>(maxprincipal);
+        xDelete<IssmDouble>(xyz_list);
         delete gauss;
 }
@@ -3098 +3114 @@
         /*Clean up and return*/
         xDelete<IssmDouble>(viscousheating);
+        xDelete<IssmDouble>(xyz_list);
         delete gauss;
 }