Changeset 18347

Timestamp:

08/08/14 09:27:46 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: finished Zalesak, still one pb with umin and umax

File:

• issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp

@@ -167 +167 @@
         IssmDouble* udot_serial = NULL;
         IssmDouble* u_serial    = NULL;
-        VecToMPISerial(&udot_serial,udot,IssmComm::GetComm());
-        VecToMPISerial(&u_serial   ,u   ,IssmComm::GetComm());
-
+        IssmDouble* ml_serial    = NULL;
+        VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
+        VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
+        VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
 
         /*Anti diffusive fluxes*/
-        Mat F    = NULL;
+        Vec Ri_plus  = NULL;
+        Vec Ri_minus = NULL;
+        double uiLmax = 3.;
+        double uiLmin = 2.;
+        VecDuplicate(u,&Ri_plus);
+        VecDuplicate(u,&Ri_minus);
+        for(int row=rstart; row<rend; row++){
+                double Pi_plus  = 0.;
+                double Pi_minus = 0.;
+                MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
+                MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+                _assert_(ncols==ncols2);
+                for(int j=0; j<ncols; j++) {
+                        _assert_(cols[j]==cols2[j]);
+                        d = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
+                        if(row!=cols[j]){
+                                if(d>0.){
+                                        Pi_plus  += d;
+                                }
+                                else{
+                                        Pi_minus += d;
+                                }
+                        }
+                }
+
+                /*Compute Qis and Ris*/
+                double Qi_plus  = ml_serial[row]/deltat*(uiLmax - u_serial[row]);
+                double Qi_minus = ml_serial[row]/deltat*(uiLmin - u_serial[row]);
+                d = 1.;
+                if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
+                VecSetValue(Ri_plus,row,(const double)d,INSERT_VALUES);
+                d = 1.;
+                if(Pi_minus!=0.) d = min(1.,Qi_minus/Pi_minus);
+                VecSetValue(Ri_minus,row,(const double)d,INSERT_VALUES);
+
+                MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
+                MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+        }
+        VecAssemblyBegin(Ri_plus);
+        VecAssemblyEnd(  Ri_plus);
+        VecAssemblyBegin(Ri_minus);
+        VecAssemblyEnd(  Ri_minus);
+
+        /*Serialize Ris*/
+        IssmDouble* Ri_plus_serial  = NULL;
+        IssmDouble* Ri_minus_serial = NULL;
+        VecToMPISerial(&Ri_plus_serial,Ri_plus,IssmComm::GetComm());
+        VecToMPISerial(&Ri_minus_serial,Ri_minus,IssmComm::GetComm());
+        VecFree(&Ri_plus);
+        VecFree(&Ri_minus);
+
+        /*Build fbar*/
         Vec Fbar = NULL;
         VecDuplicate(u,&Fbar);
-        MatDuplicate(D_petsc,MAT_SHARE_NONZERO_PATTERN,&F);
         for(int row=rstart; row<rend; row++){
                 MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
                 MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
                 _assert_(ncols==ncols2);
-                mi = 0.;
-                for(int j=0; j<ncols; j++) {
-                        _assert_(cols[j]==cols2[j]);
-                        d = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
-                        if(row!=cols[j]) mi += d;
-                        MatSetValues(F,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
-
-                }
-                VecSetValue(Fbar,row,(const double)mi,INSERT_VALUES);
+                d = 0.;
+                for(int j=0; j<ncols; j++) {
+                        _assert_(cols[j]==cols2[j]);
+                        if(row==cols[j]) continue;
+                        double f_ij = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
+                        if(f_ij>0){
+                                d += min(Ri_plus_serial[row],Ri_minus_serial[cols[j]]) * f_ij;
+                        }
+                        else{
+                                d += min(Ri_minus_serial[row],Ri_plus_serial[cols[j]]) * f_ij;
+                        }
+                }
+                VecSetValue(Fbar,row,(const double)d,INSERT_VALUES);
                 MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
                 MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
         }
-        MatAssemblyBegin(F,MAT_FINAL_ASSEMBLY);
-        MatAssemblyEnd(  F,MAT_FINAL_ASSEMBLY);
         VecAssemblyBegin(Fbar);
         VecAssemblyEnd(  Fbar);
 
         MatFree(&D_petsc);
-        MatFree(&F);
         delete Mc;
         xDelete<IssmDouble>(udot_serial);
         xDelete<IssmDouble>(u_serial);
+        xDelete<IssmDouble>(ml_serial);
+        xDelete<IssmDouble>(Ri_plus_serial);
+        xDelete<IssmDouble>(Ri_minus_serial);
 
         /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/