Changeset 18505

Timestamp:

09/11/14 08:39:35 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: finished implementing post processing... debuging to come

Location:

issm/trunk-jpl/src/c

Files:

• analyses/SeaiceAnalysis.cpp (modified) (7 diffs)
• analyses/SeaiceAnalysis.h (modified) (1 diff)
• cores/seaice_core.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp

@@ -141 +141 @@
                 for(int i=0;i<3;i++){
                         for(int j=0;j<3;j++){
-                                C[i][j] = dt*time_ratio*gauss->weight*Jdet*C[i][j];
+                                C[i][j] = dt*thickness*time_ratio*gauss->weight*Jdet*C[i][j];
                         }
                 }
@@ -320 +320 @@
 
         /*Intermediaries*/
-        int*        doflist=NULL;
-        IssmDouble* xyz_list=NULL;
+        int* doflist=NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -358 +357 @@
         xDelete<IssmDouble>(vx);
         xDelete<IssmDouble>(values);
-        xDelete<IssmDouble>(xyz_list);
         xDelete<int>(doflist);
 }/*}}}*/
@@ -416 +414 @@
         /*Intermediaries*/
         IssmDouble D[3][3];
-        IssmDouble damage,thickness,concentration;
+        IssmDouble damage,concentration;
 
         /*Get Material parameters*/
@@ -425 +423 @@
         /*Get damage input at this location*/
         Input* damage_input        = element->GetInput(DamageEnum);              _assert_(damage_input);
-        Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);     _assert_(thickness_input);
         Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
         damage_input->GetInputValue(&damage,gauss);
-        thickness_input->GetInputValue(&thickness,gauss);
         concentration_input->GetInputValue(&concentration,gauss);
 
@@ -448 +444 @@
         for(int i=0;i<3;i++){
                 for(int j=0;j<3;j++){
-                        C[i*3+j] = thickness*E*D[i][j];
+                        C[i*3+j] = E*D[i][j];
                 }
         }
@@ -623 +619 @@
         }
 }/*}}}*/
