Index: /issm/branches/trunk-dlcheng-ASE/src/c/Makefile.am =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/Makefile.am (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/Makefile.am (revision 27956) @@ -81,4 +81,5 @@ ./classes/Numberedcostfunction.cpp \ ./classes/Misfit.cpp \ + ./classes/MisfitAnnual.cpp \ ./classes/Cfsurfacesquare.cpp \ ./classes/Cfsurfacesquaretransient.cpp \ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Element.cpp (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Element.cpp (revision 27956) @@ -30,17 +30,4 @@ extern "C" void run_semic_(IssmDouble *sf_in, IssmDouble *rf_in, IssmDouble *swd_in, IssmDouble *lwd_in, IssmDouble *wind_in, IssmDouble *sp_in, IssmDouble *rhoa_in, IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *saccu_out, IssmDouble *smelt_out); - -extern "C" void run_semic_transient_(int *nx, int *ntime, int *nloop, - IssmDouble *sf_in, IssmDouble *rf_in, IssmDouble *swd_in, - IssmDouble *lwd_in, IssmDouble *wind_in, IssmDouble *sp_in, IssmDouble *rhoa_in, - IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_in, IssmDouble *qmr_in, - IssmDouble *tstic, - IssmDouble *hcrit, IssmDouble *rcrit, - IssmDouble *mask, IssmDouble *hice, IssmDouble *hsnow, - IssmDouble *albedo_in, IssmDouble *albedo_snow_in, - int *alb_scheme, IssmDouble *alb_smax, IssmDouble *alb_smin, IssmDouble *albi, IssmDouble *albl, - IssmDouble *Tamp, - IssmDouble *tmin, IssmDouble *tmax, IssmDouble *tmid, IssmDouble *mcrit, IssmDouble *wcrit, IssmDouble *tau_a, IssmDouble* tau_f, IssmDouble *afac, bool *verbose, - IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *smbi_out, IssmDouble *smbs_out, IssmDouble *saccu_out, IssmDouble *smelt_out, IssmDouble *refr_out, IssmDouble *albedo_out, IssmDouble *albedo_snow_out, IssmDouble *hsnow_out, IssmDouble *hice_out, IssmDouble *qmr_out); #endif // _HAVE_SEMIC_ @@ -57,4 +44,5 @@ this->parameters = NULL; this->element_type_list=NULL; + this->accumulator_values=NULL; }/*}}}*/ Element::~Element(){/*{{{*/ @@ -74,15 +62,11 @@ return false; }/*}}}*/ -void Element::ArmaProcess(bool isstepforarma,int arorder,int maorder,int numparams,int numbreaks,IssmDouble tstep_arma,IssmDouble* polyparams,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,IssmDouble* datebreaks,bool isfieldstochastic,int enum_type){/*{{{*/ +void Element::ArmaProcess(bool isstepforarma,int arorder,int maorder,IssmDouble telapsed_arma,IssmDouble tstep_arma,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,bool isfieldstochastic,int enum_type){/*{{{*/ const int numvertices = this->GetNumberOfVertices(); - int numperiods = numbreaks+1; - int basinid,M,N,arenum_type,maenum_type,basinenum_type,noiseenum_type,outenum_type,indperiod; - IssmDouble time,dt,starttime,noiseterm; + int basinid,M,N,arenum_type,maenum_type,basinenum_type,noiseenum_type,outenum_type; + IssmDouble time,dt,starttime,termconstant_basin,trend_basin,noiseterm; IssmDouble* arlagcoefs_basin = xNew(arorder); IssmDouble* malagcoefs_basin = xNew(maorder); - IssmDouble* datebreaks_basin = xNew(numbreaks); - IssmDouble* polyparams_basin = xNew(numperiods*numparams); IssmDouble* varlist = xNew(numvertices); - IssmDouble* sumpoly = xNewZeroInit(arorder+1); IssmDouble* valuesautoregression = NULL; IssmDouble* valuesmovingaverage = NULL; @@ -112,19 +96,5 @@ outenum_type = BasalforcingsSpatialDeepwaterMeltingRateEnum; break; - case(FrontalForcingsSubglacialDischargearmaEnum): - arenum_type = SubglacialdischargeValuesAutoregressionEnum; - maenum_type = SubglacialdischargeValuesMovingaverageEnum; - basinenum_type = FrontalForcingsBasinIdEnum; - noiseenum_type = SubglacialdischargeARMANoiseEnum; - outenum_type = FrontalForcingsSubglacialDischargeEnum; - break; - case(HydrologyarmapwEnum): - arenum_type = WaterPressureValuesAutoregressionEnum; - maenum_type = WaterPressureValuesMovingaverageEnum; - basinenum_type = HydrologyBasinsIdEnum; - noiseenum_type = FrictionWaterPressureNoiseEnum; - outenum_type = WaterPressureArmaPerturbationEnum; - break; - } + } /*Get time parameters*/ @@ -135,31 +105,8 @@ /*Get basin coefficients*/ this->GetInputValue(&basinid,basinenum_type); - for(int i=0;i0){ - for(int i=0;i0){ - /*Find index of tatstep compared to the breakpoints*/ - indperiod = 0; - for(int i=0;i=datebreaks_basin[i]) indperiod = i+1; - } - /*Compute polynomial with parameters of indperiod*/ - if(indperiod==0) telapsed_break = tatstep-starttime; - else telapsed_break = tatstep-datebreaks_basin[indperiod-1]; - for(int j=0;j(numvertices*arorder); IssmDouble* initvaluesmovingaverage = xNewZeroInit(numvertices*maorder); - for(int i=0;iinputs->SetArrayInput(arenum_type,this->lid,initvaluesautoregression,numvertices*arorder); this->inputs->SetArrayInput(maenum_type,this->lid,initvaluesmovingaverage,numvertices*maorder); @@ -195,17 +142,12 @@ /*Compute autoregressive term*/ IssmDouble autoregressionterm=0.; - for(int ii=0;ii(numvertices*arorder); @@ -228,13 +170,101 @@ xDelete(arlagcoefs_basin); xDelete(malagcoefs_basin); - xDelete(datebreaks_basin); - xDelete(polyparams_basin); - xDelete(sumpoly); xDelete(varlist); xDelete(valuesautoregression); xDelete(valuesmovingaverage); -} - -/*}}}*/ +}/*}}}*/ +void Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble tstep_ar,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* lagcoefs,bool isfieldstochastic,int enum_type){/*{{{*/ + + const int numvertices = this->GetNumberOfVertices(); + int basinid,M,N,arenum_type,basinenum_type,noiseenum_type,outenum_type; + IssmDouble time,dt,starttime,termconstant_basin,trend_basin,noiseterm; + IssmDouble* lagcoefs_basin = xNew(arorder); + IssmDouble* varlist = xNew(numvertices); + IssmDouble* valuesautoregression = NULL; + Input* noiseterm_input = NULL; + + /*Get field-specific enums*/ + switch(enum_type){ + case(SMBarmaEnum): + arenum_type = SmbValuesAutoregressionEnum; + basinenum_type = SmbBasinsIdEnum; + noiseenum_type = SmbARMANoiseEnum; + outenum_type = SmbMassBalanceEnum; + break; + case(FrontalForcingsRignotarmaEnum): + arenum_type = ThermalforcingValuesAutoregressionEnum; + basinenum_type = FrontalForcingsBasinIdEnum; + noiseenum_type = ThermalforcingARMANoiseEnum; + outenum_type = ThermalForcingEnum; + break; + case(BasalforcingsDeepwaterMeltingRatearmaEnum): + arenum_type = BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum; + basinenum_type = BasalforcingsLinearBasinIdEnum; + noiseenum_type = BasalforcingsDeepwaterMeltingRateNoiseEnum; + outenum_type = BasalforcingsSpatialDeepwaterMeltingRateEnum; + break; + } + + /*Get time parameters*/ + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + this->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + + /*Get basin coefficients*/ + this->GetInputValue(&basinid,basinenum_type); + for(int ii=0;ii(numvertices*arorder); + for(int i=0;iinputs->SetArrayInput(arenum_type,this->lid,initvaluesautoregression,numvertices*arorder); + xDelete(initvaluesautoregression); + } + + /*Get noise and autoregressive terms*/ + if(isfieldstochastic){ + noiseterm_input = this->GetInput(noiseenum_type); + Gauss* gauss = this->NewGauss(); + noiseterm_input->GetInputValue(&noiseterm,gauss); + delete gauss; + } + else noiseterm = 0.0; + this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M); + + /*If not AR model timestep: take the old values of variable*/ + if(isstepforar==false){ + for(int i=0;i(numvertices*arorder); + /*Assign newest values and shift older values*/ + for(int i=0;iinputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder); + xDelete(temparray); + } + + /*Add input to element*/ + this->AddInput(outenum_type,varlist,P1Enum); + + /*Cleanup*/ + xDelete(lagcoefs_basin); + xDelete(varlist); + xDelete(valuesautoregression); +}/*}}}*/ void Element::BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation){/*{{{*/ @@ -264,7 +294,4 @@ values[i]=deepwatermelt[i]*alpha+(1.-alpha)*upperwatermelt[basinid]; } - if(xIsNan(values[i])) _error_("NaN found in melt"); - if(xIsInf(values[i])) _error_("Inf found in melt"); - if(fabs(values[i])>1.e+10) _error_("melt exceeds 1.e+10"); } @@ -944,5 +971,5 @@ /*Get material parameters :*/ - IssmDouble rho_water = this->FindParam(MaterialsRhoFreshwaterEnum); + IssmDouble rho_water = this->FindParam(MaterialsRhoSeawaterEnum); IssmDouble rho_ice = this->FindParam(MaterialsRhoIceEnum); @@ -1903,11 +1930,8 @@ /*Are we in transient or static? */ if(M==1){ - if(N!=1) _error_("Size of Input "<SetElementInput(inputs,vector_enum,vector[0]); } - else if(M==iomodel->numberofvertices){ - if(N!=1) _error_("Size of Input "<SetElementInput(inputs,NUM_VERTICES,vertexlids,values,vector_enum); @@ -1948,5 +1972,5 @@ } else{ - _error_("Size of Input "<(evalues); @@ -1954,34 +1978,13 @@ } else{ - _error_("Size of Input "<SetElementInput(inputs,vector_enum,vector[0]); - } - else if(M==2){ - /*create transient input: */ - IssmDouble* times = xNew(N); - for(int t=0;tSetTransientInput(vector_enum,times,N); - TransientInput* transientinput = inputs->GetTransientInput(vector_enum); - - for(int t=0;tObjectEnum()){ - case TriaEnum: transientinput->AddTriaTimeInput( t,1,&(this->lid),&value,P0Enum); break; - case PentaEnum: transientinput->AddPentaTimeInput(t,1,&(this->lid),&value,P0Enum); break; - default: _error_("Not implemented yet"); - } - } - xDelete(times); - } - else if(M==iomodel->numberofelements){ - if(N!=1) _error_("Size of Input "<numberofelements){ if (code==5){ //boolean this->SetBoolInput(inputs,vector_enum,reCast(vector[this->Sid()])); @@ -2002,5 +2005,5 @@ TransientInput* transientinput = inputs->GetTransientInput(vector_enum); for(int t=0;tSid()+t]; + value=vector[N*this->Sid()+t]; switch(this->ObjectEnum()){ case TriaEnum: transientinput->AddTriaTimeInput( t,1,&(this->lid),&value,P0Enum); break; @@ -2011,8 +2014,24 @@ xDelete(times); } - - else{ - _error_("Size of Input "<(N); + if(M==1)times[0]=0; + if(M==2)for(int t=0;tSetTransientInput(vector_enum,times,N); + TransientInput* transientinput = inputs->GetTransientInput(vector_enum); + + for(int t=0;tObjectEnum()){ + case TriaEnum: transientinput->AddTriaTimeInput( t,1,&(this->lid),&value,P0Enum); break; + case PentaEnum: transientinput->AddPentaTimeInput(t,1,&(this->lid),&value,P0Enum); break; + default: _error_("Not implemented yet"); + } + } + xDelete(times); + } + else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); } else if(vector_type==3){ //Double array matrix @@ -2025,11 +2044,7 @@ xDelete(layers); } - else{ - _error_("Size of Input "<numberofelements){ - _error_("not implemented"); + if (code==5){ //boolean + _error_("not implemented"); + //datasetinput->AddInput(new BoolInput(input_enum,reCast(vector[this->Sid()])),input_id); + } + else if (code==6){ //integer + _error_("not implemented"); + //datasetinput->AddInput(new IntInput(input_enum,reCast(vector[this->Sid()])),input_id); + } + else if (code==7){ //IssmDouble + _error_("not implemented"); + //datasetinput->AddInput(new DoubleInput(input_enum,vector[this->Sid()]),input_id); + } + else _error_("could not recognize nature of vector from code " << code); } else if(M==iomodel->numberofelements+1){ _error_("not supported"); + ///*create transient input: */ + //IssmDouble* times = xNew(N); + //for(t=0;tSid()+t]; + // switch(this->ObjectEnum()){ + // case TriaEnum: transientinput->AddTimeInput(new TriaInput( input_enum,&value,P0Enum)); break; + // case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,&value,P0Enum)); break; + // case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,&value,P0Enum)); break; + // default: _error_("Not implemented yet"); + // } + //} + //xDelete(times); } else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long"); @@ -2357,17 +2400,4 @@ Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input); return (input->GetInputMax()<=0.); -} -/*}}}*/ -bool Element::IsAllMinThicknessInElement(){/*{{{*/ - - IssmDouble minthickness = this->FindParam(MasstransportMinThicknessEnum); - - Input* input=this->GetInput(ThicknessEnum); _assert_(input); - if(input->GetInputMax()<=(minthickness+0.00000001)){ - return true; - } - else{ - return false; - } } /*}}}*/ @@ -2441,23 +2471,12 @@ const int numvertices = this->GetNumberOfVertices(); bool isadjustsmb = false; - int basinid,bb1,bb2,mindex; - IssmDouble ela,refelevation_b,time,dt,fracyear,yts; - IssmDouble monthsteps[12] = {0.,1./12,2./12,3./12,4./12,5./12,6./12,7./12,8./12,9./12,10./12,11./12}; + int basinid,bb1,bb2; + IssmDouble ela,refelevation_b; IssmDouble* surfacelist = xNew(numvertices); IssmDouble* smblist = xNew(numvertices); - /* numelevbins values of lapse rates at current month */ + /* numelevbins values of lapse rates */ IssmDouble* lapserates_b = xNew(numelevbins); - /* (numelevbins-1) limits between elevation bins at current month (be cautious with indexing) */ + /* (numelevbins-1) limits between elevation bins (be cautious with indexing) */ IssmDouble* elevbins_b = xNew(numelevbins-1); - - /*Find month of current time step*/ - this->parameters->FindParam(&yts,ConstantsYtsEnum); - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - fracyear = time/yts-floor(time/yts); - for(int i=1;i<12;i++){ - if(fracyear>=monthsteps[i-1]) mindex = i-1; - } - if(fracyear>=monthsteps[11]) mindex = 11; /*Retrieve SMB values non-adjusted for SMB lapse rate*/ @@ -2468,24 +2487,21 @@ this->GetInputValue(&basinid,SmbBasinsIdEnum); refelevation_b = refelevation[basinid]; - /*Retrieve bins and laps rates for this basin at this month*/ - for(int ii=0;ii<(numelevbins-1);ii++) elevbins_b[ii] = elevbins[basinid*(numelevbins-1)*12+mindex*(numelevbins-1)+ii]; + for(int ii=0;ii<(numelevbins-1);ii++) elevbins_b[ii] = elevbins[basinid*(numelevbins-1)+ii]; for(int ii=0;ii=elevbins_b[ii]) bb1 = ii+1; - if(refelevation_b>=elevbins_b[ii]) bb2 = ii+1; + if(surfacelist[v]<=elevbins_b[ii]) bb1 = ii; + if(refelevation_b<=elevbins_b[ii]) bb2 = ii; } - + /*Check for elevations above highest bin limit */ + if(surfacelist[v]>elevbins_b[numelevbins-1-1]) bb1 = numelevbins-1; + if(refelevation_b>elevbins_b[numelevbins-1-1]) bb2 = numelevbins-1; /*Vertex and Reference elevation in same elevation bin*/ if(bb1==bb2){ @@ -2522,5 +2538,4 @@ IssmDouble deepwaterel,upperwaterel,deepwatermelt,upperwatermelt; IssmDouble base[MAXVERTICES]; - IssmDouble perturbation[MAXVERTICES]; IssmDouble values[MAXVERTICES]; IssmDouble time; @@ -2534,5 +2549,4 @@ this->GetInputListOnVertices(&base[0],BaseEnum); - this->GetInputListOnVertices(&perturbation[0],BasalforcingsPerturbationMeltingRateEnum); for(int i=0;i=upperwaterel){ @@ -2546,7 +2560,4 @@ values[i]=deepwatermelt*alpha+(1.-alpha)*upperwatermelt; } - - /*Add perturbation*/ - values[i] += perturbation[i]; } @@ -2766,98 +2777,13 @@ }/*}}}*/ -void Element::MonthlyFactorBasin(IssmDouble* monthlyfac,int enum_type){/*{{{*/ - - /*Variable declaration*/ - bool ratevariable; - const int numvertices = this->GetNumberOfVertices(); - int basinid,mindex,mindexnext,basinenum_type,varenum_type,indperiod; - IssmDouble time,dt,fracyear,fracyearnext,fracmonth,fracmonthnext,yts; - IssmDouble monthsteps[12] = {0.,1./12,2./12,3./12,4./12,5./12,6./12,7./12,8./12,9./12,10./12,11./12}; - IssmDouble* monthlyfac_b = xNew(12); - IssmDouble* monthlyrate_b = xNew(12); - IssmDouble* fracdtinmonth = xNew(12); - IssmDouble* rateinmonth = xNew(numvertices*12); - IssmDouble* varlistinput = xNew(numvertices); - IssmDouble* varlist = xNewZeroInit(numvertices); - - /*Get field-specific enums*/ - switch(enum_type){ - case(FrontalForcingsSubglacialDischargearmaEnum): - basinenum_type = FrontalForcingsBasinIdEnum; - varenum_type = FrontalForcingsSubglacialDischargeEnum; - ratevariable = true; - break; - case(HydrologyarmapwEnum): - basinenum_type = HydrologyBasinsIdEnum; - varenum_type = FrictionWaterPressureEnum; - ratevariable = false; - break; - } - - /*Evaluate the month index now and at (now-timestepjump)*/ - this->parameters->FindParam(&yts,ConstantsYtsEnum); - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); _assert_(dt=monthsteps[i-1]) mindex = i-1; - if(fracyearnext>=monthsteps[i-1]) mindexnext = i-1; - } - if(fracyear>=monthsteps[11]) mindex = 11; - if(fracyearnext>=monthsteps[11]) mindexnext = 11; - - /*Calculate fraction of the time step spent in each month*/ - for(int i=0;i<12;i++){ - if(mindexi) fracdtinmonth[i] = 1.0/dt*yts/12.0; - else if(mindexi && mindexnexti) fracdtinmonth[i] = 1.0/dt*yts/12.0; - else if(mindex>i && mindexnextGetInputValue(&basinid,basinenum_type); - for(int i=0;i<12;i++) monthlyfac_b[i] = monthlyfac[basinid*12+i]; - - /*Retrieve input*/ - this->GetInputListOnVertices(varlistinput,varenum_type); - - /*Calculate monthly rate for each month and weight-average it for application over dt*/ - for(int v=0;vAddInput(varenum_type,varlist,P1DGEnum); - - /*Clean-up*/ - xDelete(fracdtinmonth); - xDelete(rateinmonth); - xDelete(monthlyfac_b); - xDelete(monthlyrate_b); - xDelete(varlist); - xDelete(varlistinput); -}/*}}}*/ -void Element::MonthlyPiecewiseLinearEffectBasin(int nummonthbreaks,IssmDouble* monthlyintercepts,IssmDouble* monthlytrends,IssmDouble* monthlydatebreaks,int enum_type){/*{{{*/ +void Element::MonthlyEffectBasin(IssmDouble* monthlyeff, int enum_type){/*{{{*/ /*Variable declaration*/ const int numvertices = this->GetNumberOfVertices(); - int basinid,mindex,basinenum_type,varenum_type,indperiod; - int numperiods = nummonthbreaks+1; + int basinid,mindex,basinenum_type,varenum_type; IssmDouble time,fracyear,yts; IssmDouble monthsteps[12] = {0.,1./12,2./12,3./12,4./12,5./12,6./12,7./12,8./12,9./12,10./12,11./12}; - IssmDouble* datebreaks_b = xNew(nummonthbreaks); - IssmDouble* intercepts_b = xNew(numperiods*12); - IssmDouble* trends_b = xNew(numperiods*12); - IssmDouble* varlist = xNew(numvertices); + IssmDouble* monthlyeff_b = xNew(12); + IssmDouble* varlist = xNew(numvertices); /*Get field-specific enums*/ @@ -2880,26 +2806,7 @@ /*Get basin-specific parameters of the element*/ this->GetInputValue(&basinid,basinenum_type); - if(nummonthbreaks>0){ - for(int i=0;i0){ - /*Find index of time compared to the breakpoints*/ - indperiod = 0; - for(int i=0;idatebreaks_b[i]) indperiod = i+1; - } - /*Compute intercept+trend terms with parameters of indperiod*/ - if(indperiod==0) telapsed_break = time; - else telapsed_break = time-datebreaks_b[indperiod-1]; - piecewiselin = intercepts_b[indperiod]+trends_b[indperiod]*telapsed_break; - } - else piecewiselin = intercepts_b[indperiod]+trends_b[indperiod]*time; + for(int ii=0;ii<12;ii++){ + monthlyeff_b[ii] = monthlyeff[basinid*12+ii]; + } /*Retrieve values non-adjusted for monthly effects*/ @@ -2907,5 +2814,7 @@ /*Adjust values using monthly effects*/ - for(int v=0;v(datebreaks_b); - xDelete(intercepts_b); - xDelete(trends_b); + xDelete(monthlyeff_b); xDelete(varlist); -} -/*}}}*/ +}/*}}}*/ void Element::BeckmannGoosseFloatingiceMeltingRate(){/*{{{*/ @@ -3741,343 +3647,4 @@ } /*}}}*/ -void Element::SmbDebrisEvatt(){/*{{{*/ - - const int NUM_VERTICES = this->GetNumberOfVertices(); - const int NUM_VERTICES_DAYS_PER_YEAR = NUM_VERTICES * 365; // 365 FIXME - - int i,vertexlids[MAXVERTICES];; - IssmDouble* smb=xNew(NUM_VERTICES); - IssmDouble* melt=xNew(NUM_VERTICES); - IssmDouble* summermelt=xNew(NUM_VERTICES); - IssmDouble* albedo=xNew(NUM_VERTICES); - IssmDouble* summeralbedo=xNew(NUM_VERTICES); - IssmDouble* accu=xNew(NUM_VERTICES); - - // climate inputs - IssmDouble* temperature=xNew(NUM_VERTICES_DAYS_PER_YEAR); - IssmDouble* precip=xNew(NUM_VERTICES_DAYS_PER_YEAR); - IssmDouble* lw=xNew(NUM_VERTICES_DAYS_PER_YEAR); - IssmDouble* sw=xNew(NUM_VERTICES_DAYS_PER_YEAR); - IssmDouble* wind=xNew(NUM_VERTICES_DAYS_PER_YEAR); - IssmDouble* humidity=xNew(NUM_VERTICES_DAYS_PER_YEAR); - IssmDouble* yearlytemperatures=xNew(NUM_VERTICES); memset(yearlytemperatures, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* p_ampl=xNew(NUM_VERTICES); - IssmDouble* t_ampl=xNew(NUM_VERTICES); - IssmDouble* lw_ampl=xNew(NUM_VERTICES); - IssmDouble* sw_ampl=xNew(NUM_VERTICES); - IssmDouble* wind_ampl=xNew(NUM_VERTICES); - IssmDouble* humidity_ampl=xNew(NUM_VERTICES); - - IssmDouble* surface=xNew(NUM_VERTICES); - IssmDouble* s0t=xNew(NUM_VERTICES); - IssmDouble* snowheight=xNew(NUM_VERTICES); - IssmDouble* debriscover=xNew(NUM_VERTICES); - IssmDouble rho_water,rho_ice,Tf,debris,debris_here; - IssmDouble qlaps,rlaps,dsgrad,dlgrad,windspeedgrad,humiditygrad,Tm; - IssmDouble inv_twelve=1./365.; - IssmDouble time,yts,time_yr,lambda; - IssmDouble DailyMelt,CleanIceDailyMelt, CumDailyMelt=0,CleanIceMelt,CumDailySummerMelt=0; - IssmDouble MeanAlbedo=0, MeanSummerAlbedo=0; - bool isdebris,isAnderson,iscryokarst; - this->parameters->FindParam(&isdebris,TransientIsdebrisEnum); - this->parameters->FindParam(&isAnderson,SmbDebrisIsAndersonEnum); - this->parameters->FindParam(&iscryokarst,SmbDebrisIsCryokarstEnum); - IssmDouble PhiD=0.,p; - IssmDouble icealbedo=this->FindParam(SmbIcealbedoEnum); - IssmDouble snowalbedo=this->FindParam(SmbSnowalbedoEnum); - IssmDouble debrisalbedo=this->FindParam(SmbDebrisalbedoEnum); - IssmDouble Lm=this->FindParam(MaterialsLatentheatEnum); - IssmDouble D0=this->FindParam(SmbDebrisAndersonD0Enum); - int step; - this->FindParam(&step,StepEnum); - - // cryokarst - int dim=1,domaintype; - this->parameters->FindParam(&domaintype,DomainTypeEnum); - if(domaintype!=Domain2DverticalEnum){ - dim=2; - } - IssmDouble taud_plus=110e3, taud_minus=60e3; - IssmDouble taud, slope, gravity, taudx, taudy; - this->parameters->FindParam(&gravity,ConstantsGEnum); - IssmDouble* slopex = xNew(NUM_VERTICES); - IssmDouble* slopey = xNew(NUM_VERTICES); - IssmDouble* icethickness = xNew(NUM_VERTICES); - - /*Get material parameters :*/ - rho_water=this->FindParam(MaterialsRhoSeawaterEnum); - rho_ice=this->FindParam(MaterialsRhoIceEnum); - IssmDouble sconv=(rho_water/rho_ice); - Tf=this->FindParam(MaterialsMeltingpointEnum); - - /*Get parameters for height corrections*/ - qlaps=this->FindParam(SmbDesfacEnum); // comment MR; on alpine galciers we dont have the desertification effect - rlaps=this->FindParam(SmbRlapsEnum); - dsgrad=this->FindParam(SmbSWgradEnum); - dlgrad=this->FindParam(SmbLWgradEnum); - windspeedgrad=this->FindParam(SmbWindspeedgradEnum); - humiditygrad=this->FindParam(SmbHumiditygradEnum); - - /* Get time */ - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&yts,ConstantsYtsEnum); - time_yr=floor(time/yts)*yts; - - /*Get inputs*/ - DatasetInput* tempday =this->GetDatasetInput(SmbMonthlytemperaturesEnum); _assert_(tempday); - DatasetInput* precipday =this->GetDatasetInput(SmbPrecipitationEnum); _assert_(precipday); - DatasetInput* lwday =this->GetDatasetInput(SmbMonthlydlradiationEnum); _assert_(lwday); - DatasetInput* swday =this->GetDatasetInput(SmbMonthlydsradiationEnum); _assert_(swday); - DatasetInput* windday =this->GetDatasetInput(SmbMonthlywindspeedEnum); _assert_(windday); - DatasetInput* humidityday =this->GetDatasetInput(SmbMonthlyairhumidityEnum); _assert_(humidityday); - - /*loop over vertices: */ - Gauss* gauss=this->NewGauss(); - for(int month=0;month<365;month++){ - for(int iv=0;ivGaussVertex(iv); - tempday->GetInputValue(&temperature[iv*365+month],gauss,month); - temperature[iv*365+month]=temperature[iv*365+month]-Tf; // conversion from Kelvin to celcius for PDD module - precipday->GetInputValue(&precip[iv*365+month],gauss,month); - precip[iv*365+month]=precip[iv*365+month]*yts; // from m/s to m/a - lwday->GetInputValue(&lw[iv*365+month],gauss,month); - swday->GetInputValue(&sw[iv*365+month],gauss,month); - windday->GetInputValue(&wind[iv*365+month],gauss,month); - humidityday->GetInputValue(&humidity[iv*365+month],gauss,month); - } - } - - /*Recover info at the vertices: */ - GetInputListOnVertices(&surface[0],SurfaceEnum); - GetInputListOnVertices(&s0t[0],SmbS0tEnum); - GetInputListOnVertices(&snowheight[0],SmbSnowheightEnum); - GetInputListOnVertices(&debriscover[0],DebrisThicknessEnum); - GetInputListOnVertices(&t_ampl[0],SmbTemperaturesAnomalyEnum); - GetInputListOnVertices(&p_ampl[0],SmbPrecipitationsAnomalyEnum); - GetInputListOnVertices(&lw_ampl[0],SmbDsradiationAnomalyEnum); - GetInputListOnVertices(&sw_ampl[0],SmbDlradiationAnomalyEnum); - GetInputListOnVertices(&wind_ampl[0],SmbWindspeedAnomalyEnum); - GetInputListOnVertices(&humidity_ampl[0],SmbAirhumidityAnomalyEnum); - if(iscryokarst){ - GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum); - GetInputListOnVertices(&icethickness[0],ThicknessEnum); - if(dim==2){ - GetInputListOnVertices(&slopey[0],SurfaceSlopeYEnum); - } - taudx=rho_ice*gravity*icethickness[i]*slopex[i]; - if(dim==2) taudy=rho_ice*gravity*icethickness[i]*slopey[i]; - taud=sqrt(taudx*taudx+taudy*taudy); - } - IssmDouble Alphaeff,Alphaeff_cleanice; - - /*measure the surface mass balance*/ - for (int iv = 0; iv i.e - - /*Partition of precip in solid and liquid parts */ - IssmDouble temp_plus=1; - IssmDouble temp_minus=-1.; - IssmDouble frac_solid; - if(Tm>=temp_plus){ - frac_solid=0; - }else if(Tm<=temp_minus){ - frac_solid=1; - }else{ - frac_solid=1*(1-cos(PI*(temp_plus-Tm)/(temp_plus-temp_minus)))/2; - } - - /*Get yearly temperatures and accumulation */ - yearlytemperatures[iv]=yearlytemperatures[iv]+((temperature[iv*365+iday]-rlaps*st+Tf+t_ampl[iv]))/365; // Has to be in Kelvin - accu[iv]=accu[iv]+P*frac_solid; - if(yearlytemperatures[iv]>Tf) yearlytemperatures[iv]=Tf; - - CleanIceDailyMelt=((In-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ - Q*(1.-Alphaeff)+ - (Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))*Tm)/((1-PhiD)*rho_ice*Lm)/(1.+ - ((Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ - K*debris)-(Lv*Ustar*Ustar*((Humidity))*(exp(-Gamma*Xr)))/((1.-PhiD)* - rho_ice*Lm*Ustar)/(((Um - -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris))); - if(CleanIceDailyMelt<0) CleanIceDailyMelt=0.; - DailyMelt=CleanIceDailyMelt; - - if(ismb==1){ - - //snowheight[iv]=snowheight[iv]+(P-CleanIceDailyMelt*yts/365); - IssmDouble sn_prev; - sn_prev=snowheight[iv]; - snowheight[iv]=sn_prev+(-CleanIceDailyMelt*yts/365);//P - - if(snowheight[iv]<=0) snowheight[iv]=0.; - if(snowheight[iv]<=0.0001){ - p=debris_here*PhiD/(2*0.2*0.01); //Eq. 51 from Evatt et al 2015 without source term g*t - if(p>1.) p=1.; - if(p>=0.999){ - Alphaeff=debrisalbedo; - } else { - Alphaeff=Alphaeff_cleanice+p*(debrisalbedo-Alphaeff_cleanice); - } - debris=debris_here; - DailyMelt=((In-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ - Q*(1.-Alphaeff)+ - (Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))*Tm)/((1-PhiD)*rho_ice*Lm)/(1.+ - ((Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ - K*debris)-(Lv*Ustar*Ustar*((Humidity))*(exp(-Gamma*Xr)))/((1.-PhiD)* - rho_ice*Lm*Ustar)/(((Um-2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris))); - if(DailyMelt<0) DailyMelt=0.; - MeanSummerAlbedo=MeanSummerAlbedo+Alphaeff; - CumDailySummerMelt=CumDailySummerMelt+DailyMelt/365; - } - } - CumDailyMelt=CumDailyMelt+DailyMelt/365; - } - MeanAlbedo=MeanAlbedo+Alphaeff; - if(ismb==0) CleanIceMelt=CumDailyMelt; - } - - if(iscryokarst){ - if(taud>=taud_plus){ - lambda=0; - }else if(taud>=taud_minus & taudAddInput(SmbMassBalanceEnum,smb,P1Enum); - this->AddInput(SmbAccumulationEnum,accu,P1Enum); - this->AddInput(SmbMeltEnum,melt,P1Enum); - this->AddInput(SmbSummerMeltEnum,summermelt,P1Enum); - this->AddInput(SmbSnowheightEnum,snowheight,P1Enum); - this->AddInput(SmbAlbedoEnum,albedo,P1Enum); - this->AddInput(SmbSummerAlbedoEnum,summeralbedo,P1Enum); - this->AddInput(TemperaturePDDEnum,yearlytemperatures,P1Enum); // TemperaturePDD is wrong here, but don't want to create new Enum ... - - /*clean-up*/ - xDelete(temperature); - xDelete(precip); - xDelete(lw); - xDelete(sw); - xDelete(wind); - xDelete(humidity); - xDelete(smb); - xDelete(surface); - xDelete(melt); - xDelete(summermelt); - xDelete(albedo); - xDelete(summeralbedo); - xDelete(accu); - xDelete(yearlytemperatures); - xDelete(s0t); - xDelete(snowheight); - xDelete(debriscover); - xDelete(t_ampl); - xDelete(p_ampl); - xDelete(lw_ampl); - xDelete(sw_ampl); - xDelete(humidity_ampl); - xDelete(wind_ampl); - xDelete(slopex); - xDelete(slopey); - xDelete(icethickness); -} -/*}}}*/ - - - void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int* parray_size, int output_enum){/*{{{*/ @@ -4136,7 +3703,4 @@ this->CalvingRateVonmises(); break; - case CalvingVonmisesADEnum: - this->CalvingRateVonmisesAD(); - break; case CalvingCrevasseDepthEnum: this->CalvingCrevasseDepth(); @@ -4144,7 +3708,4 @@ case CalvingParameterizationEnum: this->CalvingRateParameterization(); - break; - case CalvingCalvingMIPEnum: - this->CalvingRateCalvingMIP(); break; case CalvingTestEnum: @@ -4313,12 +3874,12 @@ /* Start looping on the number of vertices: */ - Gauss* gauss=this->NewGauss(); + GaussTria gauss; for(int iv=0;ivGaussVertex(iv); + gauss.GaussVertex(iv); /* Get variables */ - bed_input->GetInputValue(&bed,gauss); - qsg_input->GetInputValue(&qsg,gauss); - TF_input->GetInputValue(&TF,gauss); + bed_input->GetInputValue(&bed,&gauss); + qsg_input->GetInputValue(&qsg,&gauss); + TF_input->GetInputValue(&TF,&gauss); if(basin_icefront_area[basinid]==0.) meltrates[iv]=0.; @@ -4329,5 +3890,5 @@ /* calculate melt rates */ meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s] - } + } if(xIsNan(meltrates[iv])) _error_("NaN found in vector"); @@ -4339,6 +3900,5 @@ /*Cleanup and return*/ - delete gauss; - xDelete(basin_icefront_area); + xDelete(basin_icefront_area); } /*}}}*/ @@ -4559,321 +4119,4 @@ xDelete(st); xDelete(s0gcm); -} -/*}}}*/ -void Element::SmbSemicTransient(){/*{{{*/ - - bool isverbose=VerboseSmb(); - if(isverbose && this->Sid()==0){ - _printf0_("smb core: initialize.\n"); - } - /*only compute SMB at the surface: */ - if (!IsOnSurface()) return; - - const int NUM_VERTICES = this->GetNumberOfVertices(); - - // daily forcing inputs - IssmDouble* dailyrainfall =xNew(NUM_VERTICES); - IssmDouble* dailysnowfall =xNew(NUM_VERTICES); - IssmDouble* dailydlradiation=xNew(NUM_VERTICES); - IssmDouble* dailydsradiation=xNew(NUM_VERTICES); - IssmDouble* dailywindspeed =xNew(NUM_VERTICES); - IssmDouble* dailypressure =xNew(NUM_VERTICES); - IssmDouble* dailyairdensity =xNew(NUM_VERTICES); - IssmDouble* dailyairhumidity=xNew(NUM_VERTICES); - IssmDouble* dailytemperature=xNew(NUM_VERTICES); - - // inputs: geometry - IssmDouble* s=xNew(NUM_VERTICES); - IssmDouble* s0gcm=xNew(NUM_VERTICES); - IssmDouble* st=xNew(NUM_VERTICES); - - // inputs - IssmDouble* tsurf_in =xNew(NUM_VERTICES); - IssmDouble* mask_in =xNew(NUM_VERTICES); - IssmDouble* Tamp_in =xNew(NUM_VERTICES); - IssmDouble* albedo_in =xNew(NUM_VERTICES); - IssmDouble* albedo_snow_in =xNew(NUM_VERTICES); - IssmDouble* hice_in =xNew(NUM_VERTICES); - IssmDouble* hsnow_in =xNew(NUM_VERTICES); - IssmDouble* qmr_in =xNew(NUM_VERTICES); - - // outputs - IssmDouble* tsurf_out =xNew(NUM_VERTICES); memset(tsurf_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* smb_out =xNew(NUM_VERTICES); memset(smb_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* smbi_out =xNew(NUM_VERTICES); memset(smb_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* smbs_out =xNew(NUM_VERTICES); memset(smb_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* saccu_out =xNew(NUM_VERTICES); memset(saccu_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* smelt_out =xNew(NUM_VERTICES); memset(smelt_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* refr_out =xNew(NUM_VERTICES); memset(refr_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* albedo_out =xNew(NUM_VERTICES); memset(albedo_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* albedo_snow_out =xNew(NUM_VERTICES); memset(albedo_snow_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* hsnow_out =xNew(NUM_VERTICES); memset(hsnow_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* hice_out =xNew(NUM_VERTICES); memset(hice_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - IssmDouble* qmr_out =xNew(NUM_VERTICES); memset(qmr_out, 0., NUM_VERTICES*sizeof(IssmDouble)); - - IssmDouble rho_water,rho_ice,desfac,desfacElev,rlaps,rdl; - IssmDouble alb_smax, alb_smin, albi, albl; - IssmDouble hcrit, rcrit; // parameters for ? and refreezing. - int alb_scheme; - // albedo parameters - slatter - IssmDouble tmin, tmax; - // albedo parameters - isba - IssmDouble mcrit, tau_a, tau_f, wcrit; - // albedo parameters - alex - IssmDouble tmid, afac; - - IssmDouble tstart, time,yts,time_yr,dt; - - /* Get time: */ - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - this->parameters->FindParam(&yts,ConstantsYtsEnum); - this->parameters->FindParam(&tstart,TimesteppingStartTimeEnum); - time_yr=floor(time/yts)*yts; - //dt = dt * yts; - - /*Get material parameters :*/ - rho_water=this->FindParam(MaterialsRhoFreshwaterEnum); - rho_ice=this->FindParam(MaterialsRhoIceEnum); - desfac=this->FindParam(SmbDesfacEnum); - desfacElev=this->FindParam(SmbDesfacElevEnum); - rlaps=this->FindParam(SmbRlapsEnum); - rdl=this->FindParam(SmbRdlEnum); - - this->FindParam(&alb_scheme,SmbAlbedoSchemeEnum); - this->FindParam(&hcrit,SmbSemicHcritEnum); - this->FindParam(&rcrit,SmbSemicRcritEnum); - alb_smax=this->FindParam(SmbAlbedoSnowMaxEnum); - alb_smin=this->FindParam(SmbAlbedoSnowMinEnum); - albi=this->FindParam(SmbAlbedoIceEnum); - albl=this->FindParam(SmbAlbedoLandEnum); - - // albedo parameters - this->FindParam(&tmid,SmbSemicTmidEnum); - this->FindParam(&tmin,SmbSemicTminEnum); - this->FindParam(&tmax,SmbSemicTmaxEnum); - this->FindParam(&mcrit,SmbSemicMcritEnum); - this->FindParam(&wcrit,SmbSemicWcritEnum); - this->FindParam(&tau_a,SmbSemicTauAEnum); - this->FindParam(&tau_f,SmbSemicTauFEnum); - this->FindParam(&afac,SmbSemicAfacEnum); - - /* Recover info at the vertices: */ - GetInputListOnVertices(&s[0],SurfaceEnum); - GetInputListOnVertices(&s0gcm[0],SmbS0gcmEnum); - - if(isverbose && this->Sid()==0){ - _printf0_("smb core: allocate inputs.\n"); - _printf0_("smb core: time_yr : " << time_yr/yts <<"\n"); - _printf0_("smb core: time : " << time <<"\n"); - _printf0_("smb core: dt : " << dt <<"\n"); - } - /* loop over vertices and days */ - Gauss* gauss=this->NewGauss(); - /* Retrieve inputs: */ - Input* dailysnowfall_input = this->GetInput(SmbDailysnowfallEnum,time); _assert_(dailysnowfall_input); - Input* dailyrainfall_input = this->GetInput(SmbDailyrainfallEnum,time); _assert_(dailyrainfall_input); - Input* dailydlradiation_input = this->GetInput(SmbDailydlradiationEnum,time); _assert_(dailydlradiation_input); - Input* dailydsradiation_input = this->GetInput(SmbDailydsradiationEnum,time); _assert_(dailydsradiation_input); - Input* dailywindspeed_input = this->GetInput(SmbDailywindspeedEnum,time); _assert_(dailywindspeed_input); - Input* dailypressure_input = this->GetInput(SmbDailypressureEnum,time); _assert_(dailypressure_input); - Input* dailyairdensity_input = this->GetInput(SmbDailyairdensityEnum,time); _assert_(dailyairdensity_input); - Input* dailyairhumidity_input = this->GetInput(SmbDailyairhumidityEnum,time); _assert_(dailyairhumidity_input); - Input* dailytemperature_input = this->GetInput(SmbDailytemperatureEnum,time); _assert_(dailytemperature_input); - - /*temporal Enum depending on time*/ - int enum_temp =TemperatureSEMICEnum; - int enum_hice =SmbHIceEnum; - int enum_hsnow =SmbHSnowEnum; - int enum_albedo =SmbAlbedoEnum; - int enum_albedo_snow=SmbAlbedoSnowEnum; - int enum_qmr =SmbSemicQmrEnum; - if (tstart+dt == time) { - /* Load inital value at first time step*/ - enum_temp=TemperatureEnum; - enum_hice=SmbHIceInitEnum; - enum_hsnow=SmbHSnowInitEnum; - enum_albedo=SmbAlbedoInitEnum; - enum_albedo_snow=SmbAlbedoSnowInitEnum; - enum_qmr =SmbSemicQmrInitEnum; - } - //if(isverbose && this->Sid()==0)_printf0_("smb core: assign temp.\n"); - Input* tsurf_input = this->GetInput(enum_temp); _assert_(tsurf_in); - //if(isverbose && this->Sid()==0)_printf0_("smb core: assign mask.\n"); - Input* mask_input = this->GetInput(SmbMaskEnum); _assert_(mask_input); - //if(isverbose && this->Sid()==0)_printf0_("smb core: assign Tamp.\n"); - Input* Tamp_input = this->GetInput(SmbTampEnum); _assert_(Tamp_input); - //if(isverbose && this->Sid()==0)_printf0_("smb core: assign albedo.\n"); - Input* albedo_input = this->GetInput(enum_albedo); _assert_(albedo_input); - Input* albedo_snow_input = this->GetInput(enum_albedo_snow); _assert_(albedo_snow_input); - Input* hice_input = this->GetInput(enum_hice); _assert_(hice_input); - Input* hsnow_input = this->GetInput(enum_hsnow); _assert_(hsnow_input); - Input* qmr_input = this->GetInput(enum_qmr); _assert_(qmr_input); - - if(isverbose && this->Sid()==0)_printf0_("smb core: assign inputs done....\n"); - for(int iv=0;ivGaussVertex(iv); - /* get forcing */ - dailyrainfall_input->GetInputValue(&dailyrainfall[iv],gauss); - dailysnowfall_input->GetInputValue(&dailysnowfall[iv],gauss); - dailydlradiation_input->GetInputValue(&dailydlradiation[iv],gauss); - dailydsradiation_input->GetInputValue(&dailydsradiation[iv],gauss); - dailywindspeed_input->GetInputValue(&dailywindspeed[iv],gauss); - dailypressure_input->GetInputValue(&dailypressure[iv],gauss); - dailyairdensity_input->GetInputValue(&dailyairdensity[iv],gauss); - dailyairhumidity_input->GetInputValue(&dailyairhumidity[iv],gauss); - dailytemperature_input->GetInputValue(&dailytemperature[iv],gauss); - tsurf_input->GetInputValue(&tsurf_in[iv],gauss); - - /* Get Albedo information */ - albedo_input->GetInputValue(&albedo_in[iv],gauss); - albedo_snow_input->GetInputValue(&albedo_snow_in[iv],gauss); - mask_input->GetInputValue(&mask_in[iv],gauss); - Tamp_input->GetInputValue(&Tamp_in[iv],gauss); - - hsnow_input->GetInputValue(&hsnow_in[iv],gauss); - hice_input->GetInputValue(&hice_in[iv],gauss); - qmr_input->GetInputValue(&qmr_in[iv],gauss); - - /* Surface temperature correction */ - st[iv]=(s[iv]-s0gcm[iv])/1000.; - dailytemperature[iv]=dailytemperature[iv]-rlaps *st[iv]; - - /* Precipitation correction (Vizcaino et al. 2010) */ - if (s0gcm[iv] < desfacElev) { - dailysnowfall[iv] = dailysnowfall[iv]*exp(desfac*(max(s[iv],desfacElev)-desfacElev)); - dailyrainfall[iv] = dailyrainfall[iv]*exp(desfac*(max(s[iv],desfacElev)-desfacElev)); - }else{ - dailysnowfall[iv] = dailysnowfall[iv]*exp(desfac*(max(s[iv],desfacElev)-s0gcm[iv])); - dailyrainfall[iv] = dailyrainfall[iv]*exp(desfac*(max(s[iv],desfacElev)-s0gcm[iv])); - } - - /* downward longwave radiation correction (Marty et al. 2002) */ - st[iv]=(s[iv]-s0gcm[iv])/1000.; - dailydlradiation[iv]=dailydlradiation[iv]+rdl*st[iv]; - } - if(isverbose && this->Sid()==0){ - _printf0_("smb core: assign tsurf_in :" << tsurf_in[0] << "\n"); - _printf0_("smb core: assign dailytemperature:" << dailytemperature[0] << "\n"); - _printf0_("smb core: assign hsnow :" << hsnow_in[0] << "\n"); - _printf0_("smb core: assign hice :" << hice_in[0] << "\n"); - _printf0_("smb core: assign mask :" << mask_in[0] << "\n"); - _printf0_("smb core: assign Tamp :" << Tamp_in[0] << "\n"); - _printf0_("smb core: assign albedo :" << albedo_in[0] << "\n"); - _printf0_("smb core: assign albedo_snow :" << albedo_snow_in[0] << "\n"); - _printf0_("smb core: assign albedo_scheme :" << alb_scheme << "\n"); - _printf0_("smb core: assign qmr :" << qmr_in[0] << "\n"); - } - - if(isverbose && this->Sid()==0)_printf0_("smb core: call run_semic_transient module.\n"); - /* call semic */ - int nx=NUM_VERTICES, ntime=1, nloop=1; - bool semic_verbose=false; //VerboseSmb(); - run_semic_transient_(&nx, &ntime, &nloop, - dailysnowfall, dailyrainfall, dailydsradiation, dailydlradiation, - dailywindspeed, dailypressure, dailyairdensity, dailyairhumidity, dailytemperature, tsurf_in, qmr_in, - &dt, - &hcrit, &rcrit, - mask_in, hice_in, hsnow_in, - albedo_in, albedo_snow_in, - &alb_scheme, &alb_smax, &alb_smin, &albi, &albl, - Tamp_in, - &tmin, &tmax, &tmid, &mcrit, &wcrit, &tau_a, &tau_f, &afac, &semic_verbose, - tsurf_out, smb_out, smbi_out, smbs_out, saccu_out, smelt_out, refr_out, albedo_out, albedo_snow_out, hsnow_out, hice_out, qmr_out); - - for (int iv = 0; iv ice - w.e. : water equivalenet. - */ - smb_out[iv] = smb_out[iv]*yts; // w.e. m/sec -> m/yr - smbi_out[iv] = smbi_out[iv]*rho_water/rho_ice; // w.e. m/sec -> ice m/yr - smbs_out[iv] = smbs_out[iv]*yts; // w.e. m/sec -> m/yr - saccu_out[iv] = saccu_out[iv]*yts; // w.e. m/sec -> m/yr - smelt_out[iv] = smelt_out[iv]*rho_water/rho_ice; // w.e. m/sec -> ice m/yr - refr_out[iv] = refr_out[iv]*rho_water/rho_ice; // w.e. m/sec -> ice m/yr - } - - if(isverbose && this->Sid()==0){ - _printf0_("smb core: tsurf_out " << tsurf_out[0] << " " << tsurf_out[1] << " " << tsurf_out[2] << "\n"); - _printf0_("smb core: hice_out " << hice_out[0] << " " << hice_out[1] << " " << hice_out[2] << "\n"); - _printf0_("smb core: hsnow_out " << hsnow_out[0] << "\n"); - _printf0_("smb core: smb_ice " << smbi_out[0]*yts << "\n"); - _printf0_("smb core: smb_ice " << albedo_out[0] <<" "<ObjectEnum()){ - case TriaEnum: - this->AddInput(TemperatureSEMICEnum, &tsurf_out[0],P1DGEnum); - // SMBout = SMB_ice + SMB_snow values. - //this->AddInput(SmbMassBalanceTotalEnum,&smb_out[0],P1DGEnum); - // water equivalent SMB ice to ice equivalent. - this->AddInput(SmbMassBalanceEnum, &smbi_out[0],P1DGEnum); - this->AddInput(SmbMassBalanceIceEnum, &smbi_out[0],P1DGEnum); - this->AddInput(SmbMassBalanceSnowEnum, &smbs_out[0],P1DGEnum); - this->AddInput(SmbMassBalanceSemicEnum,&smb_out[0],P1DGEnum); - //this->AddInput(SmbMassBalanceSnowEnum,&smbs_out[0],P1DGEnum); - // saccu - accumulation of snow. - this->AddInput(SmbAccumulationEnum,&saccu_out[0],P1DGEnum); - // smelt - this->AddInput(SmbMeltEnum, &smelt_out[0],P1DGEnum); - this->AddInput(SmbRefreezeEnum, &refr_out[0],P1DGEnum); - this->AddInput(SmbAlbedoEnum, &albedo_out[0],P1DGEnum); - this->AddInput(SmbAlbedoSnowEnum, &albedo_snow_out[0],P1DGEnum); - this->AddInput(SmbHSnowEnum, &hsnow_out[0],P1DGEnum); - this->AddInput(SmbHIceEnum, &hice_out[0],P1DGEnum); - this->AddInput(SmbSemicQmrEnum, &qmr_out[0],P1DGEnum); - break; - case PentaEnum: - // TODO - break; - case TetraEnum: - // TODO - break; - default: _error_("Not implemented yet"); - } - - /*clean-up {{{*/ - delete gauss; - xDelete(dailysnowfall); - xDelete(dailyrainfall); - xDelete(dailydlradiation); - xDelete(dailydsradiation); - xDelete(dailywindspeed); - xDelete(dailypressure); - xDelete(dailyairdensity); - xDelete(dailyairhumidity); - xDelete(dailypressure); - xDelete(dailytemperature); - - /*for outputs*/ - xDelete(tsurf_out); - xDelete(smb_out); - xDelete(smbi_out); - xDelete(smbs_out); - xDelete(saccu_out); - xDelete(smelt_out); - xDelete(refr_out); - xDelete(albedo_out); - xDelete(albedo_snow_out); - xDelete(hsnow_out); - xDelete(hice_out); - xDelete(qmr_out); - - /*for inputs*/ - xDelete(hsnow_in); - xDelete(hice_in); - xDelete(mask_in); - xDelete(Tamp_in); - xDelete(albedo_in); - xDelete(albedo_snow_in); - xDelete(tsurf_in); - xDelete(qmr_in); - - /* for inputs:geometry */ - xDelete(s); - xDelete(st); - xDelete(s0gcm); - /*}}}*/ } /*}}}*/ @@ -5420,32 +4663,4 @@ } /*}}}*/ -void Element::SubglacialWaterPressure(int output_enum){/*{{{*/ - - bool ispwHydroArma; - int M; - int numvertices = this->GetNumberOfVertices(); - IssmDouble p_water[numvertices]; - IssmDouble* perturbationvalues = xNew(numvertices); - Gauss* gauss=this->NewGauss(); - Friction* friction = new Friction(this); - /*Calculate subglacial water pressure*/ - for(int i=0;iGaussVertex(i); - p_water[i] = friction->SubglacialWaterPressure(gauss); - } - /*Add perturbation in water pressure if HydrologyIsWaterPressureArmaEnum is true*/ - this->parameters->FindParam(&ispwHydroArma,HydrologyIsWaterPressureArmaEnum); - if(ispwHydroArma){ - this->GetInputListOnVertices(perturbationvalues,WaterPressureArmaPerturbationEnum); - for(int i=0;iAddInput(output_enum,p_water,P1DGEnum); - /*Clean-up*/ - delete gauss; - delete friction; - xDelete(perturbationvalues); - -}/*}}}*/ void Element::StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ @@ -5755,30 +4970,4 @@ } /*}}}*/ -IssmDouble Element::TotalSmbMelt(IssmDouble* mask, bool scaled){/*{{{*/ - - /*Retrieve values of the mask defining the element: */ - for(int i=0;iGetNumberOfVertices();i++){ - if(mask[this->vertices[i]->Sid()]<=0.){ - return 0.; - } - } - - /*Return: */ - return this->TotalSmbMelt(scaled); -} -/*}}}*/ -IssmDouble Element::TotalSmbRefreeze(IssmDouble* mask, bool scaled){/*{{{*/ - - /*Retrieve values of the mask defining the element: */ - for(int i=0;iGetNumberOfVertices();i++){ - if(mask[this->vertices[i]->Sid()]<=0.){ - return 0.; - } - } - - /*Return: */ - return this->TotalSmbRefreeze(scaled); -} -/*}}}*/ void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Element.h =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Element.h (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Element.h (revision 27956) @@ -57,4 +57,5 @@ int* element_type_list; + IssmDouble* accumulator_values; int element_type; @@ -68,6 +69,6 @@ /*bool AnyActive(void);*/ bool AnyFSet(void); - void ArmaProcess_pre18Oct2022(bool isstepforarma,int arorder,int maorder,IssmDouble telapsed_arma,IssmDouble tstep_arma,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,bool isfieldstochastic,int enum_type); - void ArmaProcess(bool isstepforarma,int arorder,int maorder,int numparams,int numbreaks,IssmDouble tstep_arma,IssmDouble* polyparams,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,IssmDouble* datebreaks,bool isfieldstochastic,int enum_type); + void ArmaProcess(bool isstepforarma,int arorder,int maorder,IssmDouble telapsed_arma,IssmDouble tstep_arma,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,bool isfieldstochastic,int enum_type); + void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble tstep_ar,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* lagcoefs,bool isfieldstochastic,int enum_type); void BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation); void CalvingSetZeroRate(void); @@ -140,5 +141,5 @@ void InputCreateP1FromConstant(Inputs* inputs,IoModel* iomodel,IssmDouble value,int vector_enum); void ControlInputCreate(IssmDouble* doublearray,IssmDouble* independents_min,IssmDouble* independents_max,Inputs*inputs,IoModel* iomodel,int M,int N,IssmDouble scale,int input_enum,int id); - void DatasetInputAdd(int enum_type,IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int input_enum); + void DatasetInputAdd(int enum_type,IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code,int input_enum); void InputUpdateFromConstant(IssmDouble constant, int name); void InputUpdateFromConstant(int constant, int name); @@ -155,5 +156,4 @@ bool IsOceanInElement(); bool IsOceanOnlyInElement(); - bool IsAllMinThicknessInElement(); bool IsLandInElement(); void Ismip6FloatingiceMeltingRate(); @@ -165,6 +165,5 @@ void MigrateGroundingLine(IssmDouble* sheet_ungrounding); void MismipFloatingiceMeltingRate(); - void MonthlyFactorBasin(IssmDouble* monthlyfac,int enum_type); - void MonthlyPiecewiseLinearEffectBasin(int nummonthbreaks,IssmDouble* monthlyintercepts,IssmDouble* monthlytrends,IssmDouble* monthlydatebreaks,int enum_type); + void MonthlyEffectBasin(IssmDouble* monthlyeff, int enum_type); void BeckmannGoosseFloatingiceMeltingRate(); void MungsmtpParameterization(void); @@ -177,5 +176,4 @@ void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac); void PositiveDegreeDaySicopolis(bool isfirnwarming); - void SmbDebrisEvatt(); void RignotMeltParameterization(); void ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum); @@ -188,5 +186,4 @@ void SetIntInput(Inputs* inputs,int enum_in,int value); void SmbSemic(); - void SmbSemicTransient(); int Sid(); void SmbGemb(IssmDouble timeinputs, int count); @@ -199,10 +196,7 @@ void StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input); void StressMaxPrincipalCreateInput(void); - void SubglacialWaterPressure(int output_enum); IssmDouble TotalFloatingBmb(IssmDouble* mask, bool scaled); IssmDouble TotalGroundedBmb(IssmDouble* mask, bool scaled); IssmDouble TotalSmb(IssmDouble* mask, bool scaled); - IssmDouble TotalSmbMelt(IssmDouble* mask, bool scaled); - IssmDouble TotalSmbRefreeze(IssmDouble* mask, bool scaled); void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int cs_enum); void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum); @@ -241,7 +235,5 @@ virtual void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");}; virtual void CalvingRateParameterization(void){_error_("not implemented yet");}; - virtual void CalvingRateCalvingMIP(void){_error_("not implemented yet");}; virtual void CalvingRateVonmises(void){_error_("not implemented yet");}; - virtual void CalvingRateVonmisesAD(void){_error_("not implemented yet");}; virtual void CalvingRateTest(void){_error_("not implemented yet");}; virtual void CalvingCrevasseDepth(void){_error_("not implemented yet");}; @@ -332,4 +324,6 @@ virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0; virtual IssmDouble Misfit(int modelenum,int observationenum,int weightsenum)=0; + virtual void MisfitAccumulate(int modelenum,IssmDouble dt)=0; + virtual IssmDouble MisfitAnnual(int modelenum,int observationenum,int weightsenum,IssmDouble annual_dt)=0; virtual IssmDouble MisfitArea(int weightsenum)=0; virtual void MovingFrontalVelocity(void){_error_("not implemented yet");}; @@ -390,6 +384,4 @@ virtual IssmDouble TotalGroundedBmb(bool scaled)=0; virtual IssmDouble TotalSmb(bool scaled)=0; - virtual IssmDouble TotalSmbMelt(bool scaled)=0; - virtual IssmDouble TotalSmbRefreeze(bool scaled)=0; virtual void Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0; virtual void UpdateConstraintsExtrudeFromBase(void)=0; @@ -417,5 +409,5 @@ virtual void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom)=0; virtual void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0; - virtual void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount)=0; + virtual void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids)=0; virtual void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0; virtual void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0; Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Penta.cpp (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Penta.cpp (revision 27956) @@ -177,13 +177,5 @@ } } - else if(interpolation_enum==P0Enum){ - Penta* penta=this; - for(;;){ - penta->AddInput(input_enum,&values[0],interpolation_enum); - if (penta->IsOnSurface()) break; - penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id); - } - } - else _error_("Interpolation "<GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=this->GetInput(VzEnum); _assert_(vz_input); + Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=this->GetInput(VzEnum); _assert_(vz_input); /* Start looping on the number of vertices: */ @@ -1451,15 +1443,103 @@ IssmDouble Penta::GetIcefrontArea(){/*{{{*/ - /*We need to be on base and cross the levelset*/ + IssmDouble bed[NUMVERTICES]; //basinId[NUMVERTICES]; + IssmDouble Haverage,frontarea; + IssmDouble x1,y1,x2,y2,distance; + IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES]; + int* indices=NULL; + IssmDouble* H=NULL;; + int nrfrontbed,numiceverts; + + if(!this->IsOnBase()) return 0; if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0; - if(!this->IsOnBase()) return 0; - - /*Spawn Tria element from the base of the Penta: */ - Tria* tria=(Tria*)SpawnTria(0,1,2); - IssmDouble frontarea = tria->GetIcefrontArea(); - delete tria->material; delete tria; - + + /*Retrieve all inputs and parameters*/ + Element::GetInputListOnVertices(&bed[0],BedEnum); + Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum); + Element::GetInputListOnVertices(&bases[0],BaseEnum); + Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum); + + nrfrontbed=0; + for(int i=0;iGetLevelsetIntersectionBase(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.); + _assert_(numiceverts); + + /*3 Write coordinates*/ + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + int counter = 0; + if((numiceverts>0) && (numiceverts(numthk); + for(int iv=0;iv(indices); + xDelete(H); + + _assert_(frontarea>0); return frontarea; -}/*}}}*/ +} +/*}}}*/ void Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ @@ -3411,4 +3491,7 @@ } + /*Get out if this is not an element input*/ + if(!IsInputEnum(control_enum)) return; + /*Prepare index list*/ ElementInput* input=this->inputs->GetControlInputData(control_enum,"value"); _assert_(input); @@ -4270,140 +4353,4 @@ /*Return: */ return Total_Smb; -} -/*}}}*/ -IssmDouble Penta::TotalSmbMelt(bool scaled){/*{{{*/ - - /*The smbmelt[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */ - IssmDouble base,smbmelt,rho_ice,scalefactor; - IssmDouble Total_SmbMelt=0; - IssmDouble lsf[NUMVERTICES]; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Get material parameters :*/ - rho_ice=FindParam(MaterialsRhoIceEnum); - - if(!IsIceInElement() || !IsOnSurface()) return 0.; - - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*First calculate the area of the base (cross section triangle) - * http://en.wikipedia.org/wiki/Triangle - * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ - base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); - - /*Now get the average SMB over the element*/ - Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum); - if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){ - /*Partially ice-covered element*/ - bool mainlyice; - int point; - IssmDouble* smbmelt_vertices = xNew(NUMVERTICES); - IssmDouble weights[NUMVERTICES2D]; - IssmDouble lsf2d[NUMVERTICES2D]; - IssmDouble f1,f2,phi; - Element::GetInputListOnVertices(&smbmelt_vertices[0],SmbMeltEnum); - for(int i=0;i(NUMVERTICES); - Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum); - /*Compute loop only over lower vertices: i(scalefactor_vertices); - } - else scalefactor = 1.0; - - /*Cleanup*/ - xDelete(smbmelt_vertices); - } - - else{ - /*Fully ice-covered element*/ - Input* smbmelt_input = this->GetInput(SmbMeltEnum); _assert_(smbmelt_input); - smbmelt_input->GetInputAverage(&smbmelt); - - if(scaled==true){ - Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input); - scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element - } - else scalefactor=1.0; - } - - Total_SmbMelt=rho_ice*base*smbmelt*scalefactor;// smbmelt on element in kg s-1 - - /*Return: */ - return Total_SmbMelt; -} -/*}}}*/ -IssmDouble Penta::TotalSmbRefreeze(bool scaled){/*{{{*/ - - /*The smbrefreeze[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */ - IssmDouble base,smbrefreeze,rho_ice,scalefactor; - IssmDouble Total_SmbRefreeze=0; - IssmDouble lsf[NUMVERTICES]; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Get material parameters :*/ - rho_ice=FindParam(MaterialsRhoIceEnum); - - if(!IsIceInElement() || !IsOnSurface()) return 0.; - - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*First calculate the area of the base (cross section triangle) - * http://en.wikipedia.org/wiki/Triangle - * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ - base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); - - /*Now get the average SMB over the element*/ - Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum); - if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){ - /*Partially ice-covered element*/ - bool mainlyice; - int point; - IssmDouble* smbrefreeze_vertices = xNew(NUMVERTICES); - IssmDouble weights[NUMVERTICES2D]; - IssmDouble lsf2d[NUMVERTICES2D]; - IssmDouble f1,f2,phi; - Element::GetInputListOnVertices(&smbrefreeze_vertices[0],SmbRefreezeEnum); - for(int i=0;i(NUMVERTICES); - Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum); - /*Compute loop only over lower vertices: i(scalefactor_vertices); - } - else scalefactor = 1.0; - - /*Cleanup*/ - xDelete(smbrefreeze_vertices); - } - - else{ - /*Fully ice-covered element*/ - Input* smbrefreeze_input = this->GetInput(SmbRefreezeEnum); _assert_(smbrefreeze_input); - smbrefreeze_input->GetInputAverage(&smbrefreeze); - - if(scaled==true){ - Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input); - scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element - } - else scalefactor=1.0; - } - - Total_SmbRefreeze=rho_ice*base*smbrefreeze*scalefactor;// smbrefreeze on element in kg s-1 - - /*Return: */ - return Total_SmbRefreeze; } /*}}}*/ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Penta.h =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Penta.h (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Penta.h (revision 27956) @@ -64,5 +64,4 @@ void ComputeSigmaNN(){_error_("not implemented yet");}; void ComputeStressTensor(); - //void ComputeMeanEla(IssmDouble* paltitude, int* pcounter); void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters,Inputs* inputsin); void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp); @@ -138,5 +137,7 @@ IssmDouble MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id); IssmDouble MinEdgeLength(IssmDouble* xyz_list); - IssmDouble Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");}; + IssmDouble Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");}; + void MisfitAccumulate(int modelenum,IssmDouble dt){_error_("not implemented yet");}; + IssmDouble MisfitAnnual(int modelenum,int observationenum,int weightsenum,IssmDouble annual_dt){_error_("not implemented yet");}; IssmDouble MisfitArea(int weightsenum){_error_("not implemented yet");}; void MovingFrontalVelocity(void); @@ -201,6 +202,4 @@ IssmDouble TotalGroundedBmb(bool scaled); IssmDouble TotalSmb(bool scaled); - IssmDouble TotalSmbMelt(bool scaled); - IssmDouble TotalSmbRefreeze(bool scaled); void Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement); void UpdateConstraintsExtrudeFromBase(void); @@ -230,5 +229,5 @@ void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");}; void SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");}; - void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids,int* vcount){_error_("not implemented yet");}; + void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");}; void SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");}; void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");}; Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/PentaRef.cpp (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/PentaRef.cpp (revision 27956) @@ -307,5 +307,4 @@ basis[6]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+zeta)*(1.-zeta); return; - #ifndef _HAVE_AD_ case P2xP1Enum: /*Corner nodes*/ @@ -460,5 +459,4 @@ basis[14]=gauss->coord3*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.); return; - #endif default: _error_("Element type "<coord1*gauss->coord2*gauss->coord3*zeta; return; - #ifndef _HAVE_AD_ case P2xP1Enum: /*Nodal function 1*/ @@ -1060,5 +1057,4 @@ return; - #endif default: _error_("Element type "<vertices = NULL; this->material = NULL; + this->accumulator_values = xNew(NUMVERTICES); + for(int i=0;iaccumulator_values[i] = 0; if(nummodels>0){ this->element_type_list=xNew(nummodels); @@ -105,4 +107,5 @@ else tria->element_type_list = NULL; tria->element_type=this->element_type; + tria->accumulator_values=this->accumulator_values; tria->numanalyses=nanalyses; @@ -345,68 +348,4 @@ smax_fl_input->GetInputValue(&sigma_max_floating,&gauss); smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss); - sl_input->GetInputValue(&sealevel,&gauss); - - /*Tensile stress threshold*/ - if(groundedice<0) - sigma_max = sigma_max_floating; - else - sigma_max = sigma_max_grounded; - - /*Assign values*/ - if(bed>sealevel){ - calvingrate[iv] = 0.; - } - else{ - calvingrate[iv] = sqrt(vx*vx+vy*vy)*sigma_vm/sigma_max; - } - } - - /*Add input*/ - this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - this->CalvingRateToVector(); -} -/*}}}*/ -void Tria::CalvingRateVonmisesAD(){/*{{{*/ - - /*First, compute Von Mises Stress*/ - this->ComputeSigmaVM(); - - /*Now compute calving rate*/ - IssmDouble calvingrate[NUMVERTICES]; - IssmDouble sigma_vm,vx,vy; - IssmDouble sigma_max,sigma_max_floating,sigma_max_grounded,n; - IssmDouble groundedice,bed,sealevel; - int M; - int basinid; - IssmDouble* sigma_max_floating_basin=NULL; - IssmDouble* sigma_max_grounded_basin=NULL; - - /*Retrieve all inputs and parameters we will need*/ - Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input); - Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); - Input* bs_input = this->GetInput(BaseEnum); _assert_(bs_input); - Input* sl_input = this->GetInput(SealevelEnum); _assert_(sl_input); - Input* sigma_vm_input = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input); - - this->Element::GetInputValue(&basinid,CalvingBasinIdEnum); - - parameters->FindParam(&sigma_max_floating_basin,&M,CalvingADStressThresholdFloatingiceEnum); - parameters->FindParam(&sigma_max_grounded_basin,&M,CalvingADStressThresholdGroundediceEnum); - - sigma_max_floating = sigma_max_floating_basin[basinid]; - sigma_max_grounded = sigma_max_grounded_basin[basinid]; - - /* Start looping on the number of vertices: */ - GaussTria gauss; - for(int iv=0;ivGetInputValue(&sigma_vm,&gauss); - vx_input->GetInputValue(&vx,&gauss); - vy_input->GetInputValue(&vy,&gauss); - gr_input->GetInputValue(&groundedice,&gauss); - bs_input->GetInputValue(&bed,&gauss); sl_input->GetInputValue(&sealevel,&gauss); @@ -633,11 +572,4 @@ IssmDouble ds, db, da, dt, dw, r, R; - if(!IsIceInElement()){ - for (int iv=0;ivAddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - this->CalvingRateToVector(); - return; - } - /*Retrieve all inputs and parameters we will need*/ IssmDouble rc = FindParam(CalvingRcEnum); @@ -689,6 +621,6 @@ /*3. "Additional" crevasse opening*/ - vel = sqrt(vx*vx + vy*vy); - da = H* max(0., log(vel*365*24*3600/1600.))/log(1.2); + vel = sqrt(vx*vx + vy*vy)/(365.25*24*3600); + da = H* max(0., log(vel/1600.))/log(1.2); /*4. deal with shallow ice*/ @@ -698,5 +630,5 @@ dw = 0.; R = smb*365.25*24*3600; //convert from m/s to m/yr - if(R>1.5 && R<=3.){ + if(R>1.5 && R<=3.){ dw = 4*1.5*(R - 1.5); } @@ -707,8 +639,4 @@ /*Total calving rate*/ r = (ds+db+da+dt+dw)/H; - //if(this->Id()==1){ - // printf("rc = %g\n",rc); - // printf("ds = %g\n",ds); - //} calvingrate[iv]= mig_max * max(0., min(1., (r - rc)/(1 - rc))); //P&DC: mig_max = 3000 m/yr _assert_(!xIsNan(calvingrate[iv])); @@ -1120,84 +1048,4 @@ } /*}}}*/ -void Tria::CalvingRateCalvingMIP(){/*{{{*/ - - IssmDouble calvingrate[NUMVERTICES]; - IssmDouble calvingratex[NUMVERTICES]; - IssmDouble calvingratey[NUMVERTICES]; - int experiment = 1; /* exp:1 by default */ - int dim, domaintype; - IssmDouble vx, vy, vel, c, wrate; - IssmDouble time, groundedice, yts; - - /*Get problem dimension and whether there is moving front or not*/ - this->FindParam(&domaintype,DomainTypeEnum); - this->FindParam(&time,TimeEnum); - this->FindParam(&yts,ConstantsYtsEnum); - - switch(domaintype){ - case Domain2DverticalEnum: dim = 1; break; - case Domain2DhorizontalEnum: dim = 2; break; - case Domain3DEnum: dim = 2; break; - default: _error_("mesh "<GetInput(VxEnum); _assert_(vx_input); - Input *vy_input = this->GetInput(VyEnum); _assert_(vy_input); - Input *wrate_input = this->GetInput(CalvingAblationrateEnum); _assert_(wrate_input); - Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); - - /* Use which experiment: use existing Enum */ - this->FindParam(&experiment, CalvingUseParamEnum); - - /* Start looping on the number of vertices: */ - GaussTria gauss; - for(int iv=0;ivGetInputValue(&vx,&gauss); - vy_input->GetInputValue(&vy,&gauss); - vel=sqrt(vx*vx+vy*vy)+1.e-14; - - /* no calving for grounded ice in EXP4 */ - gr_input->GetInputValue(&groundedice,&gauss); - - switch (experiment) { - case 1: - case 3: - /* Exp 1 and 3: set c=v-wrate, wrate=0, so that w=0 */ - wrate = 0.0; - break; - case 2: - /* Exp 2: set c=v-wrate(given)*/ - wrate_input->GetInputValue(&wrate,&gauss); - break; - case 4: - /* Exp 4: set c=v-wrate(given), for the first 500 years, then c=0 for the second 500 years*/ - if((groundedice<0) && (time<=500.0*yts)) { - // wrate_input->GetInputValue(&wrate,&gauss); - wrate = -750*sin(2.0*M_PI*time/yts/1000)/yts; // m/a -> m/s - } - else { - /* no calving on the grounded ice*/ - wrate = vel; - } - break; - default: - _error_("The experiment is not supported yet!"); - } - - calvingrate[iv] = vel - wrate; - calvingratex[iv] = vx - wrate*vx/vel; - calvingratey[iv] = vy - wrate*vy/vel; - } - /*Add input*/ - this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); - this->AddInput(CalvingratexEnum,&calvingratex[0],P1DGEnum); - this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); -} -/*}}}*/ IssmDouble Tria::CharacteristicLength(void){/*{{{*/ @@ -2353,8 +2201,4 @@ //compute sea level load weights this->GetFractionGeometry(loadweights,&phi,&point1,&fraction1,&fraction2,&istrapeze1,levelset); - - //failsafe for phi so small that GetFractionGeometry returns 0 - if (phi==0) phi=1e-16; - for (int i=0;iGetBarycenterFromLevelset(platbar,plongbar, phi, fraction1, fraction2, late, longe, point1,istrapeze1,planetradius); @@ -2772,13 +2616,13 @@ IssmDouble Tria::GetIcefrontArea(){/*{{{*/ - IssmDouble bed[NUMVERTICES]; + IssmDouble bed[NUMVERTICES]; //basinId[NUMVERTICES]; IssmDouble Haverage,frontarea; IssmDouble x1,y1,x2,y2,distance; IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES]; int* indices=NULL; - - /*Return if no ice front present*/ + IssmDouble* H=NULL;; + int nrfrontbed,numiceverts; + if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0; - //if(!this->IsIcefront()) return 0.; /*Retrieve all inputs and parameters*/ @@ -2788,87 +2632,81 @@ Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum); - /*Only continue if all 3 vertices are below sea level*/ - for(int i=0;i=0.) return 0.; - - /*2. Find coordinates of where levelset crosses 0*/ - int numiceverts; - IssmDouble s[2],x[2],y[2]; - this->GetLevelsetIntersection(&indices, &numiceverts, &s[0],MaskIceLevelsetEnum,0.); - _assert_(numiceverts); - if(numiceverts>2){ - Input* ls_input = this->GetInput(MaskIceLevelsetEnum); - ls_input->Echo(); - } - - /*3 Write coordinates*/ - IssmDouble xyz_list[NUMVERTICES][3]; - ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); - int counter = 0; - if((numiceverts>0) && (numicevertsGetLevelsetIntersection(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.); + _assert_(numiceverts); + + /*3 Write coordinates*/ + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES); + int counter = 0; + if((numiceverts>0) && (numiceverts(numthk); + for(int iv=0;iv(indices); - - switch(numiceverts){ - case 1: // average over triangle - H[0]=Haux[0]; - H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]); - H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]); - Haverage=(H[1]+H[2])/2; - break; - case 2: // average over quadrangle - H[0]=Haux[0]; - H[1]=Haux[1]; - H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]); - H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]); - Haverage=(H[2]+H[3])/2; - break; - case 3: - if(counter==1) distance = 0; //front has 0 width on this element because levelset is 0 at a single vertex - else if(counter==2){ //two vertices with levelset=0: averaging ice front depth over both - Haverage = 0; - for(int i=0;i(H); _assert_(frontarea>0); @@ -4456,4 +4294,8 @@ IssmDouble Jdet; IssmDouble Jelem = 0; + IssmDouble Jelem1 = 0; + IssmDouble Jelem2 = 0; + IssmDouble Jelem3 = 0; + IssmDouble scaling = 0; IssmDouble xyz_list[NUMVERTICES][3]; GaussTria *gauss = NULL; @@ -4461,4 +4303,5 @@ /*If on water, return 0: */ if(!IsIceInElement())return 0; + /*Retrieve all inputs we will be needing: */ @@ -4481,6 +4324,98 @@ /*compute misfit between model and observation */ - Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->weight; - } + scaling=Jdet*weight*weight*gauss->weight; + Jelem+=0.5*(model-observation)*(model-observation); + //Jelem+=1e-19*log(fabs((model+1e-16)/(observation+1e-16)))*log(fabs((model+1e-16)/(observation+1e-16))); + } + Jelem*=scaling; + //Jelem1*=scaling; + //Jelem2=Jelem+Jelem1; + //_printf0_(" Tria.cpp:4327 Jtotal: "<GetInput(modelenum); _assert_(model_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + while(gauss->next()){ + /*Get parameters at gauss point*/ + model_input->GetInputValue(&model,gauss); + if (model != model) { + _printf0_(" Tria.cpp:4358 accumalator_value: " << this->accumulator_values[i] << " model: " << model << " dt: " << dt << "\n"); + } else { + //_printf0_(" Tria.cpp:4361 accumalator_value: " << this->accumulator_values[i] << " model: " << model << " dt: " << dt << "\n"); + this->accumulator_values[i++] += model * dt; + } + } + delete gauss; +} +/*}}}*/ +IssmDouble Tria::MisfitAnnual(int modelenum,int observationenum,int weightsenum, IssmDouble annual_dt){/*{{{*/ + /*Intermediaries*/ + IssmDouble model,observation,weight; + IssmDouble Jdet; + IssmDouble Jelem = 0; + IssmDouble Jelem1 = 0; + IssmDouble Jelem2 = 0; + IssmDouble Jelem3 = 0; + IssmDouble scaling = 0; + IssmDouble xyz_list[NUMVERTICES][3]; + GaussTria *gauss = NULL; + int i = 0; + + /*If on water, return 0: */ + if(!IsIceInElement())return 0; + + /*Retrieve all inputs we will be needing: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* model_input=this->GetInput(modelenum); _assert_(model_input); + Input* observation_input=this->GetInput(observationenum);_assert_(observation_input); + Input* weights_input =this->GetInput(weightsenum); _assert_(weights_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + while(gauss->next()){ + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Get parameters at gauss point*/ + observation_input->GetInputValue(&observation,gauss); + weights_input->GetInputValue(&weight,gauss); + + /*compute misfit between model and observation */ + scaling=Jdet*weight*gauss->weight; + model = this->accumulator_values[i] / annual_dt; + if (model != model ) { + _printf0_(" Tria.cpp:4405 accumalator_value: " << this->accumulator_values[i] << " model: " << model << " annual_dt: " << annual_dt << "\n"); + model = 0; + } else { + //_printf0_(" Tria.cpp:4407 accumalator_value: " << this->accumulator_values[i] << " model: " << model << " annual_dt: " << annual_dt << "\n"); + } + this->accumulator_values[i++] = 0; + Jelem+=0.5*(model-observation)*(model-observation); + //Jelem+=1e-19*log(fabs((model+1e-16)/(observation+1e-16)))*log(fabs((model+1e-16)/(observation+1e-16))); + } + Jelem*=scaling; + //Jelem1*=scaling; + //Jelem2=Jelem+Jelem1; + //_printf0_(" Tria.cpp:4327 Jtotal: "<GetInput(CalvingratexEnum); _assert_(calvingratex_input); calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); @@ -4610,10 +4543,8 @@ case DefaultCalvingEnum: case CalvingVonmisesEnum: - case CalvingVonmisesADEnum: case CalvingTestEnum: case CalvingParameterizationEnum: case CalvingLevermannEnum: case CalvingPollardEnum: - case CalvingCalvingMIPEnum: calvingratex_input->GetInputValue(&c[0],&gauss); calvingratey_input->GetInputValue(&c[1],&gauss); @@ -4729,11 +4660,6 @@ /*Do we assume that the calving front does not move if MICI is not engaged?*/ - bool regrowth = false; - bool apply_as_retreat = true; - if(!regrowth){ - movingfrontvx[iv] = 0.; - movingfrontvy[iv] = 0.; - } - + movingfrontvx[iv] = 0.; + movingfrontvy[iv] = 0.; //movingfrontvx[iv] = -2000./(365*24*3600.)*dlsf[0]/norm_dlsf; //movingfrontvy[iv] = -2000./(365*24*3600.)*dlsf[1]/norm_dlsf; @@ -4747,8 +4673,4 @@ } else if (MICI==2 && Hc>135. && bed<0. && fabs(ls)<100.e3){ // Crawford et all - - /*if 1: RETREAT rate - *if 0: calving rate*/ - if(0) v[0]=0.; v[1]=0.; /*5C Bn (worst case scenario)*/ @@ -4760,10 +4682,4 @@ movingfrontvx[iv] = v[0] -C*dlsf[0]/norm_dlsf; movingfrontvy[iv] = v[1] -C*dlsf[1]/norm_dlsf; - - /*disable regrowth if calving rate is too low*/ - if(!regrowth && C(NUMVERTICES); - IssmDouble* smbmelt_vertices = xNew(NUMVERTICES); - IssmDouble f1,f2,phi; - - Element::GetInputListOnVertices(&smbmelt_vertices[0],SmbMeltEnum); - GetFractionGeometry(weights,&phi,&point,&f1,&f2,&mainlyice,lsf); - smbmelt = 0.0; - for(int i=0;i(NUMVERTICES); - Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum); - scalefactor = 0.0; - for(int i=0;i(scalefactor_vertices); - } - else scalefactor = 1.0; - - /*Cleanup*/ - xDelete(weights); - xDelete(smbmelt_vertices); - } - else{ - /*Fully ice-covered element*/ - Input* smbmelt_input = this->GetInput(SmbMeltEnum); _assert_(smbmelt_input); - smbmelt_input->GetInputAverage(&smbmelt); // average smbmelt on element in m ice s-1 - - if(scaled==true){ - Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input); - scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element - } - else scalefactor=1.0; - } - - Total_Melt=rho_ice*base*smbmelt*scalefactor; // smbmelt on element in kg s-1 - - /*Return: */ - return Total_Melt; -} -/*}}}*/ -IssmDouble Tria::TotalSmbRefreeze(bool scaled){/*{{{*/ - - /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/ - IssmDouble base,smbrefreeze,rho_ice,scalefactor; - IssmDouble Total_Refreeze=0; - IssmDouble lsf[NUMVERTICES]; - IssmDouble xyz_list[NUMVERTICES][3]; - - /*Get material parameters :*/ - rho_ice=FindParam(MaterialsRhoIceEnum); - - if(!IsIceInElement())return 0; - - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*First calculate the area of the base (cross section triangle) - * http://en.wikipedia.org/wiki/Triangle - * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/ - base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2 - - /*Now get the average SMB over the element*/ - Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum); - if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){ - /*Partially ice-covered element*/ - bool mainlyice; - int point; - IssmDouble* weights = xNew(NUMVERTICES); - IssmDouble* smbrefreeze_vertices = xNew(NUMVERTICES); - IssmDouble f1,f2,phi; - - Element::GetInputListOnVertices(&smbrefreeze_vertices[0],SmbRefreezeEnum); - GetFractionGeometry(weights,&phi,&point,&f1,&f2,&mainlyice,lsf); - smbrefreeze = 0.0; - for(int i=0;i(NUMVERTICES); - Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum); - scalefactor = 0.0; - for(int i=0;i(scalefactor_vertices); - } - else scalefactor = 1.0; - - /*Cleanup*/ - xDelete(weights); - xDelete(smbrefreeze_vertices); - } - else{ - /*Fully ice-covered element*/ - Input* smbrefreeze_input = this->GetInput(SmbRefreezeEnum); _assert_(smbrefreeze_input); - smbrefreeze_input->GetInputAverage(&smbrefreeze); // average smbrefreeze on element in m ice s-1 - - if(scaled==true){ - Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input); - scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element - } - else scalefactor=1.0; - } - - Total_Refreeze=rho_ice*base*smbrefreeze*scalefactor; // smbrefreeze on element in kg s-1 - - /*Return: */ - return Total_Refreeze; } /*}}}*/ @@ -6718,5 +6505,5 @@ } /*}}}*/ -void Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* n_activevertices){ /*{{{*/ +void Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){ /*{{{*/ /*Declarations:{{{*/ @@ -6868,15 +6655,11 @@ } - AlphaIndex=xNewZeroInit(n_activevertices[this->lid]*nel); - if (horiz) AzimuthIndex=xNewZeroInit(n_activevertices[this->lid]*nel); + AlphaIndex=xNewZeroInit(3*nel); + if (horiz) AzimuthIndex=xNewZeroInit(3*nel); int intmax=pow(2,16)-1; - int* activevertices=xNew(n_activevertices[this->lid]); - - int av=0; for (int i=0;i<3;i++){ if(lids[this->vertices[i]->lid]==this->lid){ - activevertices[av]=i; for(int e=0;e(intmax*(atan2(dy,dx)/2/M_PI)); + AzimuthIndex[i*nel+e]=reCast(intmax*(atan2(dy,dx)/2/M_PI)); } - AlphaIndex[av*nel+e]=index; - } - av+=1; + AlphaIndex[i*nel+e]=index; + } } //for (int i=0;i<3;i++) } //for(int e=0;einputs->SetIntArrayInput(SealevelchangeConvolutionVerticesEnum,this->lid,activevertices,n_activevertices[this->lid]); - this->inputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*n_activevertices[this->lid]); - if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*n_activevertices[this->lid]); + this->inputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*3); + if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*3); /*}}}*/ @@ -7032,5 +6813,4 @@ /*Free allocations:{{{*/ #ifdef _HAVE_RESTRICT_ - delete activevertices; delete AlphaIndex; if(horiz) AzimuthIndex; @@ -7046,5 +6826,4 @@ #else - xDelete(activevertices); xDelete(AlphaIndex); if(horiz){ @@ -7078,6 +6857,4 @@ IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim; IssmDouble xyz_list[NUMVERTICES][3]; - int* activevertices = NULL; - int n_activevertices, av; #ifdef _HAVE_RESTRICT_ @@ -7154,53 +6931,48 @@ if(horiz) AzimIndex=xNew(SLGEOM_NUMLOADS); - this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices); - // 0<=n_activevertices<=3 is the number of vertices this element is in charge of computing fields in during the sea level convolutions - // activevertices contains the vertex indices (1,2 and/or 3) in case debugging is required, they are supposed to appear in the same order as slgeom->lids //Allocate: for(int l=0;lnbar[l]; - AlphaIndex[l]=xNewZeroInit(n_activevertices*nbar); - if(horiz) AzimIndex[l]=xNewZeroInit(n_activevertices*nbar); - - //av=0; - //for (int i=0;i<3;i++){ - for (int av=0;avlids[this->vertices[i]->lid]==this->lid){ - int i=activevertices[av]; - for(int e=0;elatbarycentre[l][e]; - IssmDouble longe= slgeom->longbarycentre[l][e]; - late=late/180*M_PI; - longe=longe/180*M_PI; - lati=latitude[i]; - longi=longitude[i]; + AlphaIndex[l]=xNewZeroInit(3*nbar); + if(horiz) AzimIndex[l]=xNewZeroInit(3*nbar); + + + for (int i=0;i<3;i++){ + if(slgeom->lids[this->vertices[i]->lid]==this->lid){ + for(int e=0;elatbarycentre[l][e]; + IssmDouble longe= slgeom->longbarycentre[l][e]; + late=late/180*M_PI; + longe=longe/180*M_PI; + + lati=latitude[i]; + longi=longitude[i]; + if(horiz){ - /*Compute azimuths*/ + /*Compute azimuths*/ dx=cos(lati)*sin(late)-sin(lati)*cos(late)*cos(longe-longi); dy=sin(longe-longi)*cos(late); //angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax] - AzimIndex[l][av*nbar+e]=reCast(intmax*(atan2(dy,dx)/2/M_PI)); + AzimIndex[l][i*nbar+e]=reCast(intmax*(atan2(dy,dx)/2/M_PI)); } - /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */ - delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda; - alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0))); - doubleindex=alpha/M_PI*reCast(M-1); //maps 0(doubleindex); //truncates doubleindex to integer part - - if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour - if (index==M-1){ //avoids out of bound case - index-=1; - lincoef=1; + /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */ + delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda; + alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0))); + doubleindex=alpha/M_PI*reCast(M-1); //maps 0(doubleindex); //truncates doubleindex to integer part + + if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour + if (index==M-1){ //avoids out of bound case + index-=1; + lincoef=1; + } + AlphaIndex[l][i*nbar+e]=index; } - AlphaIndex[l][av*nbar+e]=index; - //} - //av+=1; - } - + } } } @@ -7208,6 +6980,6 @@ /*Save all these arrayins for each element:*/ for (int l=0;linputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*n_activevertices); - if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*n_activevertices); + this->inputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*3); + if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*3); } /*}}}*/ @@ -7513,5 +7285,5 @@ barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp); - /*Free resources*/ + /*Free ressources*/ xDelete(areae); @@ -7555,5 +7327,4 @@ for(int i=0;iLoadWeigths[SLGEOM_OCEAN][i][this->lid]=loadweightsocean[i]; - #ifdef _ISSM_DEBUG_ /*{{{*/ /*Inform mask: */ @@ -7666,19 +7437,16 @@ oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]; } - + #ifdef _ISSM_DEBUG_ + this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum); + #endif oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid]; - oceanaverage*=rho_water*oceanarea; /*add ocean average in the global sealevelloads vector:*/ if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){ int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid]; - loads->vsubsealevelloads->SetValue(intj,oceanaverage,INS_VAL); + loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water*oceanarea,INS_VAL); loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL); } - else loads->vsealevelloads->SetValue(this->sid,oceanaverage,INS_VAL); - - #ifdef _ISSM_DEBUG_ - this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum); - #endif + else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water*oceanarea,INS_VAL); /*add ocean area into a global oceanareas vector:*/ @@ -7699,7 +7467,7 @@ IssmDouble* G=NULL; IssmDouble* Grot=NULL; - IssmDouble* rslfield=NULL; DoubleVecParam* parameter; bool computefuture=false; + int spatial_component=0; bool sal = false; @@ -7720,9 +7488,9 @@ parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL); + this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size); + for (int l=0;linputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size); if (rotation) this->inputs->GetArrayPtr(SealevelchangeGrotEnum,this->lid,&Grot,&size); - rslfield=this->SealevelchangeGxL(G,Grot,loads,polarmotionvector,slgeom,nel,computefuture=false); - this->SealevelchangeCollectGrdfield(sealevelpercpu,rslfield,slgeom,nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false); - + this->SealevelchangeGxL(sealevelpercpu, spatial_component=0,AlphaIndex, AlphaIndexsub, NULL, NULL, G, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false); } @@ -7738,4 +7506,8 @@ int nel,nbar; bool sal = false; + int* AlphaIndex=NULL; + int* AzimIndex=NULL; + int* AlphaIndexsub[SLGEOM_NUMLOADS]; + int* AzimIndexsub[SLGEOM_NUMLOADS]; int spatial_component=0; IssmDouble* G=NULL; @@ -7746,5 +7518,4 @@ IssmDouble* GNrot=NULL; IssmDouble* GErot=NULL; - IssmDouble* grdfield=NULL; DoubleVecParam* parameter; @@ -7767,4 +7538,7 @@ if(sal){ + this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size); + for (int l=0;linputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size); + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL); @@ -7775,4 +7549,7 @@ if(horiz){ + this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size); + for (int l=0;linputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size); + parameter = static_cast(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&GH,NULL); @@ -7787,21 +7564,12 @@ } } - //Relative sea level convolution - grdfield=this->SealevelchangeGxL(G,Grot,loads,polarmotionvector,slgeom,nel,computefuture=true); - this->SealevelchangeCollectGrdfield(&RSLGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true); + + this->SealevelchangeGxL(&RSLGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL,G, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true); if(elastic){ - //Bedrock Uplift - grdfield=this->SealevelchangeGxL(GU,GUrot,loads,polarmotionvector,slgeom,nel,computefuture=true); - this->SealevelchangeCollectGrdfield(&UGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true); - + this->SealevelchangeGxL(&UGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL, GU, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true); if(horiz){ - //Bedrock North displacement - grdfield=this->SealevelchangeHorizGxL(spatial_component=1,GH,GNrot,loads,polarmotionvector,slgeom,nel,computefuture=true); - this->SealevelchangeCollectGrdfield(&NGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true); - - //Bedrock East displacement - grdfield=this->SealevelchangeHorizGxL(spatial_component=2,GH,GErot,loads,polarmotionvector,slgeom,nel,computefuture=true); - this->SealevelchangeCollectGrdfield(&EGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true); + this->SealevelchangeGxL(&NGrd[0],spatial_component=1,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true); + this->SealevelchangeGxL(&EGrd[0],spatial_component=2,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true); } } @@ -7828,15 +7596,18 @@ } /*}}}*/ -IssmDouble* Tria::SealevelchangeGxL(IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture) { /*{{{*/ +void Tria::SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/ //This function performs the actual convolution between Green functions and surface Loads for a particular grd field - int* AlphaIndex=NULL; - int* AlphaIndexsub[SLGEOM_NUMLOADS]; - int* activevertices=NULL; + IssmDouble* grdfield=NULL; - int i,e,l,t,a, index, nbar, size, av,ae,b,c; + int i,e,l,t,a, index, nbar; bool rotation=false; + IssmDouble* Centroid_loads=NULL; + IssmDouble* Centroid_loads_copy=NULL; + IssmDouble* Subelement_loads[SLGEOM_NUMLOADS]; + IssmDouble* Subelement_loads_copy[SLGEOM_NUMLOADS]; + IssmDouble* horiz_projection=NULL; + IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS]; int nt=1; //important, ensures there is a defined value if computeviscous is false - int n_activevertices=0; //viscous @@ -7844,239 +7615,4 @@ int viscousindex=0; //important int viscousnumsteps=1; //important - - this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum); - this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); - - //Get green functions indexing & geometry - this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices); //the order in which the vertices appear here should be the same as in slgeom->lids - this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size); - for (int l=0;linputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size); - - if(computeviscous){ - this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum); - this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum); - if(computefuture) { - nt=viscousnumsteps-viscousindex+2; //number of time steps remaining to reach final_time, +1 is sufficient with no adaptative time stepping, +2 necessary otherwise; we assume the safe choice here for the sake of simplicity - if (nt>viscousnumsteps) nt=viscousnumsteps; - } - else nt=1; - } - //allocate - grdfield=xNewZeroInit(3*nt); - - //early return - if (n_activevertices==0) return grdfield; - - if(rotation){ //add rotational feedback - for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; - b=i*nt; - for (int m=0;m<3;m++){ //polar motion components - for(t=0;tlids[this->vertices[i]->lid]!=this->lid)continue; - b=i*nt; - c=av*nel; - for(ae=0;aenactiveloads;ae++){ - e=loads->combined_loads_index[ae]; - a=AlphaIndex[c+e]*viscousnumsteps; - for(t=0;tcombined_loads[ae]; - } - } - for(l=0;lnbar[l]; - c=av*nbar; - for (ae=0;aenactivesubloads[l];ae++){ - e=loads->combined_subloads_index[l][ae]; - a=AlphaIndexsub[l][c+e]*viscousnumsteps; - for(t=0;tcombined_subloads[l][ae]; - } - } - } - //av+=1; - } /*}}}*/ - - return grdfield; - -} /*}}}*/ -IssmDouble* Tria::SealevelchangeHorizGxL(int spatial_component, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture) { /*{{{*/ - - //This function performs the actual convolution between Green functions and surface Loads for a particular grd field - int* AlphaIndex=NULL; - int* AzimIndex=NULL; - int* AlphaIndexsub[SLGEOM_NUMLOADS]; - int* AzimIndexsub[SLGEOM_NUMLOADS]; - int* activevertices = NULL; - IssmDouble* grdfield=NULL; - int i,e,l,t,a,b,c, index, nbar, av, ae,n_activevertices, size; - bool rotation=false; - IssmDouble* projected_loads=NULL; - IssmDouble* projected_subloads[SLGEOM_NUMLOADS]; - IssmDouble* horiz_projection=NULL; - IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS]; - int nt=1; //important, ensures there is a defined value if computeviscous is false - - //viscous - bool computeviscous=false; - int viscousindex=0; //important - int viscousnumsteps=1; //important - - //Get green functions indexing & geometry - this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices); - this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size); - for (int l=0;linputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size); - this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size); - for (int l=0;linputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size); - - //First, figure out how many time steps to compute grdfield for - this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum); - this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); - if(computeviscous){ - this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum); - this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum); - if(computefuture) { - nt=viscousnumsteps-viscousindex+2; //number of time steps remaining to reach final_time, +1 is sufficient with no adaptative time stepping, +2 necessary otherwise; we assume the safe choice here for the sake of simplicity - if (nt>viscousnumsteps) nt=viscousnumsteps; - } - else nt=1; - } - //allocate - grdfield=xNewZeroInit(3*nt); - if (n_activevertices==0) return grdfield; - - if(rotation){ //add rotational feedback - for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; - for (int m=0;m<3;m++){ //polar motion components - for(t=0;t(loads->nactiveloads); - projected_loads=xNewZeroInit(loads->nactiveloads); - for(l=0;lnbar[l]; - projected_subloads[l]=xNewZeroInit(loads->nactivesubloads[l]); - horiz_projectionsub[l]=xNewZeroInit(loads->nactivesubloads[l]); - } - - - //Convolution - //av=0; - for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; - b=i*nt; - - //GxL needs to be projected on the right axis before summation into the grd field - //here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency - - //get projection - if (spatial_component==1){ //north - for(ae=0;aenactiveloads;ae++){ - e=loads->combined_loads_index[ae]; - horiz_projection[ae]=cos(2.0*M_PI*reCast(AzimIndex[av*nel+e])/65535.0); // 65535=2^16-1 = max value of 16 bits unsigned int - } - for(l=0;lnbar[l]; - for(ae=0;aenactivesubloads[l];ae++){ - e=loads->combined_subloads_index[l][ae]; - horiz_projectionsub[l][ae]=cos(2.0*M_PI*reCast(AzimIndexsub[l][av*nbar+e])/65535.0); - } - } - } - else if (spatial_component==2){ //east - for(ae=0;aenactiveloads;ae++){ - e=loads->combined_loads_index[ae]; - horiz_projection[ae]=sin(2.0*M_PI*reCast(AzimIndex[av*nel+e])/65535.0); - } - for(l=0;lnbar[l]; - for(ae=0;aenactivesubloads[l];ae++){ - e=loads->combined_subloads_index[l][ae]; - horiz_projectionsub[l][ae]=sin(2.0*M_PI*reCast(AzimIndexsub[l][av*nbar+e])/65535.0); - } - } - } - - //project load in the right direction - for (ae=0;aenactiveloads;ae++){ - projected_loads[ae]=loads->combined_loads[ae]*horiz_projection[ae]; - } - for(l=0;lnbar[l]; - for(ae=0;aenactivesubloads[l];ae++){ - projected_subloads[l][ae]=loads->combined_subloads[l][ae]*horiz_projectionsub[l][ae]; - } - } - - //do the convolution - c=av*nel; - for(ae=0;aenactiveloads;ae++){ - e=loads->combined_loads_index[ae]; - a=AlphaIndex[c+e]*viscousnumsteps; - for(t=0;tnbar[l]; - c=av*nbar; - for(ae=0;aenactivesubloads[l];ae++){ - e=loads->combined_subloads_index[l][ae]; - a=AlphaIndexsub[l][c+e]*viscousnumsteps; - for(t=0;t(horiz_projection); - xDelete(projected_loads); - for(l=0;l(projected_subloads[l]); - xDelete(horiz_projectionsub[l]); - } - return grdfield; - -} /*}}}*/ - - -void Tria::SealevelchangeCollectGrdfield(IssmDouble* grdfieldout, IssmDouble* grdfield, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/ - - //This function aligns grdfield with the requested output format: in a size 3 vector or in a size numberofvertices vector - // if compute viscous is on, we also interpolate the field timewise given the current timestepping as well as collect viscous deformation and update the viscous deformation time series for future time steps - int i,e,l,t,a, index, nbar, av, n_activevertices; - int nt=1; - - - //viscous - bool computeviscous=false; - int viscousindex=0; //important - int viscousnumsteps=1; //important - int* activevertices = NULL; IssmDouble* viscousfield=NULL; IssmDouble* grdfieldinterp=NULL; @@ -8086,8 +7622,6 @@ IssmDouble timeacc; - //parameters & initialization this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum); - this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices); - + this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); if(computeviscous){ this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum); @@ -8104,88 +7638,200 @@ else nt=1; } - - - if(!computeviscous){ /*{{{*/ - /*elastic or self attraction only case - store values computed on vertices, but don't repeat the computation if another element already computed it!:*/ - if(percpu){ - for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; - grdfieldout[this->vertices[i]->lid]=grdfield[i]; - //} - } - } - else{ - for(i=0;i(grdfield); - return; - } - else { //viscous case + //allocate + grdfield=xNewZeroInit(3*nt); + + if(rotation){ //add rotational feedback + for(int i=0;ilids[this->vertices[i]->lid]==this->lid){ + for (int m=0;m<3;m++){ //polar motion components + for(t=0;t(nel); + for (e=0;eloads[e]; + } + for(l=0;lnbar[l]; + Subelement_loads[l]=xNewZeroInit(nbar); + for (e=0;esubloads[l][e]); + } + } + if(loads->sealevelloads){ + for (e=0;esealevelloads[e]); + } + nbar=slgeom->nbar[SLGEOM_OCEAN]; + for (e=0;esubsealevelloads[e]); + } + } + + //Copy loads if dealing with a horizontal component: the result will need to be projected against the North or East axis for each vertex + if (spatial_component!=0){ + horiz_projection=xNewZeroInit(3*nel); + Centroid_loads_copy=xNewZeroInit(nel); + for (e=0;enbar[l]; + Subelement_loads_copy[l]=xNewZeroInit(nbar); + horiz_projectionsub[l]=xNewZeroInit(3*nbar); + for (e=0;elids[this->vertices[i]->lid]!=this->lid)continue; + + if (spatial_component!=0){ //horizontals /*{{{*/ + //GxL needs to be projected on the right axis before summation into the grd field + //here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency + if (spatial_component==1){ //north + for (e=0;e(AzimIndex[i*nel+e])/65535.0); // 65535=2^16-1 = max value of 16 bits unsigned int + } + for(l=0;lnbar[l]; + for (e=0;e(AzimIndexsub[l][i*nbar+e])/65535.0);; + } + } + } + else if (spatial_component==2){ //east + for (e=0;e(AzimIndex[i*nel+e])/65535.0); + } + for(l=0;lnbar[l]; + for (e=0;e(AzimIndexsub[l][i*nbar+e])/65535.0);; + } + } + } + for (e=0;enbar[l]; + for (e=0;enbar[l]; + for (e=0;elids[this->vertices[i]->lid]!=this->lid)continue; - grdfield[i*nt+0]+=viscousfield[i*viscousnumsteps+viscousindex]; - } - /* Map new grdfield generated by present-day loads onto viscous time vector*/ if(computefuture){ //viscousindex time and first time step of grdfield coincide, so just copy that value - for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; - grdfieldinterp[i*viscousnumsteps+viscousindex]=grdfield[i*nt+0]; + for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; + grdfieldinterp[i*viscousnumsteps+viscousindex]= grdfield[i*nt+0]; } if(viscoustimes[viscousindex]lids[this->vertices[i]->lid]!=this->lid)continue; - int i_time1= i*nt-viscousindex; - int i_time2= i*viscousnumsteps; - for(int t=viscousindex+1;tlids[this->vertices[i]->lid]!=this->lid)continue; + grdfieldinterp[i*viscousnumsteps+t] = (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)] + +lincoeff*grdfield[i*nt+(t-viscousindex)]; } } } + } + + /*update grdfield at present time using viscous stack at present time: */ + for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; + grdfield[i*nt+0]+=viscousfield[i*viscousnumsteps+viscousindex]; + } + + /*update viscous stack with future deformation from present load: */ + if(computefuture){ + for(int t=viscousnumsteps-1;t>=viscousindex;t--){ //we need to go backwards so as not to zero out viscousfield[i*viscousnumsteps+viscousindex] until the end + for(int i=0;ilids[this->vertices[i]->lid]!=this->lid)continue; + //offset viscousfield to remove all deformation that has already been added + viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t] + -grdfieldinterp[i*viscousnumsteps+viscousindex] + -viscousfield[i*viscousnumsteps+viscousindex]; + } + } /*Save viscous stack now that we updated the values:*/ this->inputs->SetArrayInput(viscousenum,this->lid,viscousfield,3*viscousnumsteps); - } - - /*store values computed on vertices*/ - if(percpu){ - for(av=0;avlids[this->vertices[i]->lid]!=this->lid)continue; + + /*Free resources:*/ + xDelete(grdfieldinterp); + } + } + /*}}}*/ + + /*store values computed on vertices, but don't repeat the computation if another element already computed it!:*/ + if(percpu){ + for(i=0;ilids[this->vertices[i]->lid]==this->lid){ grdfieldout[this->vertices[i]->lid]=grdfield[i*nt+0]; - //} - } - } - else{ - for(i=0;i(grdfield); + } + } + } + else{ + for(i=0;i(grdfield); + xDelete(Centroid_loads); + for(l=0;l(Subelement_loads[l]); + if (spatial_component!=0){ + xDelete(horiz_projection); + xDelete(Centroid_loads_copy); + for(l=0;l(Subelement_loads_copy[l]); + xDelete(horiz_projectionsub[l]); + } + } + if (computeviscous){ xDelete(viscoustimes); if (computefuture){ xDelete(grdfieldinterp); } - /*}}}*/ - } + } + } /*}}}*/ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Tria.h =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Tria.h (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Tria.h (revision 27956) @@ -55,5 +55,4 @@ void AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part); void CalvingRateVonmises(); - void CalvingRateVonmisesAD(); void CalvingRateTest(); void CalvingCrevasseDepth(); @@ -63,5 +62,4 @@ void CalvingMeltingFluxLevelset(); void CalvingRateParameterization(); - void CalvingRateCalvingMIP(); IssmDouble CharacteristicLength(void); void ComputeBasalStress(void); @@ -128,4 +126,6 @@ void MaterialUpdateFromTemperature(void){_error_("not implemented yet");}; IssmDouble Misfit(int modelenum,int observationenum,int weightsenum); + void MisfitAccumulate(int modelenum,IssmDouble dt); + IssmDouble MisfitAnnual(int modelenum,int observationenum,int weightsenum,IssmDouble annual_dt); IssmDouble MisfitArea(int weightsenum); int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum); @@ -158,6 +158,4 @@ IssmDouble TotalGroundedBmb(bool scaled); IssmDouble TotalSmb(bool scaled); - IssmDouble TotalSmbMelt(bool scaled); - IssmDouble TotalSmbRefreeze(bool scaled); void Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement); int UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf); @@ -176,5 +174,5 @@ #ifdef _HAVE_SEALEVELCHANGE_ void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); - void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount); + void SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids); void SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom); void SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae); @@ -249,7 +247,5 @@ void UpdateConstraintsExtrudeFromBase(void); void UpdateConstraintsExtrudeFromTop(void); - IssmDouble* SealevelchangeGxL(IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture); - IssmDouble* SealevelchangeHorizGxL(int spatial_component, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture); - void SealevelchangeCollectGrdfield(IssmDouble* grdfieldout, IssmDouble* grdfield, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture); + void SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture); /*}}}*/ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ArrayInput.cpp =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ArrayInput.cpp (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ArrayInput.cpp (revision 27956) @@ -44,6 +44,8 @@ ArrayInput* output = new ArrayInput(this->numberofelements_local); + output->N = xNew(this->numberofelements_local); xMemCpy(output->N,this->N,this->numberofelements_local); + output->values = xNew(this->numberofelements_local); for(int i=0;inumberofelements_local;i++){ if(this->values[i]){ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ControlInput.cpp =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ControlInput.cpp (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ControlInput.cpp (revision 27956) @@ -227,8 +227,4 @@ } /*}}}*/ -void ControlInput::AverageAndReplace(void){/*{{{*/ - this->values->AverageAndReplace(); -} -/*}}}*/ TriaInput* ControlInput::GetTriaInput(){/*{{{*/ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ControlInput.h =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ControlInput.h (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ControlInput.h (revision 27956) @@ -48,5 +48,4 @@ TriaInput* GetTriaInput(); PentaInput* GetPentaInput(); - void AverageAndReplace(void); }; #endif /* _CONTROLINPUT_H */ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Input.h =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Input.h (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Input.h (revision 27956) @@ -48,5 +48,4 @@ virtual void Pow(IssmDouble scale_factor){_error_("Not implemented yet");}; virtual void Scale(IssmDouble scale_factor){_error_("Not implemented yet");}; - virtual void AverageAndReplace(void){_error_("Not implemented yet");}; virtual int GetResultArraySize(void){_error_("Not implemented yet");}; Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Inputs.cpp =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Inputs.cpp (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Inputs.cpp (revision 27956) @@ -302,5 +302,5 @@ } /*}}}*/ -void Inputs::Shift(int xenum, IssmDouble alpha){/*{{{*/ +void Inputs::Shift(int xenum, IssmDouble alpha){/*{{{*/ _assert_(this); @@ -314,15 +314,4 @@ /*Shift: */ this->inputs[index_x]->Shift(alpha); -} -/*}}}*/ -void Inputs::AverageAndReplace(int inputenum){/*{{{*/ - - _assert_(this); - - /*Get indices from enums*/ - int index = EnumToIndex(inputenum); - if(!this->inputs[index]) _error_("Input "<inputs[index]->AverageAndReplace(); } /*}}}*/ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Inputs.h =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Inputs.h (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Inputs.h (revision 27956) @@ -50,5 +50,4 @@ void AXPY(IssmDouble alpha, int xenum, int yenum); void Shift(int inputenum, IssmDouble alpha); - void AverageAndReplace(int inputenum); void DeepEcho(void); void DeepEcho(int enum_in); Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/IntArrayInput.cpp =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/IntArrayInput.cpp (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/IntArrayInput.cpp (revision 27956) @@ -44,6 +44,8 @@ IntArrayInput* output = new IntArrayInput(this->numberofelements_local); + output->N = xNew(this->numberofelements_local); xMemCpy(output->N,this->N,this->numberofelements_local); + output->values = xNew(this->numberofelements_local); for(int i=0;inumberofelements_local;i++){ if(this->values[i]){ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/TransientInput.cpp =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/TransientInput.cpp (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/TransientInput.cpp (revision 27956) @@ -435,5 +435,4 @@ /*First, recover current time from parameters: */ bool linear_interp,average,cycle; - int timestepping; IssmDouble dt; this->parameters->FindParam(&linear_interp,TimesteppingInterpForcingEnum); @@ -441,20 +440,16 @@ this->parameters->FindParam(&cycle,TimesteppingCycleForcingEnum); this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*transient core time step*/ - this->parameters->FindParam(×tepping,TimesteppingTypeEnum); + + /*Change input time if we cycle through the forcing*/ + IssmDouble time0 = this->timesteps[0]; + IssmDouble time1 = this->timesteps[this->numtimesteps - 1]; + + /*We need the end time to be the last timestep that would be taken*/ + /* i.e., the case where GEMB has time stamps (finer timestep) after the last timestep */ + IssmDouble nsteps = reCast(time1/dt); + if (reCast(nsteps)timesteps[0]; - IssmDouble time1 = this->timesteps[this->numtimesteps - 1]; - - if(timestepping!=AdaptiveTimesteppingEnum){ - /*We need the end time to be the last timestep that would be taken*/ - /* i.e., the case where GEMB has time stamps (finer timestep) after the last timestep */ - /* warning: this assumes dt = constant!!*/ - IssmDouble nsteps = reCast(time1/dt); - if (reCast(nsteps)M!=this->numberofelements_local) _error_("not implemented for P1"); - - /*Get local sum and local size*/ - IssmDouble sum = 0.; - int weight; - for(int i=0;iM*this->N;i++) sum += this->values[i]; - weight = this->M*this->N; - - /*Get sum across all procs*/ - IssmDouble all_sum; - int all_weight; - ISSM_MPI_Allreduce((void*)&sum,(void*)&all_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); - ISSM_MPI_Allreduce((void*)&weight,(void*)&all_weight,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm()); - - /*Divide by number of procs*/ - IssmDouble newvalue = all_sum/reCast(all_weight); - - /*Now replace existing input*/ - this->Reset(P0Enum); - for(int i=0;iM*this->N;i++) this->values[i] = newvalue; -} -/*}}}*/ /*Object functions*/ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/TriaInput.h =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/TriaInput.h (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/TriaInput.h (revision 27956) @@ -41,5 +41,4 @@ void AXPY(Input* xinput,IssmDouble scalar); void Shift(IssmDouble scalar); - void AverageAndReplace(void); void PointWiseMult(Input* xinput); void Serve(int numindices,int* indices); Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/MisfitAnnual.cpp =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/MisfitAnnual.cpp (revision 27956) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/MisfitAnnual.cpp (revision 27956) @@ -0,0 +1,329 @@ +/*!\file MisfitAnnual.cpp + * \brief: MisfitAnnual object + */ + +/*Headers:*/ +/*{{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include /* DBL_EPSILON */ +#include "./classes.h" +#include "./ExternalResults/ExternalResult.h" +#include "./ExternalResults/Results.h" +#include "../datastructures/datastructures.h" +#include "./Elements/Element.h" +#include "./Elements/Elements.h" +#include "./FemModel.h" +#include "../modules/SurfaceAreax/SurfaceAreax.h" +#include "../classes/Params/Parameters.h" +#include "../classes/gauss/Gauss.h" +/*}}}*/ + +/*MisfitAnnual constructors, destructors :*/ +MisfitAnnual::MisfitAnnual(){/*{{{*/ + + this->definitionenum = -1; + this->name = NULL; + this->model_enum = UNDEF; + this->observation_enum = UNDEF; + this->weights_enum = UNDEF; + this->timeinterpolation=NULL; + this->local=1; + this->misfit=0; + this->annual_dt=0; + this->initialized=0; + this->lock=0; + +} +/*}}}*/ +MisfitAnnual::MisfitAnnual(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, int in_local, int in_weights_enum){/*{{{*/ + + this->definitionenum=in_definitionenum; + + this->name = xNew(strlen(in_name)+1); + xMemCpy(this->name,in_name,strlen(in_name)+1); + + this->timeinterpolation = xNew(strlen(in_timeinterpolation)+1); + xMemCpy(this->timeinterpolation,in_timeinterpolation,strlen(in_timeinterpolation)+1); + + this->model_enum=in_model_enum; + this->observation_enum=in_observation_enum; + this->weights_enum=in_weights_enum; + this->local=in_local; + + this->misfit=0; + this->annual_dt=0; + this->initialized=0; + this->lock=0; +} +/*}}}*/ +MisfitAnnual::~MisfitAnnual(){/*{{{*/ + if(this->name)xDelete(this->name); + if(this->timeinterpolation)xDelete(this->timeinterpolation); + this->misfit=0; + this->annual_dt=0; + this->initialized=0; + this->lock=0; +} +/*}}}*/ +/*Object virtual function resolutoin: */ +Object* MisfitAnnual::copy() {/*{{{*/ + MisfitAnnual* mf = new MisfitAnnual(this->name,this->definitionenum, this->model_enum,this->observation_enum,this->timeinterpolation,this->local,this->weights_enum); + mf->misfit=this->misfit; + mf->annual_dt=this->annual_dt; + mf->lock=this->lock; + return (Object*) mf; +} +/*}}}*/ +void MisfitAnnual::DeepEcho(void){/*{{{*/ + this->Echo(); +} +/*}}}*/ +void MisfitAnnual::Echo(void){/*{{{*/ + _printf_(" MisfitAnnual: " << name << " " << this->definitionenum << "\n"); + _printf_(" model_enum: " << model_enum << " " << EnumToStringx(model_enum) << "\n"); + _printf_(" observation_enum: " << observation_enum << " " << EnumToStringx(observation_enum) << "\n"); + _printf_(" weights_enum: " << weights_enum << " " << EnumToStringx(weights_enum) << "\n"); + _printf_(" timeinterpolation: " << timeinterpolation << "\n"); + _printf_(" local: " << local << "\n"); +} +/*}}}*/ +int MisfitAnnual::Id(void){/*{{{*/ + return -1; +} +/*}}}*/ +void MisfitAnnual::Marshall(MarshallHandle* marshallhandle){/*{{{*/ + _error_("not implemented yet!"); +} +/*}}}*/ +int MisfitAnnual::ObjectEnum(void){/*{{{*/ + return MisfitannualEnum; +} +/*}}}*/ +/*Definition virtual function resolutoin: */ +int MisfitAnnual::DefinitionEnum(){/*{{{*/ + return this->definitionenum; +} +/*}}}*/ +char* MisfitAnnual::Name(){/*{{{*/ + char* name2=xNew(strlen(this->name)+1); + xMemCpy(name2,this->name,strlen(this->name)+1); + + return name2; +} +/*}}}*/ +IssmDouble MisfitAnnual::Response(FemModel* femmodel){/*{{{*/ + + /*diverse: */ + IssmDouble time,starttime,finaltime; + IssmDouble dt; + + /*recover time parameters: */ + femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); + femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + + if (this->local==1){ /*area integration using elements: {{{*/ + /*Do we lock? i.e. are we at final_time? :*/ + _printf0_(" MisfitAnnual::133 lock: " << this->lock << "\n"); + if(this->lock>=1) { + IssmDouble misfit_t=0.; + IssmDouble all_misfit_t=0.; + IssmDouble area_t=0.; + IssmDouble all_area_t=0.; + IssmDouble area_e = 0.; + + + /*parameters: */ + IssmDouble starttime,finaltime,dt,yts; + bool iscontrol,isautodiff; + int timestepping; + int output_frequency,checkpoint_frequency; + int amr_frequency; + char **requested_outputs = NULL; + IssmDouble* element_accum_values = xNewZeroInit(3); + + /*intermediary: */ + int step; + IssmDouble time; + IssmDouble annual_dt; + + /*first, figure out if there was a check point, if so, do a reset of the FemModel* femmodel structure. */ + // femmodel->parameters->FindParam(&checkpoint_frequency,SettingsCheckpointFrequencyEnum); + // if(checkpoint_frequency) femmodel->Restart(); + + /*then recover parameters common to all solutions*/ + // femmodel->parameters->FindParam(&step,StepEnum); + // femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); + femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); + femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum); + femmodel->parameters->FindParam(×tepping,TimesteppingTypeEnum); + + /*initialize: */ + step = 0; + time = 0; + annual_dt = 0; + misfit = 0; + time = starttime; + femmodel->parameters->SetParam(starttime,TimeEnum); + femmodel->parameters->SetParam(step,StepEnum); + while(time < finaltime + (yts*DBL_EPSILON)){ //make sure we run up to finaltime. + /*Time Increment*/ + switch(timestepping){ + case AdaptiveTimesteppingEnum: + femmodel->TimeAdaptx(&dt); + if(time+dt>finaltime) dt=finaltime-time; + femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum); + break; + case FixedTimesteppingEnum: + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + break; + default: + _error_("Time stepping \""<parameters->SetParam(time,TimeEnum); + femmodel->parameters->SetParam(step,StepEnum); + + _printf0_(" MisfitAnnual::197 accumulating...\n"); + //Accumulate mean misfit on each element + for(Object* & object : femmodel->elements->objects){ + Element* element=xDynamicCast(object); + element->MisfitAccumulate(model_enum,dt); + } + + if (annual_dt >= 0) { + //if (annual_dt >= (yts - yts*DBL_EPSILON)) { + //Return mean misfit on each element + misfit_t = 0; + area_t = 0; + + for(Object* & object : femmodel->elements->objects){ + Element* element=xDynamicCast(object); + area_e=element->MisfitArea(weights_enum); + misfit_t+=element->MisfitAnnual(model_enum,observation_enum,weights_enum,annual_dt); + area_t+=area_e; + } + ISSM_MPI_Allreduce ( (void*)&misfit_t,(void*)&all_misfit_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + ISSM_MPI_Allreduce ( (void*)&area_t,(void*)&all_area_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + area_t=all_area_t; + misfit_t=all_misfit_t; + + /*Divide by surface area if not nill!: */ + if (area_t!=0) misfit_t=misfit_t/area_t; + + /*Add this time's contribution to curent misfit: */ + misfit+=misfit_t/annual_dt/(finaltime-starttime); + _printf0_(" MisfitAnnual::216 step: " << step << " misfit: " << misfit << " misfit_contrib:" << misfit_t*annual_dt/(finaltime-starttime) << "\n"); + annual_dt = 0; + } + } + _printf0_(" MisfitAnnual::227 returning misfit: " << misfit << "\n"); + return misfit; + } + if(time==finaltime)this->lock+=1; + return misfit; + } /*}}}*/ + else if (this->local==2){ /*vertex by vertex computation: {{{*/ + + IssmDouble* model = NULL; + IssmDouble* observation= NULL; + IssmDouble* weights= NULL; + int msize,osize,wsize; + + /*Are we transient?:*/ + if (time==0){ + IssmDouble misfit_t=0.; + + /*get global vectors: */ + GetVectorFromInputsx(&model,&msize,femmodel,model_enum); + GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize); + GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize); + + int count=0; + for (int i=0;i(model); + xDelete(observation); + xDelete(weights); + + /*return value: */ + return misfit; + } + else{ + + IssmDouble misfit_t=0.; + IssmDouble all_misfit_t=0.; + + /*If we are locked, return time average: */ + if(this->lock)return misfit/(time-starttime); + + /*get global vectors: */ + GetVectorFromInputsx(&model,&msize,femmodel,model_enum); + GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize); + GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize); + + int count=0; + for (int i=0;ilock=1; + + /*Free resources:*/ + xDelete(model); + xDelete(observation); + xDelete(weights); + + /*What we return is the value of misfit / time: */ + return misfit/(time-starttime); + } + + } /*}}}*/ + else{ /*global computation: {{{ */ + + IssmDouble model, observation; + + /*If we are locked, return time average: */ + if(this->lock) return misfit/(time-starttime); + + /*First, the global model response: */ + model=OutputDefinitionsResponsex(femmodel,this->model_enum); + /*Now, the observation is buried inside the elements, go fish it in the first element (cludgy, needs fixing): */ + Element* element = (Element*)femmodel->elements->GetObjectByOffset(0); _assert_(element); + Input* input = element->GetInput(observation_enum); _assert_(input); + input->GetInputAverage(&observation); + + /*Add this time's contribution to curent misfit: */ + misfit+=dt*(model-observation); + + /*Do we lock? i.e. are we at final_time? :*/ + if(time==finaltime)this->lock=1; + + /*What we return is the value of misfit / time: */ + return misfit/(time-starttime); + } /*}}}*/ + +} + /*}}}*/ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/MisfitAnnual.h =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/MisfitAnnual.h (revision 27956) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/MisfitAnnual.h (revision 27956) @@ -0,0 +1,50 @@ +/*!\file MisfitAnnual.h + * \brief: header file for MisfitAnnual object + */ + +#ifndef _MISFITANNUAL_H_ +#define _MISFITANNUAL_H_ + +/*Headers:*/ +#include "./Definition.h" +#include "./FemModel.h" + +IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,int output_enum); +void GetVectorFromInputsx( IssmDouble** pvector, int* pvector_size, FemModel* femmodel,int name); + +class MisfitAnnual: public Object, public Definition{ + + public: + + int definitionenum; + int local; + int model_enum; + char* name; + int observation_enum; + char* timeinterpolation; + int weights_enum; + + int lock; // if lock is on, we just return the value stored in "misfit". this is used so we don't compute misfit past the final_time + IssmDouble misfit; //value carried over in time. + IssmDouble annual_dt; + IssmDouble initialized; + + /*MisfitAnnual constructors, destructors :*/ + MisfitAnnual(); + MisfitAnnual(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, int in_local, int in_weights_enum); + ~MisfitAnnual(); + + /*Object virtual function resolutoin: */ + Object* copy(); + void DeepEcho(void); + void Echo(void); + int Id(void); + void Marshall(MarshallHandle* marshallhandle); + int ObjectEnum(void); + + /*Definition virtual function resolutoin: */ + int DefinitionEnum(); + char* Name(); + IssmDouble Response(FemModel* femmodel); +}; +#endif /* _MISFIT_H_ */ Index: /issm/branches/trunk-dlcheng-ASE/src/c/classes/classes.h =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/classes/classes.h (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/classes/classes.h (revision 27956) @@ -18,4 +18,5 @@ #include "./Massfluxatgate.h" #include "./Misfit.h" +#include "./MisfitAnnual.h" #include "./SealevelGeometry.h" #include "./GrdLoads.h" @@ -24,9 +25,5 @@ #include "./Numberedcostfunction.h" #include "./Cfsurfacesquare.h" -#include "./Cfsurfacesquaretransient.h" #include "./Cfdragcoeffabsgrad.h" -#include "./Cfdragcoeffabsgradtransient.h" -#include "./Cfrheologybbarabsgrad.h" -#include "./Cfrheologybbarabsgradtransient.h" #include "./Cfsurfacelogvel.h" #include "./Cflevelsetmisfit.h" @@ -93,5 +90,4 @@ #include "./Params/GenericParam.h" #include "./Params/BoolParam.h" -#include "./Params/ControlParam.h" #include "./Params/DoubleMatParam.h" #include "./Params/DoubleTransientMatParam.h" Index: /issm/branches/trunk-dlcheng-ASE/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 27956) @@ -154,4 +154,96 @@ /*}}}*/ } + else if (output_definition_enums[i]==MisfitannualEnum){ + /*Deal with misfits: {{{*/ + + /*misfit variables: */ + int nummisfits; + char** misfit_name_s = NULL; + char** misfit_definitionstring_s = NULL; + char** misfit_model_string_s = NULL; + IssmDouble** misfit_observation_s = NULL; + char** misfit_observation_string_s = NULL; + int* misfit_observation_M_s = NULL; + int* misfit_observation_N_s = NULL; + int* misfit_local_s = NULL; + char** misfit_timeinterpolation_s = NULL; + IssmDouble** misfit_weights_s = NULL; + int* misfit_weights_M_s = NULL; + int* misfit_weights_N_s = NULL; + char** misfit_weights_string_s = NULL; + + /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/misfit.m): */ + iomodel->FetchMultipleData(&misfit_name_s,&nummisfits, "md.misfitannual.name"); + iomodel->FetchMultipleData(&misfit_definitionstring_s,&nummisfits, "md.misfitannual.definitionstring"); + iomodel->FetchMultipleData(&misfit_model_string_s,&nummisfits, "md.misfitannual.model_string"); + iomodel->FetchMultipleData(&misfit_observation_s,&misfit_observation_M_s,&misfit_observation_N_s,&nummisfits, "md.misfitannual.observation"); + iomodel->FetchMultipleData(&misfit_observation_string_s,&nummisfits, "md.misfitannual.observation_string"); + iomodel->FetchMultipleData(&misfit_timeinterpolation_s,&nummisfits, "md.misfitannual.timeinterpolation"); + iomodel->FetchMultipleData(&misfit_local_s,&nummisfits, "md.misfitannual.local"); + iomodel->FetchMultipleData(&misfit_weights_s,&misfit_weights_M_s,&misfit_weights_N_s,&nummisfits, "md.misfitannual.weights"); + iomodel->FetchMultipleData(&misfit_weights_string_s,&nummisfits, "md.misfitannual.weights_string"); + + for(j=0;jnumberofvertices) || (misfit_observation_M_s[j]==iomodel->numberofvertices+1)){ + obs_vector_type=1; + } + else if ((misfit_observation_M_s[j]==iomodel->numberofelements) || (misfit_observation_M_s[j]==iomodel->numberofelements+1)){ + obs_vector_type=2; + } + else + _error_("misfit observation size not supported yet"); + + int weight_vector_type=0; + if ((misfit_weights_M_s[j]==iomodel->numberofvertices) || (misfit_weights_M_s[j]==iomodel->numberofvertices+1)){ + weight_vector_type=1; + } + else if ((misfit_weights_M_s[j]==iomodel->numberofelements) || (misfit_weights_M_s[j]==iomodel->numberofelements+1)){ + weight_vector_type=2; + } + else + _error_("misfit weight size not supported yet"); + + /*First create a misfit object for that specific string (misfit_model_string_s[j]):*/ + output_definitions->AddObject(new MisfitAnnual(misfit_name_s[j],StringToEnumx(misfit_definitionstring_s[j]),StringToEnumx(misfit_model_string_s[j]),StringToEnumx(misfit_observation_string_s[j]),misfit_timeinterpolation_s[j],misfit_local_s[j],StringToEnumx(misfit_weights_string_s[j]))); + + /*Now, for this particular misfit object, make sure we plug into the elements: the observation, and the weights.*/ + for(Object* & object : elements->objects){ + Element* element=xDynamicCast(object); + element->InputCreate(misfit_observation_s[j],inputs,iomodel,misfit_observation_M_s[j],misfit_observation_N_s[j],obs_vector_type,StringToEnumx(misfit_observation_string_s[j]),7); + element->InputCreate(misfit_weights_s[j],inputs,iomodel,misfit_weights_M_s[j],misfit_weights_N_s[j],weight_vector_type,StringToEnumx(misfit_weights_string_s[j]),7); + } + + } + + /*Free resources:*/ + for(j=0;j(string); + string = misfit_observation_string_s[j]; xDelete(string); + string = misfit_model_string_s[j]; xDelete(string); + string = misfit_weights_string_s[j]; xDelete(string); + string = misfit_name_s[j]; xDelete(string); + string = misfit_timeinterpolation_s[j]; xDelete(string); + matrix = misfit_observation_s[j]; xDelete(matrix); + matrix = misfit_weights_s[j]; xDelete(matrix); + } + xDelete(misfit_name_s); + xDelete(misfit_model_string_s); + xDelete(misfit_definitionstring_s); + xDelete(misfit_observation_s); + xDelete(misfit_observation_string_s); + xDelete(misfit_observation_M_s); + xDelete(misfit_observation_N_s); + xDelete(misfit_local_s); + xDelete(misfit_timeinterpolation_s); + xDelete(misfit_weights_s); + xDelete(misfit_weights_M_s); + xDelete(misfit_weights_N_s); + xDelete(misfit_weights_string_s); + /*}}}*/ + } else if (output_definition_enums[i]==CfsurfacesquareEnum){ /*Deal with cfsurfacesquare: {{{*/ @@ -205,11 +297,11 @@ /*First create a cfsurfacesquare object for that specific string (cfsurfacesquare_model_string_s[j]):*/ - output_definitions->AddObject(new Cfsurfacesquare(cfsurfacesquare_name_s[j],StringToEnumx(cfsurfacesquare_definitionstring_s[j]),StringToEnumx(cfsurfacesquare_model_string_s[j]),cfsurfacesquare_datatime_s[j])); + output_definitions->AddObject(new Cfsurfacesquare(cfsurfacesquare_name_s[j],StringToEnumx(cfsurfacesquare_definitionstring_s[j]),StringToEnumx(cfsurfacesquare_model_string_s[j]),StringToEnumx(cfsurfacesquare_observation_string_s[j]),StringToEnumx(cfsurfacesquare_weights_string_s[j]),cfsurfacesquare_datatime_s[j],false)); /*Now, for this particular cfsurfacesquare object, make sure we plug into the elements: the observation, and the weights.*/ for(Object* & object : elements->objects){ Element* element=xDynamicCast(object); - element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_observation_s[j],inputs,iomodel,cfsurfacesquare_observation_M_s[j],cfsurfacesquare_observation_N_s[j],obs_vector_type,StringToEnumx(cfsurfacesquare_observation_string_s[j]),SurfaceObservationEnum); - element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_weights_s[j],inputs,iomodel,cfsurfacesquare_weights_M_s[j],cfsurfacesquare_weights_N_s[j],weight_vector_type,StringToEnumx(cfsurfacesquare_weights_string_s[j]),WeightsSurfaceObservationEnum); + element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_observation_s[j],inputs,iomodel,cfsurfacesquare_observation_M_s[j],cfsurfacesquare_observation_N_s[j],obs_vector_type,StringToEnumx(cfsurfacesquare_observation_string_s[j]),7,SurfaceObservationEnum); + element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_weights_s[j],inputs,iomodel,cfsurfacesquare_weights_M_s[j],cfsurfacesquare_weights_N_s[j],weight_vector_type,StringToEnumx(cfsurfacesquare_weights_string_s[j]),7,WeightsSurfaceObservationEnum); } @@ -244,71 +336,4 @@ /*}}}*/ } - else if (output_definition_enums[i]==CfsurfacesquaretransientEnum){ - /*Deal with cfsurfacesquaretransient: {{{*/ - - /*cfsurfacesquaretransient variables: */ - int num_cfsurfacesquaretransients,test; - char **cfssqt_name_s = NULL; - char **cfssqt_definitionstring_s = NULL; - char **cfssqt_model_string_s = NULL; - IssmDouble **cfssqt_observations_s = NULL; - int *cfssqt_observations_M_s = NULL; - int *cfssqt_observations_N_s = NULL; - IssmDouble **cfssqt_weights_s = NULL; - int *cfssqt_weights_M_s = NULL; - int *cfssqt_weights_N_s = NULL; - - /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfsurfacesquaretransient.m): */ - iomodel->FetchMultipleData(&cfssqt_name_s,&num_cfsurfacesquaretransients,"md.cfsurfacesquaretransient.name"); - iomodel->FetchMultipleData(&cfssqt_definitionstring_s,&test,"md.cfsurfacesquaretransient.definitionstring"); _assert_(test==num_cfsurfacesquaretransients); - iomodel->FetchMultipleData(&cfssqt_model_string_s,&test,"md.cfsurfacesquaretransient.model_string"); _assert_(test==num_cfsurfacesquaretransients); - iomodel->FetchMultipleData(&cfssqt_observations_s,&cfssqt_observations_M_s,&cfssqt_observations_N_s,&test, "md.cfsurfacesquaretransient.observations"); _assert_(test==num_cfsurfacesquaretransients); - iomodel->FetchMultipleData(&cfssqt_weights_s,&cfssqt_weights_M_s,&cfssqt_weights_N_s, &test,"md.cfsurfacesquaretransient.weights"); _assert_(test==num_cfsurfacesquaretransients); - - for(j=0;jnumberofvertices+1)) _error_("observations should be a P1 time series"); - if (cfssqt_weights_M_s[j]!=iomodel->numberofvertices+1) _error_("weights should be a P1 time series"); - _assert_(cfssqt_observations_N_s[j]>0); - - /*extract data times from last row of observations*/ - IssmDouble *datatimes = xNew(cfssqt_observations_N_s[j]); - for(int k=0;kAddObject(new Cfsurfacesquaretransient(cfssqt_name_s[j], StringToEnumx(cfssqt_definitionstring_s[j]), StringToEnumx(cfssqt_model_string_s[j]), cfssqt_observations_N_s[j],datatimes )); - xDelete(datatimes); - - /*Now, for this particular cfsurfacesquaretransient object, make sure we plug into the elements: the observation, and the weights.*/ - for(Object* & object : elements->objects){ - Element* element=xDynamicCast(object); - element->DatasetInputAdd(StringToEnumx(cfssqt_definitionstring_s[j]),cfssqt_observations_s[j],inputs,iomodel,cfssqt_observations_M_s[j],cfssqt_observations_N_s[j],1,SurfaceObservationEnum,SurfaceObservationEnum); - element->DatasetInputAdd(StringToEnumx(cfssqt_definitionstring_s[j]),cfssqt_weights_s[j],inputs,iomodel,cfssqt_weights_M_s[j],cfssqt_weights_N_s[j],1,WeightsSurfaceObservationEnum,WeightsSurfaceObservationEnum); - - } - } - - /*Free resources:*/ - for(j=0;j(string); - string = cfssqt_model_string_s[j]; xDelete(string); - string = cfssqt_name_s[j]; xDelete(string); - matrix = cfssqt_observations_s[j]; xDelete(matrix); - matrix = cfssqt_weights_s[j]; xDelete(matrix); - } - xDelete(cfssqt_name_s); - xDelete(cfssqt_model_string_s); - xDelete(cfssqt_definitionstring_s); - xDelete(cfssqt_observations_s); - xDelete(cfssqt_observations_M_s); - xDelete(cfssqt_observations_N_s); - xDelete(cfssqt_weights_s); - xDelete(cfssqt_weights_M_s); - xDelete(cfssqt_weights_N_s); - /*}}}*/ - } else if (output_definition_enums[i]==CfdragcoeffabsgradEnum){ /*Deal with cfdragcoeffabsgrad: {{{*/ @@ -342,5 +367,5 @@ /*First create a cfdragcoeffabsgrad object for that specific string (cfdragcoeffabsgrad_model_string_s[j]):*/ - output_definitions->AddObject(new Cfdragcoeffabsgrad(cfdragcoeffabsgrad_name_s[j],StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]))); + output_definitions->AddObject(new Cfdragcoeffabsgrad(cfdragcoeffabsgrad_name_s[j],StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),false)); /*Now, for this particular cfdragcoeffabsgrad object, make sure we plug into the elements: the observation, and the weights.*/ @@ -349,5 +374,5 @@ Element* element=xDynamicCast(object); - element->DatasetInputAdd(StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),cfdragcoeffabsgrad_weights_s[j],inputs,iomodel,cfdragcoeffabsgrad_weights_M_s[j],cfdragcoeffabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),WeightsSurfaceObservationEnum); + element->DatasetInputAdd(StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),cfdragcoeffabsgrad_weights_s[j],inputs,iomodel,cfdragcoeffabsgrad_weights_M_s[j],cfdragcoeffabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),7,WeightsSurfaceObservationEnum); } @@ -373,173 +398,4 @@ /*}}}*/ } - else if (output_definition_enums[i]==CfdragcoeffabsgradtransientEnum){ - /*Deal with cfdragcoeffabsgradtransient: {{{*/ - - /*cfdragcoeffabsgrad variables: */ - int num_cfdragcoeffabsgradtransients, test; - char** cfdraggradt_name_s = NULL; - char** cfdraggradt_definitionstring_s = NULL; - IssmDouble** cfdraggradt_weights_s = NULL; - int* cfdraggradt_weights_M_s = NULL; - int* cfdraggradt_weights_N_s = NULL; - - /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfdragcoeffabsgradtransient.m): */ - iomodel->FetchMultipleData(&cfdraggradt_name_s,&num_cfdragcoeffabsgradtransients, "md.cfdragcoeffabsgradtransient.name"); - iomodel->FetchMultipleData(&cfdraggradt_definitionstring_s,&num_cfdragcoeffabsgradtransients, "md.cfdragcoeffabsgradtransient.definitionstring"); - iomodel->FetchMultipleData(&cfdraggradt_weights_s,&cfdraggradt_weights_M_s,&cfdraggradt_weights_N_s,&test, "md.cfdragcoeffabsgradtransient.weights"); - - for(j=0;jnumberofvertices+1) _error_("weights should be a P1 time series"); - - /*extract data times from last row of observations*/ - IssmDouble *datatimes = xNew(cfdraggradt_weights_N_s[j]); - for(int k=0;kAddObject(new Cfdragcoeffabsgradtransient(cfdraggradt_name_s[j],StringToEnumx(cfdraggradt_definitionstring_s[j]), cfdraggradt_weights_N_s[j], datatimes)); - - /*Now, for this particular cfdragcoeffabsgrad object, make sure we plug into the elements: the observation, and the weights.*/ - for(Object* & object : elements->objects){ - - Element* element=xDynamicCast(object); - - element->DatasetInputAdd(StringToEnumx(cfdraggradt_definitionstring_s[j]),cfdraggradt_weights_s[j],inputs,iomodel,cfdraggradt_weights_M_s[j],cfdraggradt_weights_N_s[j],1,WeightsSurfaceObservationEnum,WeightsSurfaceObservationEnum); - - } - } - - /*Free resources:*/ - for(j=0;j(string); - string = cfdraggradt_name_s[j]; xDelete(string); - matrix = cfdraggradt_weights_s[j]; xDelete(matrix); - } - xDelete(cfdraggradt_name_s); - xDelete(cfdraggradt_definitionstring_s); - xDelete(cfdraggradt_weights_s); - xDelete(cfdraggradt_weights_M_s); - xDelete(cfdraggradt_weights_N_s); - /*}}}*/ - } - else if (output_definition_enums[i]==CfrheologybbarabsgradEnum){ - /*Deal with cfrheologybbarabsgrad: {{{*/ - - /*cfrheologybbarabsgrad variables: */ - int num_cfrheologybbarabsgrads; - char** cfrheologybbarabsgrad_name_s = NULL; - char** cfrheologybbarabsgrad_definitionstring_s = NULL; - IssmDouble** cfrheologybbarabsgrad_weights_s = NULL; - int* cfrheologybbarabsgrad_weights_M_s = NULL; - int* cfrheologybbarabsgrad_weights_N_s = NULL; - char** cfrheologybbarabsgrad_weights_string_s = NULL; - - /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfrheologybbarabsgrad.m): */ - iomodel->FetchMultipleData(&cfrheologybbarabsgrad_name_s,&num_cfrheologybbarabsgrads, "md.cfrheologybbarabsgrad.name"); - iomodel->FetchMultipleData(&cfrheologybbarabsgrad_definitionstring_s,&num_cfrheologybbarabsgrads, "md.cfrheologybbarabsgrad.definitionstring"); - iomodel->FetchMultipleData(&cfrheologybbarabsgrad_weights_s,&cfrheologybbarabsgrad_weights_M_s,&cfrheologybbarabsgrad_weights_N_s,&num_cfrheologybbarabsgrads, "md.cfrheologybbarabsgrad.weights"); - iomodel->FetchMultipleData(&cfrheologybbarabsgrad_weights_string_s,&num_cfrheologybbarabsgrads, "md.cfrheologybbarabsgrad.weights_string"); - - for(j=0;jnumberofvertices) || (cfrheologybbarabsgrad_weights_M_s[j]==iomodel->numberofvertices+1)){ - weight_vector_type=1; - } - else if ((cfrheologybbarabsgrad_weights_M_s[j]==iomodel->numberofelements) || (cfrheologybbarabsgrad_weights_M_s[j]==iomodel->numberofelements+1)){ - weight_vector_type=2; - } - else - _error_("cfrheologybbarabsgrad weight size not supported yet"); - - /*First create a cfrheologybbarabsgrad object for that specific string (cfrheologybbarabsgrad_model_string_s[j]):*/ - output_definitions->AddObject(new Cfrheologybbarabsgrad(cfrheologybbarabsgrad_name_s[j],StringToEnumx(cfrheologybbarabsgrad_definitionstring_s[j]),StringToEnumx(cfrheologybbarabsgrad_weights_string_s[j]))); - - /*Now, for this particular cfrheologybbarabsgrad object, make sure we plug into the elements: the observation, and the weights.*/ - for(Object* & object : elements->objects){ - - Element* element=xDynamicCast(object); - - element->DatasetInputAdd(StringToEnumx(cfrheologybbarabsgrad_definitionstring_s[j]),cfrheologybbarabsgrad_weights_s[j],inputs,iomodel,cfrheologybbarabsgrad_weights_M_s[j],cfrheologybbarabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfrheologybbarabsgrad_weights_string_s[j]),WeightsSurfaceObservationEnum); - - } - - } - /*Free resources:*/ - for(j=0;j(string); - string = cfrheologybbarabsgrad_weights_string_s[j]; xDelete(string); - string = cfrheologybbarabsgrad_name_s[j]; xDelete(string); - matrix = cfrheologybbarabsgrad_weights_s[j]; xDelete(matrix); - } - xDelete(cfrheologybbarabsgrad_name_s); - xDelete(cfrheologybbarabsgrad_definitionstring_s); - xDelete(cfrheologybbarabsgrad_weights_s); - xDelete(cfrheologybbarabsgrad_weights_M_s); - xDelete(cfrheologybbarabsgrad_weights_N_s); - xDelete(cfrheologybbarabsgrad_weights_string_s); - /*}}}*/ - } - else if (output_definition_enums[i]==CfrheologybbarabsgradtransientEnum){ - /*Deal with cfrheologybbarabsgradtransient: {{{*/ - - /*cfrheologybbarabsgrad variables: */ - int num_cfrheologybbarabsgradtransients, test; - char** cfrheogradt_name_s = NULL; - char** cfrheogradt_definitionstring_s = NULL; - IssmDouble** cfrheogradt_weights_s = NULL; - int* cfrheogradt_weights_M_s = NULL; - int* cfrheogradt_weights_N_s = NULL; - char** cfrheogradt_weights_string_s = NULL; - - /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfrheologybbarabsgradtransient.m): */ - iomodel->FetchMultipleData(&cfrheogradt_name_s,&num_cfrheologybbarabsgradtransients, "md.cfrheologybbarabsgradtransient.name"); - iomodel->FetchMultipleData(&cfrheogradt_definitionstring_s,&num_cfrheologybbarabsgradtransients, "md.cfrheologybbarabsgradtransient.definitionstring"); - iomodel->FetchMultipleData(&cfrheogradt_weights_s,&cfrheogradt_weights_M_s,&cfrheogradt_weights_N_s,&test, "md.cfrheologybbarabsgradtransient.weights"); - - for(j=0;jnumberofvertices+1) _error_("weights should be a P1 time series"); - - /*extract data times from last row of observations*/ - IssmDouble *datatimes = xNew(cfrheogradt_weights_N_s[j]); - for(int k=0;kAddObject(new Cfrheologybbarabsgradtransient(cfrheogradt_name_s[j],StringToEnumx(cfrheogradt_definitionstring_s[j]), cfrheogradt_weights_N_s[j], datatimes)); - - /*Now, for this particular cfrheologybbarabsgrad object, make sure we plug into the elements: the observation, and the weights.*/ - for(Object* & object : elements->objects){ - - Element* element=xDynamicCast(object); - - element->DatasetInputAdd(StringToEnumx(cfrheogradt_definitionstring_s[j]),cfrheogradt_weights_s[j],inputs,iomodel,cfrheogradt_weights_M_s[j],cfrheogradt_weights_N_s[j],1,WeightsSurfaceObservationEnum,WeightsSurfaceObservationEnum); - - } - } - - /*Free resources:*/ - for(j=0;j(string); - string = cfrheogradt_name_s[j]; xDelete(string); - matrix = cfrheogradt_weights_s[j]; xDelete(matrix); - } - xDelete(cfrheogradt_name_s); - xDelete(cfrheogradt_definitionstring_s); - xDelete(cfrheogradt_weights_s); - xDelete(cfrheogradt_weights_M_s); - xDelete(cfrheogradt_weights_N_s); - /*}}}*/ - } else if (output_definition_enums[i]==CfsurfacelogvelEnum){ /*Deal with cfsurfacelogvel: {{{*/ @@ -595,5 +451,5 @@ /*First create a cfsurfacelogvel object for that specific string (cfsurfacelogvel_modeltring[j]):*/ - output_definitions->AddObject(new Cfsurfacelogvel(cfsurfacelogvel_name[j],StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_datatime[j])); + output_definitions->AddObject(new Cfsurfacelogvel(cfsurfacelogvel_name[j],StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_datatime[j],false)); /*Now, for this particular cfsurfacelogvel object, make sure we plug into the elements: the observation, and the weights.*/ @@ -602,7 +458,7 @@ Element* element=xDynamicCast(object); - element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vxobs[j],inputs,iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vxobs_string[j]),VxObsEnum); - element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vyobs[j],inputs,iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vyobs_string[j]),VyObsEnum); - element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_weights[j],inputs,iomodel,cfsurfacelogvel_weights_M[j],cfsurfacelogvel_weights_N[j],weight_vector_type,StringToEnumx(cfsurfacelogvel_weightstring[j]),WeightsSurfaceObservationEnum); + element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vxobs[j],inputs,iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vxobs_string[j]),7,VxObsEnum); + element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vyobs[j],inputs,iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vyobs_string[j]),7,VyObsEnum); + element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_weights[j],inputs,iomodel,cfsurfacelogvel_weights_M[j],cfsurfacelogvel_weights_N[j],weight_vector_type,StringToEnumx(cfsurfacelogvel_weightstring[j]),7,WeightsSurfaceObservationEnum); } @@ -689,11 +545,11 @@ /*First create a cflevelsetmisfit object for that specific string (cflevelsetmisfit_model_string_s[j]):*/ - output_definitions->AddObject(new Cflevelsetmisfit(cflevelsetmisfit_name_s[j],StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),StringToEnumx(cflevelsetmisfit_model_string_s[j]),cflevelsetmisfit_datatime_s[j])); + output_definitions->AddObject(new Cflevelsetmisfit(cflevelsetmisfit_name_s[j],StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),StringToEnumx(cflevelsetmisfit_model_string_s[j]),StringToEnumx(cflevelsetmisfit_observation_string_s[j]),StringToEnumx(cflevelsetmisfit_weights_string_s[j]),cflevelsetmisfit_datatime_s[j],false)); /*Now, for this particular cflevelsetmisfit object, make sure we plug into the elements: the observation, and the weights.*/ for(Object* & object : elements->objects){ Element* element=xDynamicCast(object); - element->DatasetInputAdd(StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),cflevelsetmisfit_observation_s[j],inputs,iomodel,cflevelsetmisfit_observation_M_s[j],cflevelsetmisfit_observation_N_s[j],obs_vector_type,StringToEnumx(cflevelsetmisfit_observation_string_s[j]),LevelsetObservationEnum); - element->DatasetInputAdd(StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),cflevelsetmisfit_weights_s[j],inputs,iomodel,cflevelsetmisfit_weights_M_s[j],cflevelsetmisfit_weights_N_s[j],weight_vector_type,StringToEnumx(cflevelsetmisfit_weights_string_s[j]),WeightsLevelsetObservationEnum); + element->DatasetInputAdd(StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),cflevelsetmisfit_observation_s[j],inputs,iomodel,cflevelsetmisfit_observation_M_s[j],cflevelsetmisfit_observation_N_s[j],obs_vector_type,StringToEnumx(cflevelsetmisfit_observation_string_s[j]),7,LevelsetObservationEnum); + element->DatasetInputAdd(StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),cflevelsetmisfit_weights_s[j],inputs,iomodel,cflevelsetmisfit_weights_M_s[j],cflevelsetmisfit_weights_N_s[j],weight_vector_type,StringToEnumx(cflevelsetmisfit_weights_string_s[j]),7,WeightsLevelsetObservationEnum); } } Index: /issm/branches/trunk-dlcheng-ASE/src/m/boundaryconditions/SetMarineIceSheetBC.m =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/m/boundaryconditions/SetMarineIceSheetBC.m (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/m/boundaryconditions/SetMarineIceSheetBC.m (revision 27956) @@ -33,8 +33,14 @@ else %Guess where the ice front is - pos=find(sum(md.mask.ocean_levelset(md.mesh.elements)<0.,2) >0.); - vertexonfloatingice=zeros(md.mesh.numberofvertices,1); - vertexonfloatingice(md.mesh.elements(pos,:))=1.; - vertexonicefront=double(md.mesh.vertexonboundary & vertexonfloatingice); + %pos=find(sum(md.mask.ocean_levelset(md.mesh.elements)<0.,2) >0.); + %vertexonfloatingice=zeros(md.mesh.numberofvertices,1); + %vertexonfloatingice(md.mesh.elements(pos,:))=1.; + %vertexonicefront=double(md.mesh.vertexonboundary & vertexonfloatingice); + + only_ocean_levelset=-double(md.mask.ocean_levelset<=0. & md.mask.ice_levelset>=0.); + only_ocean_levelset_sum=sum(only_ocean_levelset(md.mesh.elements)<0.,2); + pos=find(only_ocean_levelset_sum >0. & only_ocean_levelset_sum < 3.); + vertexonicefront=zeros(md.mesh.numberofvertices,1); + vertexonicefront(md.mesh.elements(pos,:))=1.; end pos=find(md.mesh.vertexonboundary & ~vertexonicefront); @@ -59,8 +65,14 @@ error('mesh type not supported yet'); end -segmentsfront=md.mask.ice_levelset(md.mesh.segments(:,1:numbernodesfront))==0; -segments=find(sum(segmentsfront,2)~=numbernodesfront); +%segmentsfront=md.mask.ice_levelset(md.mesh.segments(:,1:numbernodesfront))==0; +%segments=find(sum(segmentsfront,2)~=numbernodesfront); %Find all nodes for these segments and spc them -pos=md.mesh.segments(segments,1:numbernodesfront); +%pos=md.mesh.segments(segments,1:numbernodesfront); + +numbernodesfront=3; +elementsfront=md.mask.ice_levelset(md.mesh.elements(:,1:numbernodesfront))==0; +elements=find(sum(elementsfront,2)>0); +pos=md.mesh.elements(elements,1:numbernodesfront); + md.stressbalance.spcvx(pos(:))=0; md.stressbalance.spcvy(pos(:))=0; Index: /issm/branches/trunk-dlcheng-ASE/src/m/exp/exptool.m =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/m/exp/exptool.m (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/m/exp/exptool.m (revision 27956) @@ -16,5 +16,5 @@ % - markersize (default=7) % - markeredgecolor (default='r') -% - nofigurecopy (default=0) do not copy current figure, this is needed on some platform to avoid an offset in the figure +% - nofigurecopy (default=0) do not copy current figure, this is needed on some plateform to avoid an offset in the figure % % Usage: @@ -139,5 +139,5 @@ disableDefaultInteractivity(gca); %disables the built-in interactions for the specified axes -%Build backup structure for do and redo +%Build backup structre for do and redo backup=cell(1,3); backup{1,1}=A; Index: /issm/branches/trunk-dlcheng-ASE/src/m/solve/waitonlock.m =================================================================== --- /issm/branches/trunk-dlcheng-ASE/src/m/solve/waitonlock.m (revision 27955) +++ /issm/branches/trunk-dlcheng-ASE/src/m/solve/waitonlock.m (revision 27956) @@ -58,4 +58,6 @@ command = [command ' "[ -f ' lockfilename ' ] && [ -f ' logfilename ' ]" 2>/dev/null']; end + disp(command) + disp(num2str(cluster.port)) end