Changeset 24208 for issm/trunk-jpl/src


Ignore:
Timestamp:
10/04/19 15:17:14 (5 years ago)
Author:
wchu28
Message:

remove from Penta

Location:
issm/trunk-jpl/src/c/classes/Elements
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r24205 r24208  
    219219                virtual void       ComputeStressTensor(void)=0;
    220220                virtual void       ComputeEsaStrainAndVorticity(void)=0;
    221            virtual void          ComputeRadarAttenuation(void){_error_("not implemented yet");};
    222                 virtual void       ComputeRadarPower(void)=0;           
    223221                virtual void       Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
    224222                virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r24205 r24208  
    784784        /*Clean up and return*/
    785785        delete gauss;
    786 }
    787 /*}}}*/
    788 void       Penta::ComputeRadarAttenuation(){/*{{{*/
    789    IssmDouble  xyz_list[NUMVERTICES][3];
    790         IssmDouble  eps0=8.85418782e-12;
    791         IssmDouble  eps_ice=3.17;
    792         IssmDouble  e=2.7183;
    793         IssmDouble  c=299792458;
    794         IssmDouble  k=1.3806488e-23;
    795         IssmDouble  eVtoJ=1.6e-19;
    796         IssmDouble  Tr_M07=252.1500;
    797         IssmDouble  Tr_W97=258.1500;
    798         IssmDouble  sig_ice_M07=9.2;
    799         IssmDouble      sig_ice_W97=9;
    800         IssmDouble  cond_H_M07=3.2;
    801         IssmDouble  cond_H_W97=4;
    802         IssmDouble  cond_Cl_M07=0.43;
    803         IssmDouble  cond_Cl_W97=0.55;
    804         IssmDouble  cond_NH_M07=0.8;
    805         IssmDouble  cond_NH_W97=1;
    806         IssmDouble  E_M07=8.1600e-20;
    807         IssmDouble  E_W97=9.2800e-20;
    808         IssmDouble  E_H_M07=3.2000e-20;
    809         IssmDouble  E_H_W97=3.3600e-20;
    810         IssmDouble  E_Cl_M07=3.0400e-20;
    811         IssmDouble  E_Cl_W97=3.6800e-20;
    812         IssmDouble  E_NH=3.6800e-20;
    813         /*IssmDouble  mol_H_hol=1.6;
    814         IssmDouble  mol_Cl_hol=0.4;
    815         IssmDouble  mol_NH_hol=0.5;
    816         IssmDouble  mol_H_lgp=0.2;
    817         IssmDouble  mol_Cl_lgp=1.8;
    818         IssmDouble  mol_NH_lgp=0.4; */
    819         IssmDouble  radar_ice_period;
    820         IssmDouble  mol_H, mol_Cl, mol_NH;
    821         IssmDouble  attenuation_rate_macgregor[NUMVERTICES];
    822         IssmDouble  attenuation_rate_wolff[NUMVERTICES];
    823         IssmDouble  temperature, attenuation_rate_M07_pureice, attenuation_rate_M07_H, attenuation_rate_M07_Cl, attenuation_rate_M07_NH;
    824         IssmDouble      attenuation_rate_W97_pureice, attenuation_rate_W97_H, attenuation_rate_W97_Cl, attenuation_rate_W97_NH;
    825         IssmDouble m1, m2, m3, m4, m5, m6, m7, m8;
    826         IssmDouble  w2, w3, w4, w5, w6, w7, w8;
    827         GaussPenta* gauss=NULL;
    828 
    829    /* Get node coordinates and dof list: */
    830         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    831 
    832         /*Retrieve all inputs we will be needing: */
    833         Input* temp_input=inputs->GetInput(TemperatureEnum); _assert_(temp_input);
    834         //Input* ice_period_input=inputs->GetInput(RadarIcePeriodEnum); _assert_(ice_period_input);
    835        
    836         /* Start looping on the number of vertices: */
    837         gauss=new GaussPenta();
    838         for (int iv=0;iv<NUMVERTICES;iv++){
    839                 gauss->GaussVertex(iv);
    840 
    841    /*Get ice temperature: */
    842                 temp_input->GetInputValue(&temperature,gauss);
    843 //              ice_period_input->GetInputValue(&radar_ice_period,gauss);
    844 
    845         /*Ice period condition: */
    846         //      if (radar_ice_period>0.){
    847                         mol_H=1.6;
    848                         mol_Cl=0.4;
    849                         mol_NH=0.5;
    850         //      }
    851         //      else{
    852         //              mol_H=0.2;
    853         //              mol_Cl=1.8;
    854         //              mol_NH=0.4;
    855         //      }
    856 
    857         /*Compute M07 radar conductivity constant: */
    858                 m1=(10*log10(e))/(1000*eps0*sqrt(eps_ice)*c);
    859                 m2=E_M07/k;
    860                 m3=cond_H_M07*mol_H;
    861                 m4=E_H_M07/k;
    862                 m5=cond_Cl_M07*mol_Cl;
    863                 m6=E_Cl_M07/k;
    864                 m7=cond_NH_M07*mol_NH;
    865                 m8=E_NH/k;
    866 
    867         /*Compute MacGregor attenuation rate: */
    868                 attenuation_rate_M07_pureice=m1*sig_ice_M07*exp(m2*((1/Tr_M07)-(1/temperature)));
    869                 attenuation_rate_M07_H=m3*exp(m4*((1/Tr_M07)-(1/temperature)));
    870                 attenuation_rate_M07_Cl=m5*exp(m6*((1/Tr_M07)-(1/temperature)));
    871                 attenuation_rate_M07_NH=m7*exp(m8*((1/Tr_M07)-(1/temperature)));
    872                 attenuation_rate_macgregor[iv]=attenuation_rate_M07_pureice+attenuation_rate_M07_H+attenuation_rate_M07_Cl+attenuation_rate_M07_NH;
    873 
    874         /*Compute W97 radar conductivity constant: */
    875                 w2=E_W97/k;
    876                 w3=cond_H_W97*mol_H;
    877                 w4=E_H_W97/k;
    878                 w5=cond_Cl_W97*mol_Cl;
    879                 w6=E_Cl_W97/k;
    880                 w7=cond_NH_W97*mol_NH;
    881                 w8=E_NH/k;
    882 
    883         /*Compute Wolff attenuation rate: */
    884     attenuation_rate_W97_pureice=m1*sig_ice_W97*exp(w2*((1/Tr_W97)-(1/temperature)));
    885          attenuation_rate_W97_H=w3*exp(w4*((1/Tr_W97)-(1/temperature)));
    886          attenuation_rate_W97_Cl=w5*exp(w6*((1/Tr_W97)-(1/temperature)));
    887          attenuation_rate_W97_NH=w7*exp(w8*((1/Tr_W97)-(1/temperature)));
    888          attenuation_rate_wolff[iv]=attenuation_rate_W97_pureice+attenuation_rate_W97_H+attenuation_rate_W97_Cl+attenuation_rate_W97_NH;
    889         }
    890 
    891    /*Add Attenuation rate results into inputs*/
    892         this->inputs->AddInput(new PentaInput(RadarAttenuationMacGregorEnum,&attenuation_rate_macgregor[0],P1Enum));
    893         this->inputs->AddInput(new PentaInput(RadarAttenuationWolffEnum,&attenuation_rate_wolff[0],P1Enum));
    894 
    895    /*Clean up and return*/
    896         delete gauss;
    897 }
    898 /*}}}*/
    899 void       Penta::ComputeRadarPower(){/*{{{*/   
    900         IssmDouble  xyz_list[NUMVERTICES][3];
    901         IssmDouble  power_M07[NUMVERTICES];
    902         IssmDouble  power_W97[NUMVERTICES];
    903         IssmDouble      depth[NUMVERTICES];
    904         IssmDouble  aircraft_elev=0.5;
    905         IssmDouble  eps_ice=3.15;
    906         IssmDouble  t_tp=273.15;                        /* triple point temperature [K] */
    907         IssmDouble  p_tp=611.73;                        /* water pressure [Pa] */
    908         IssmDouble  gamma=7.4200e-07;   /* Clausius-Clapeyron constant [K/kPa] */
    909         IssmDouble  attenuation_rate_macgregor, attenuation_rate_wolff, attenuation_total_M07, attenuation_total_W97;
    910         IssmDouble  thickness, surface, z, temperature, geometric_loss, reflectivity;
    911         IssmDouble  rho_ice, gravity, pressure, pressure_melting_pt, frozen_temp, basal_temp, basal_pmp; 
    912         GaussPenta* gauss=NULL;
    913 
    914    /* Get node coordinates and dof list: */
    915         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    916 
    917         /*Retrieve all inputs we will be needing: */
    918    Input* atten_input_M07=inputs->GetInput(RadarAttenuationMacGregorEnum); _assert_(atten_input_M07);
    919         Input* atten_input_W97=inputs->GetInput(RadarAttenuationWolffEnum); _assert_(atten_input_W97);
    920         Input* surf_input=inputs->GetInput(SurfaceEnum); _assert_(surf_input);
    921         Input* thick_input=inputs->GetInput(ThicknessEnum); _assert_(thick_input);
    922         Input* temp_input=inputs->GetInput(TemperatureEnum); _assert_(temp_input);
    923 
    924         /* Start looping on the number of vertices: */
    925    gauss=new GaussPenta();
    926         for (int iv=0;iv<NUMVERTICES;iv++){
    927                 gauss->GaussVertex(iv);
    928        
    929         /*Get all the inputs: */
    930                 atten_input_M07->GetInputValue(&attenuation_rate_macgregor,gauss);
    931                 atten_input_W97->GetInputValue(&attenuation_rate_wolff,gauss);
    932                 thick_input->GetInputValue(&thickness,gauss);
    933                 temp_input->GetInputValue(&temperature,gauss);
    934                 surf_input->GetInputValue(&surface,gauss);
    935 
    936         /*Compute depth below the ice surface: */
    937                 z=xyz_list[iv][2];
    938                 depth[iv]=(surface-z)/1e3;
    939        
    940         /*Compute total attenuation: */
    941                 attenuation_total_M07=attenuation_rate_macgregor*depth[iv];
    942                 attenuation_total_W97=attenuation_rate_wolff*depth[iv];
    943        
    944         /*Compute geometric loss: */
    945                 geometric_loss=10*log10((depth[iv]+aircraft_elev)/sqrt(eps_ice));
    946        
    947         /*Compute radar power: */
    948                 power_M07[iv]=-geometric_loss-attenuation_total_M07;
    949                 power_W97[iv]=-geometric_loss-attenuation_total_W97;
    950        
    951                 /*Identify basal elements: */
    952                 if(this->IsOnBase() && iv<NUMVERTICES2D){
    953 
    954                         /*Compute pressure melting point: */
    955                         rho_ice=FindParam(MaterialsRhoIceEnum);
    956                         gravity=FindParam(ConstantsGEnum);
    957                         pressure=rho_ice*gravity*thickness;     
    958                         pressure_melting_pt=t_tp-gamma*(pressure-p_tp);
    959                         if((temperature-pressure_melting_pt)<=-1){
    960                                 reflectivity=-40;
    961                         }
    962                         else if((temperature-pressure_melting_pt)>-1 && (temperature-pressure_melting_pt)<0){
    963                                 reflectivity=0;
    964                         }
    965                         else{
    966                                 reflectivity=70;
    967                         }
    968                         power_M07[iv]=power_M07[iv]+reflectivity;
    969                         power_W97[iv]=power_W97[iv]+reflectivity;
    970                 }
    971         }
    972          
    973          /*Add power results into inputs*/
    974                 this->inputs->AddInput(new PentaInput(RadarPowerMacGregorEnum,&power_M07[0],P1Enum));
    975                 this->inputs->AddInput(new PentaInput(RadarPowerWolffEnum,&power_W97[0],P1Enum));
    976 
    977                 /*Clean up and return*/
    978                 delete gauss;
    979786}
    980787/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r24205 r24208  
    5858                void           ComputeSigmaNN(){_error_("not implemented yet");};
    5959                void           ComputeStressTensor();
    60                 void                            ComputeRadarAttenuation();
    61                 void           ComputeRadarPower();
    6260                void           Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    6361                void           ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r24205 r24208  
    4949                void        ComputeSigmaNN(){_error_("not implemented yet");};
    5050                void        ComputeStressTensor(){_error_("not implemented yet");};
    51                 void        ComputeRadarAttenuation(){_error_("not implemented yet");};
    52                 void        ComputeRadarPower(){_error_("not implemented yet");};
    5351                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
    5452                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r24205 r24208  
    4949                void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
    5050                void        ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");};
    51                 void        ComputeRadarAttenuation(){_error_("not implemented yet!");};
    52                 void        ComputeRadarPower(){_error_("not implemented yet");};
    5351                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
    5452                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r24205 r24208  
    6363                void        ComputeStressTensor();
    6464                void        ComputeSurfaceNormalVelocity();
    65                 void        ComputeRadarAttenuation(){_error_("not implemented yet");};
    66                 void        ComputeRadarPower(){_error_("not implemented yet");};               
    6765                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
    6866                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N, int M);
Note: See TracChangeset for help on using the changeset viewer.