+void SeaiceAnalysis::PostProcess(FemModel* femmodel){/*{{{*/
+
+        /*Intermediaties*/
+        IssmDouble *xyz_list = NULL;
+        IssmDouble *xyz_list_old = NULL;
+        IssmDouble *xyz_list_new = NULL;
+        IssmDouble  dt,meshx,meshy,vx,vy,eps_xx,eps_yy,eps_zz,area_old,area_new;
+        IssmDouble  time_relaxation_stress,time_relaxation_damage,damage,concentration,thickness;
+        IssmDouble  sigma_dot_xx,sigma_dot_yy,sigma_dot_xy;
+        IssmDouble  sigma_pred_xx,sigma_pred_yy,sigma_pred_xy;
+        IssmDouble  sigma_xx,sigma_yy,sigma_xy;
+        IssmDouble  min_c,min_h,max_h;
+        IssmDouble  epsilon[3]; /*eps_xx,eps_yy,eps_xy*/
+        IssmDouble  C[3][3];
+
+        /*Fetch the parameters once for all*/
+        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+        femmodel->parameters->FindParam(&min_c,SeaiceMinConcentrationEnum);
+        femmodel->parameters->FindParam(&min_h,SeaiceMinThicknessEnum);
+        femmodel->parameters->FindParam(&max_h,SeaiceMaxThicknessEnum);
+
+        /*Loop over the elements of this partition and update accordingly*/
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+
+                /*Get some inputs needed for the update*/
+                element->GetVerticesCoordinates(&xyz_list);
+                time_relaxation_stress = element->GetMaterialParameter(MaterialsTimeRelaxationStressEnum);
+                time_relaxation_damage = element->GetMaterialParameter(MaterialsTimeRelaxationDamageEnum);
+                Input* meshx_input = element->GetInput(MeshXEnum); _assert_(meshx_input);
+                Input* meshy_input = element->GetInput(MeshYEnum); _assert_(meshy_input);
+                Input* vx_input    = element->GetInput(VxEnum);    _assert_(vx_input);
+                Input* vy_input    = element->GetInput(VyEnum);    _assert_(vy_input);
+
+                /*Preallocate future inputs*/
+                int numvertices=element->GetNumberOfVertices();
+                IssmDouble* meshx_new_list = xNew<IssmDouble>(numvertices);
+                IssmDouble* meshy_new_list = xNew<IssmDouble>(numvertices);
+                xyz_list_old = xNew<IssmDouble>(numvertices*3);
+                xyz_list_new = xNew<IssmDouble>(numvertices*3);
+
+                /*1. update vertex positions (using a Gauss object for convenience)*/
+                Gauss* gauss=element->NewGauss();
+                for (int iv=0;iv<numvertices;iv++){
+                        gauss->GaussVertex(iv);
+
+                        meshx_input->GetInputValue(&meshx,gauss);
+                        meshy_input->GetInputValue(&meshy,gauss);
+                        vx_input->GetInputValue(&vx,gauss);
+                        vy_input->GetInputValue(&vy,gauss);
+
+                        meshx_new_list[iv] = meshx + vx*dt;
+                        meshy_new_list[iv] = meshy + vy*dt;
+
+                        xyz_list_old[iv*3+0] = meshx;
+                        xyz_list_old[iv*3+1] = meshy;
+                        xyz_list_old[iv*3+2] = 0.;
+
+                        xyz_list_new[iv*3+0] = meshx_new_list[iv];
+                        xyz_list_new[iv*3+1] = meshy_new_list[iv];
+                        xyz_list_new[iv*3+2] = 0.;
+                }
+                element->AddInput(MeshXEnum,meshx_new_list,P1Enum);
+                element->AddInput(MeshYEnum,meshy_new_list,P1Enum);
+
+                /*Now we are going to use a point in the center of the element, because we are P0 for the stress*/
+                gauss->GaussNode(P0Enum,0);
+
+                /*Calculate sigma_dot*/
+                sigma_dot_xx=sigma_dot_yy=sigma_dot_xy=0.;
+                this->CreateCTensor(&C[0][0],element,gauss);
+                element->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+                for(int j=0;j<3;j++){
+                        sigma_dot_xx += C[0][j]*epsilon[0];
+                        sigma_dot_yy += C[1][j]*epsilon[1];
+                        sigma_dot_xy += C[2][j]*epsilon[2];
+                }
+
+                /*Get predicted stress state*/
+                Input* sigma_xx_input = element->GetInput(StressTensorPredictorxxEnum); _assert_(sigma_xx_input);
+                Input* sigma_yy_input = element->GetInput(StressTensorPredictoryyEnum); _assert_(sigma_yy_input);
+                Input* sigma_xy_input = element->GetInput(StressTensorPredictorxyEnum); _assert_(sigma_xy_input);
+                sigma_xx_input->GetInputAverage(&sigma_xx);
+                sigma_yy_input->GetInputAverage(&sigma_yy);
+                sigma_xy_input->GetInputAverage(&sigma_xy);
+
+                /*Calculate new stress and push to element*/
+                sigma_xx = sigma_xx + dt*(sigma_dot_xx - sigma_xx/time_relaxation_stress);
+                sigma_yy = sigma_yy + dt*(sigma_dot_yy - sigma_yy/time_relaxation_stress);
+                sigma_xy = sigma_xy + dt*(sigma_dot_xy - sigma_xy/time_relaxation_stress);
+                element->AddInput(StressTensorxxEnum,&sigma_xx,P0Enum);
+                element->AddInput(StressTensoryyEnum,&sigma_yy,P0Enum);
+                element->AddInput(StressTensorxyEnum,&sigma_xy,P0Enum);
+
+                /*Update Damage According*/
+                Input* damage_input = element->GetInput(DamageEnum); _assert_(damage_input);
+                damage_input->GetInputAverage(&damage);
+                damage = damage*(1.-dt/time_relaxation_damage);
+                element->AddInput(DamageEnum,&damage,P0Enum);
+
+                /*Prepare new predictor*/
+                this->CreateCTensor(&C[0][0],element,gauss); /*C is updated now that there is a new damage in inputs*/
+                for(int j=0;j<3;j++){
+                        sigma_pred_xx += C[0][j]*epsilon[0];
+                        sigma_pred_yy += C[1][j]*epsilon[1];
+                        sigma_pred_xy += C[2][j]*epsilon[2];
+                }
+                element->AddInput(StressTensorPredictorxxEnum,&sigma_pred_xx,P0Enum);
+                element->AddInput(StressTensorPredictoryyEnum,&sigma_pred_yy,P0Enum);
+                element->AddInput(StressTensorPredictorxyEnum,&sigma_pred_xy,P0Enum);
+
+                /*Calculate Old and new area (FIXME: for now we assume trianlges...)*/
+                element->JacobianDeterminant(&area_old,xyz_list_old,gauss);
+                element->JacobianDeterminant(&area_new,xyz_list_new,gauss);
+                area_old = area_old/SQRT3;
+                area_new = area_new/SQRT3;
+
+                /*Update ice thickness and concentration using element distortion*/
+                Input* thickness_input     = element->GetInput(SeaiceThicknessEnum);     _assert_(thickness_input);
+                Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input);
+                thickness_input->GetInputAverage(&thickness);
+                concentration_input->GetInputAverage(&concentration);
+                IssmDouble ice_area   = concentration*area_old;
+                IssmDouble ice_volume = thickness*area_old;
+                if(concentration>min_c){
+                        concentration = ice_area  /area_new;
+                        thickness     = ice_volume/area_new;
+
+                        /* lower bounds */
+                        concentration = ((concentration>min_c)?(concentration):(min_c));
+                        thickness     = ((thickness    >min_h)?(thickness    ):(min_h));
+
+                        /* upper bounds (only for the concentration) */
+                        concentration    = ((concentration<1.   )?(concentration):(1));
+                        thickness        = ((thickness    <max_h)?(thickness    ):(max_h));
+                }
+                element->AddInput(SeaiceConcentrationEnum,&concentration,P0Enum);
+                element->AddInput(SeaiceThicknessEnum,&thickness,P0Enum);
+
+                /*Clean up*/
+                xDelete<IssmDouble>(meshx_new_list);
+                xDelete<IssmDouble>(meshy_new_list);
+                xDelete<IssmDouble>(xyz_list);
+                xDelete<IssmDouble>(xyz_list_old);
+                xDelete<IssmDouble>(xyz_list_new);
+        }
+
+}/*}}}*/

issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.h

@@ -36 +36 @@
                 void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
                 void GetM(IssmDouble* M,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+                void PostProcess(FemModel* femmodel);
 };
 #endif

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

@@ -13 +13 @@
 
         /*parameters: */
-        bool save_results;
+        bool           save_results;
+        SeaiceAnalysis analysis = SeaiceAnalysis();
 
         /*recover parameters: */
@@ -20 +21 @@
         /*Launch solution sequence for the only analysis we are interested in*/
         femmodel->SetCurrentConfiguration(SeaiceAnalysisEnum);
-        SeaiceAnalysis* analysis = new SeaiceAnalysis();
-        analysis->UpdateDamageAndStress(femmodel);
+        analysis.UpdateDamageAndStress(femmodel);
         solutionsequence_linear(femmodel);
+        analysis.PostProcess(femmodel);
 
         if(save_results){
                 if(VerboseSolution()) _printf0_("   saving results\n");
-                int outputs[6] = {VxEnum,VyEnum,VelEnum,SeaiceConcentrationEnum,SeaiceThicknessEnum,DamageEnum};
-                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],6);
+                int outputs[14] = {VxEnum,VyEnum,VelEnum,
+                        SeaiceConcentrationEnum,SeaiceThicknessEnum,DamageEnum,
+                        StressTensorPredictorxxEnum,StressTensorPredictoryyEnum,StressTensorPredictorxyEnum,
+                        StressTensorxxEnum,StressTensoryyEnum,StressTensorxyEnum,
+                        MeshXEnum,MeshYEnum};
+                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],14);
         }