Changeset 27490

Timestamp:

01/01/23 08:43:06 (2 years ago)

Author:

Mathieu Morlighem

Message:

CHG: deal with issolid and rheologymodel: they should be bool* and int* respectively, otherwise AD gives us a lot of warnings... and some compilers too!

Location:

issm/trunk-jpl/src/c

Files:

• classes/IoModel.cpp (modified) (1 diff)
• classes/IoModel.h (modified) (1 diff)
• classes/Materials/Matlitho.cpp (modified) (4 diffs)
• classes/Materials/Matlitho.h (modified) (1 diff)
• cores/love_core.cpp (modified) (9 diffs)
• modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/classes/IoModel.cpp

@@ -1124 +1124 @@
 }
 /*}}}*/
+void  IoModel::FetchData(bool** pmatrix,int* pM,int* pN,const char* data_name){/*{{{*/
+
+        /*output: */
+        int M,N;
+        IssmPDouble* matrix=NULL;
+        bool*        bool_matrix=NULL;
+        int code=0;
+
+        /*recover my_rank:*/
+        int my_rank=IssmComm::GetRank();
+
+        /*Set file pointer to beginning of the data: */
+        fid=this->SetFilePointerToData(&code,NULL,data_name);
+
+        if(code!=5 && code!=6 && code!=7)_error_("expecting a IssmDouble, integer or boolean matrix for \""<<data_name<<"\""<<" (Code is "<<code<<")");
+
+        /*Now fetch: */
+
+        /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
+        /*numberofelements: */
+        if(my_rank==0){
+                if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
+        }
+
+        ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+        if(my_rank==0){
+                if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
+        }
+        ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+        /*Now allocate matrix: */
+        if(M*N){
+                matrix=xNew<IssmPDouble>(M*N);
+
+                /*Read matrix on node 0, then broadcast: */
+                if(my_rank==0){
+                        if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
+                }
+
+                ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+        }
+
+        /*Now cast to bool: */
+        if(M*N){
+                bool_matrix=xNew<bool>(M*N);
+                for(int i=0;i<M;i++){
+                        for(int j=0;j<N;j++){
+                                bool_matrix[i*N+j]=(bool)matrix[i*N+j];
+                        }
+                }
+        }
+        else{
+                bool_matrix=NULL;
+        }
+        /*Free resources:*/
+        xDelete<IssmPDouble>(matrix);
+
+        /*Assign output pointers: */
+        *pmatrix=bool_matrix;
+        if (pM)*pM=M;
+        if (pN)*pN=N;
+
+}
+/*}}}*/
 void  IoModel::FetchData(int** pmatrix,int* pM,int* pN,const char* data_name){/*{{{*/
         int i,j;

issm/trunk-jpl/src/c/classes/IoModel.h

@@ -136 +136 @@
                 void        FetchData(char*** pstrings,int* pnumstrings,const char* data_name);
                 void        FetchData(int** pmatrix,int* pM,int* pN,const char* data_name);
+                void        FetchData(bool**  pboolmatrix,int* pM,int* pN,const char* data_name);
                 void        FetchData(IssmDouble**  pscalarmatrix,int* pM,int* pN,const char* data_name);
 #if _HAVE_AD_  && !defined(_WRAPPERS_)

issm/trunk-jpl/src/c/classes/Materials/Matlitho.cpp

@@ -14 +14 @@
 /*Matlitho constructors and destructor*/
 Matlitho::Matlitho(){/*{{{*/
-        this->numlayers=0;
-        this->radius=NULL;
-        this->viscosity=NULL;
-        this->lame_lambda=NULL;
-        this->lame_mu=NULL;
-        this->burgers_viscosity=NULL;
-        this->burgers_mu=NULL;
-        this->ebm_alpha=NULL;
-        this->ebm_delta=NULL;   
-        this->ebm_taul=NULL;
-        this->ebm_tauh=NULL;
-        this->density=NULL;
-        this->rheologymodel=NULL;
-        this->issolid=NULL;
-        return;
-}
-/*}}}*/
-Matlitho::Matlitho(int matlitho_mid, IoModel* iomodel){/*{{{*/
-
-        IssmDouble* rheologymodeld=NULL;
-        IssmDouble* issolidd=NULL;
+        this->numlayers         = 0;
+        this->radius            = NULL;
+        this->viscosity         = NULL;
+        this->lame_lambda       = NULL;
+        this->lame_mu           = NULL;
+        this->burgers_viscosity = NULL;
+        this->burgers_mu        = NULL;
+        this->ebm_alpha         = NULL;
+        this->ebm_delta         = NULL;
+        this->ebm_taul          = NULL;
+        this->ebm_tauh          = NULL;
+        this->density           = NULL;
+        this->rheologymodel     = NULL;
+        this->issolid           = NULL;
+        return;
+}
+/*}}}*/
+Matlitho::Matlitho(int matlitho_mid, IoModel* iomodel, bool* issolid_in, int* rheo_in){/*{{{*/
 
         this->mid=matlitho_mid;
@@ -72 +69 @@
         xMemCpy<IssmDouble>(this->density, iomodel->Data("md.materials.density"),this->numlayers);
 
-        this->rheologymodel=xNew<IssmDouble>(this->numlayers);
-        xMemCpy<IssmDouble>(this->rheologymodel, iomodel->Data("md.materials.rheologymodel"),this->numlayers);
-
-        this->issolid=xNew<IssmDouble>(this->numlayers);
-        xMemCpy<IssmDouble>(this->issolid, iomodel->Data("md.materials.issolid"),this->numlayers);
-
-        /*rheologymodeld= xNew<IssmDouble>(this->numlayers);
-        this->rheologymodel=xNew<bool>(this->numlayers);
-        xMemCpy<IssmDouble>(rheologymodeld, iomodel->Data("md.materials.rheologymodel"),this->numlayers);
-        for (int i=0;i<this->numlayers;i++)this->rheologymodel[i]=reCast<bool,IssmDouble>(rheologymodeld[i]);
-
-        issolidd= xNew<IssmDouble>(this->numlayers);
+        this->rheologymodel=xNew<int>(this->numlayers);
+        xMemCpy<int>(this->rheologymodel, rheo_in, this->numlayers);
+
         this->issolid=xNew<bool>(this->numlayers);
-        xMemCpy<IssmDouble>(issolidd, iomodel->Data("md.materials.issolid"),this->numlayers);
-        for (int i=0;i<this->numlayers;i++)this->issolid[i]=reCast<bool,IssmDouble>(issolidd[i]);*/
-
-        /*Free resources: */
-        xDelete<IssmDouble>(rheologymodeld);
-        xDelete<IssmDouble>(issolidd);
+        xMemCpy<bool>(this->issolid, issolid_in, this->numlayers);
 }
 /*}}}*/
@@ -106 +89 @@
         xDelete<IssmDouble>(ebm_tauh);
         xDelete<IssmDouble>(density);
-        xDelete<IssmDouble>(rheologymodel);
-        xDelete<IssmDouble>(issolid);
+        xDelete<int>(rheologymodel);
+        xDelete<bool>(issolid);
 
         return;
@@ -141 +124 @@
                 matlitho->ebm_tauh=xNew<IssmDouble>(this->numlayers); xMemCpy<IssmDouble>(matlitho->ebm_tauh, this->ebm_tauh,this->numlayers);
                 matlitho->density=xNew<IssmDouble>(this->numlayers); xMemCpy<IssmDouble>(matlitho->density, this->density,this->numlayers);
-                matlitho->rheologymodel=xNew<IssmDouble>(this->numlayers); xMemCpy<IssmDouble>(matlitho->rheologymodel, this->rheologymodel,this->numlayers);
-                matlitho->issolid=xNew<IssmDouble>(this->numlayers); xMemCpy<IssmDouble>(matlitho->issolid, this->issolid,this->numlayers);
-                /*matlitho->rheologymodel=xNew<bool>(this->numlayers); for(int i=0;i<this->numlayers;i++)matlitho->rheologymodel[i]=this->rheologymodel[i]; 
-                matlitho->issolid=xNew<bool>(this->numlayers); for(int i=0;i<this->numlayers;i++)matlitho->issolid[i]=this->issolid[i];*/ 
+                matlitho->rheologymodel=xNew<int>(this->numlayers); xMemCpy<int>(matlitho->rheologymodel, this->rheologymodel,this->numlayers);
+                matlitho->issolid=xNew<bool>(this->numlayers); xMemCpy<bool>(matlitho->issolid, this->issolid,this->numlayers);
         }
 

