Changeset 26101


Ignore:
Timestamp:
03/16/21 08:36:00 (4 years ago)
Author:
Mathieu Morlighem
Message:

BUG: trying to fix AD and rewrote some portions of the code

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp

    r21577 r26101  
    4444void GiaDeflectionCorex( IssmDouble* pwi, IssmDouble* pdwidt, GiaDeflectionCoreArgs* arguments){
    4545
    46         /*intermediary: */
    47         int i;
    48 
    49         /*output: */
    50         IssmDouble wi=0;
    51         IssmDouble dwidt=0;
    52 
    53         /*inputs: {{{*/
    54         /*constant: */
    55         int idisk=1; // disk #
    56         IssmDouble yts;
    57 
    58         /*coming from the model structure, runtime configurable:*/
    59         /*gia solution parameters: */
    60         int iedge=0;
    61 
    62         /*gia inputs: */
    63         IssmDouble ri; //radial distance from center of disk to vertex  i
    64         IssmDouble re; //radius of disk
    65         IssmDouble* hes; //loading history (in ice thickness)
    66         IssmDouble* times; //loading history times
    67         int numtimes; //loading history length
    68         IssmDouble currenttime;
    69         int Ntime; // number of times with load history
    70         int Ntimm; // Ntime-1 : for slope/y-cept of load segments
    71         int Ntimp; // Ntime+1 : for evaluation time 
    72         IssmDouble* blockt_time=NULL;
    73         IssmDouble* blockt_bi=NULL;
    74         IssmDouble* blockt_dmi=NULL;
    75         IssmDouble* blocky_zhload=NULL;
    76 
    77         /*gia material parameters: */
    78         IssmDouble lithosphere_shear_modulus;
    79         IssmDouble lithosphere_density;
    80         IssmDouble mantle_shear_modulus;
    81         IssmDouble mantle_viscosity;
    82         IssmDouble mantle_density;
    83         IssmDouble lithosphere_thickness;
    84 
    85         /*ice properties: */
    86         IssmDouble rho_ice;
    87 
    88         /*some debug info: */
    89         int disk_id;
    90 
    91 /*}}}*/
    92 
    9346        /*Recover material parameters and loading history: see GiaDeflectionCoreArgs for more details {{{*/
    94         ri                        = arguments->ri;
    95         re                        = arguments->re;
    96         hes                       = arguments->hes;
    97         times                     = arguments->times;
    98         numtimes                  = arguments->numtimes;
    99         currenttime               = arguments->currenttime;
    100         lithosphere_shear_modulus = arguments->lithosphere_shear_modulus;
    101         lithosphere_density       = arguments->lithosphere_density;
    102         mantle_shear_modulus      = arguments->mantle_shear_modulus;
    103         mantle_viscosity          = arguments->mantle_viscosity;
    104         mantle_density            = arguments->mantle_density;
    105         lithosphere_thickness     = arguments->lithosphere_thickness;
    106         rho_ice                   = arguments->rho_ice;
    107         disk_id                   = arguments->idisk;
    108         iedge                     = arguments->iedge;
    109         yts                       = arguments->yts;
     47        IssmDouble  ri                        = arguments->ri;                        //radial distance from center of disk to vertex i
     48        IssmDouble  re                        = arguments->re;                        //radius of disk
     49        IssmDouble *hes                       = arguments->hes;                       //loading history (in ice thickness)
     50        IssmDouble *times                     = arguments->times;                     //loading history times
     51        int         numtimes                  = arguments->numtimes;                  //loading history length
     52        IssmDouble  currenttime               = arguments->currenttime;
     53        IssmDouble  lithosphere_shear_modulus = arguments->lithosphere_shear_modulus;
     54        IssmDouble  lithosphere_density       = arguments->lithosphere_density;
     55        IssmDouble  mantle_shear_modulus      = arguments->mantle_shear_modulus;
     56        IssmDouble  mantle_viscosity          = arguments->mantle_viscosity;
     57        IssmDouble  mantle_density            = arguments->mantle_density;
     58        IssmDouble  lithosphere_thickness     = arguments->lithosphere_thickness;
     59        IssmDouble  rho_ice                   = arguments->rho_ice;
     60        int         disk_id                   = arguments->idisk;
     61        int         iedge                     = arguments->iedge;
     62        IssmDouble  yts                       = arguments->yts;
    11063
    11164        /*}}}*/
    11265
    11366        /*Modify inputs to match naruse code: */
    114         Ntime=numtimes;
    115         Ntimm=Ntime-1;
    116         Ntimp=Ntime+1;
     67        int Ntime=numtimes; // number of times with load history
     68        int Ntimm=Ntime-1; // Ntime-1 : for slope/y-cept of load segments
     69        int Ntimp=Ntime+1; // Ntime+1 : for evaluation time
    11770
    11871        /*Prepare block inputs for fortran distme and what0 routines of the naruse code: {{{*/
    11972        /*Now, let's set pset from the data that we got in input to GiaDeflectionCorex: */
    120         blockp_.pset[0]=lithosphere_thickness;
    121         blockp_.pset[1]=mantle_viscosity;
    122         blockp_.pset[2]=lithosphere_shear_modulus;
    123         blockp_.pset[3]=mantle_shear_modulus;
    124         blockp_.pset[4]=lithosphere_density;
    125         blockp_.pset[5]=mantle_density;
    126         blockp_.pset[6]=re;
    127         blocko_.rhoi=rho_ice;
     73        blockp_.pset[0]=reCast<IssmPDouble>(lithosphere_thickness);
     74        blockp_.pset[1]=reCast<IssmPDouble>(mantle_viscosity);
     75        blockp_.pset[2]=reCast<IssmPDouble>(lithosphere_shear_modulus);
     76        blockp_.pset[3]=reCast<IssmPDouble>(mantle_shear_modulus);
     77        blockp_.pset[4]=reCast<IssmPDouble>(lithosphere_density);
     78        blockp_.pset[5]=reCast<IssmPDouble>(mantle_density);
     79        blockp_.pset[6]=reCast<IssmPDouble>(re);
     80        blocko_.rhoi=reCast<IssmPDouble>(rho_ice);
     81   blockrad_.distrad=reCast<IssmPDouble>(ri)/1000.0; // in km
    12882
    12983        /*loading history: */
    130         blocky_zhload=xNew<IssmDouble>(Ntime);
    131         for(i=0;i<Ntime;i++){
    132         blocky_zhload[i]=hes[i];
     84        IssmPDouble* blocky_zhload=xNew<IssmPDouble>(Ntime);
     85        for(int i=0;i<Ntime;i++) blocky_zhload[i]=reCast<IssmPDouble>(hes[i]);
     86
     87        /*times in kyr: */
     88        IssmPDouble* blockt_time=xNew<IssmPDouble>(Ntimp);
     89        for(int i=0;i<Ntimp;i++){
     90                blockt_time[i]=times[i]/1000.0/yts;
     91                if(i==numtimes-1) blockt_time[i]=reCast<IssmPDouble>(times[numtimes-1])/1000.0/yts; // final loading time, same as evaluation time
     92                if(i==numtimes)   blockt_time[i]=reCast<IssmPDouble>(times[numtimes-1])/1000.0/yts;   // evaluation time
    13393        }
    13494
    135         /*times in kyr: */
    136         blockt_time=xNew<IssmDouble>(Ntimp);
    137         for (i=0;i<Ntimp;i++){
    138                 blockt_time[i]=times[i]/1000.0/yts;
    139                 if(i==numtimes-1)blockt_time[i]=times[numtimes-1]/1000.0/yts; // final loading time, same as evaluation time
    140                 if(i==numtimes)blockt_time[i]=times[numtimes-1]/1000.0/yts;   // evaluation time
    141         }
    142 
    143         /*bi: */
    144         blockt_bi=xNew<IssmDouble>(Ntimm);
    145 
    146         /*dmi: */
    147         blockt_dmi=xNew<IssmDouble>(Ntimm);
    148 
    149         /*radial distance of i-th element: */
    150         blockrad_.distrad=ri/1000.0; // in km
    151         /*}}}*/
     95        IssmPDouble* blockt_bi=xNew<IssmPDouble>(Ntimm);
     96        IssmPDouble* blockt_dmi=xNew<IssmPDouble>(Ntimm);
    15297
    15398        /*Call distme driver: */
    … …  
    158103
    159104        /*output solution: */
    160         wi = blocks_.aswokm_w;
    161         dwidt = blocks_.aswokm_dwdt;
    162         *pwi=wi;
    163         *pdwidt=dwidt;
     105        *pwi=reCast<IssmDouble>(blocks_.aswokm_w);
     106        *pdwidt=reCast<IssmDouble>(blocks_.aswokm_dwdt);
    164107
    165108        /*Free ressources: */
    166         xDelete<IssmDouble>(blockt_time);
    167         xDelete<IssmDouble>(blockt_bi);
    168         xDelete<IssmDouble>(blockt_dmi);
    169         xDelete<IssmDouble>(blocky_zhload);
     109        xDelete<IssmPDouble>(blockt_time);
     110        xDelete<IssmPDouble>(blockt_bi);
     111        xDelete<IssmPDouble>(blockt_dmi);
     112        xDelete<IssmPDouble>(blocky_zhload);
    170113
    171114}
Note: See TracChangeset for help on using the changeset viewer.