issm/trunk-jpl/src/c/classes/Materials/Matlitho.h

@@ -28 +28 @@
                 IssmDouble*  ebm_tauh;
                 IssmDouble*  density;
-                IssmDouble*  rheologymodel;
-                IssmDouble*  issolid;
+                int*         rheologymodel;
+                bool*        issolid;
 
                 Matlitho();
-                Matlitho(int matlitho_id, IoModel* iomodel);
+                Matlitho(int matlitho_id, IoModel* iomodel, bool* issolid_in, int* rheo_in);
                 ~Matlitho();
                 void SetMid(int matlitho_mid);

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

@@ -589 +589 @@
         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
+        if(vi!=0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid
                 doubletype ka=la0 + 2.0/3.0*mu0; //Bulk modulus
-                if (rheo==2){//EBM
+                if(rheo==2){//EBM
                         mu=muEBM<doubletype>(layer_index, omega, matlitho, femmodel);
                         la=ka-2.0/3.0*mu;
                 } 
-                else if (rheo==1){//Burgers
+                else if(rheo==1){//Burgers
                         doubletype vi2=matlitho->burgers_viscosity[layer_index];
                         doubletype mu2=matlitho->burgers_mu[layer_index];
@@ -674 +674 @@
         doubletype frh,frhg0,fgr0,fgr,fn,rm0,rlm,flm;
         doubletype xmin,xmax,x,dr;
-        doubletype g,ro, issolid;
+        doubletype g,ro;
+        bool       issolid;
 
         if (pomega) { //frequency and degree dependent terms /*{{{*/
@@ -690 +691 @@
                         ro=matlitho->density[layer_index];
                         issolid=matlitho->issolid[layer_index];
-                        if(issolid!=0){
+                        if(issolid){
                                 //GetEarthRheology<doubletype>(&la, &mu,layer_index,omega,matlitho);   
                                 mu=vars->mu[layer_index*vars->nfreq+vars->ifreq];
@@ -769 +770 @@
                                 g=GetGravity<doubletype>(x*ra,layer_index,femmodel,matlitho,vars);
 
-                                if(issolid==1){
+                                if(issolid){
                                         yi_prefactor[nindex+ 1*6+3]= fn/x;                  // in dy[1*6+3]
                                         yi_prefactor[nindex+ 5*6+2]= -(fgr/g0*fn)/x;        // in dy[5*6+2]
@@ -802 +803 @@
                                 g=GetGravity<doubletype>(x*ra,layer_index,femmodel,matlitho,vars);
                                 nindex=nsteps*36+n*36;
-                                if(issolid==1){
+                                if(issolid){
                                         yi_prefactor[nindex+ 1*6+5]= -frhg0;       // in dy[1*6+5]
                                         yi_prefactor[nindex+ 2*6+0]= -1.0/x;       // in dy[2*6+0]
@@ -825 +826 @@
         //computes yi derivatives at r=radius[layer_index]+ n/nstep*(radius[layer_index+1]-radius[layer_index])
 
-        int issolid=matlitho->issolid[layer_index];
+        bool issolid=matlitho->issolid[layer_index];
         int iy,id,ny, nindex, nstep, nsteps;
         //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
@@ -842 +843 @@
                            fn=(deg*(deg+1.0));
 
-                           if(issolid==1){
+                           if(issolid){
                            ny = 6;
 
@@ -899 +900 @@
         nindex=nsteps*36+n*36;
 
-        if(issolid==1){
+        if(issolid){
                 ny = 6;
         } else {
@@ -1235 +1236 @@
 
                 // Boundary Condition matrix - solid regions
-                if (matlitho->issolid[i]){
+                if(matlitho->issolid[i]){
                         one = -1.0;
-                        if (i>0) if (!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;

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

@@ -197 +197 @@
                                                 break;
                                         case MatlithoEnum: { /*{{{*/
-                                                        iomodel->FetchData(13,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.ebm_alpha","md.materials.ebm_delta","md.materials.ebm_taul","md.materials.ebm_tauh","md.materials.rheologymodel","md.materials.issolid","md.materials.density");
-                                                        materials->AddObject(new Matlitho(iomodel->numberofelements+1,iomodel));
-                                                        iomodel->DeleteData(13,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.ebm_alpha","md.materials.ebm_delta","md.materials.ebm_taul","md.materials.ebm_tauh","md.materials.rheologymodel","md.materials.issolid","md.materials.density");
+                                                        bool* issolid = NULL;
+                                                        int*  rheologymodel = NULL;
+                                                        iomodel->FetchData(&issolid, NULL, NULL, "md.materials.issolid");
+                                                        iomodel->FetchData(&rheologymodel, NULL, NULL, "md.materials.rheologymodel");
+                                                        iomodel->FetchData(11,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.ebm_alpha","md.materials.ebm_delta","md.materials.ebm_taul","md.materials.ebm_tauh","md.materials.density");
+                                                        materials->AddObject(new Matlitho(iomodel->numberofelements+1, iomodel, issolid, rheologymodel));
+                                                        iomodel->DeleteData(11,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.ebm_alpha","md.materials.ebm_delta","md.materials.ebm_taul","md.materials.ebm_tauh","md.materials.density");
+                                                        xDelete<bool>(issolid);
+                                                        xDelete<int>(rheologymodel);
                                                 }
                                                 /*}}}*/