#########################################
#Jan 2016: Dr. Satkumar Tomer's 1-D model copied
####GW model 1-D GW model function
gw_model <- function(pars, hmin.under.flow, pd.under.flow, rain){
n <- length(rain)
net_recharge <- rf*rain
h <- array(0,n)
discharge <- array(0,n)
time = array(0,n)
h[1] <- hini
time[1] = 1
for (i in 1:(n-1)){
discharge[i] <- sy*(1-pd.under.flow)*(h[i]-hmin.under.flow+ net_recharge[i]/sy)
if (discharge[i]<0){
discharge[i] <- 0
}
h[i+1] <- h[i] + (net_recharge[i]-discharge[i])/sy
time[i+1] = time[i] + 1
}
return(data.frame(time=time,h=h,discharge=discharge))
}
#synthetic forcing
n = 12*20 # temporal length in months
set.seed(123) # to have the repeatibility in results
rain <- rgamma(n, shape=8, rate = 0.1) # synthetic rainfall
#barplot(rain,xlab="Time",ylab="Rainfall (mm)",main="Monthly rainfall")
rain=rain/1000#mm to m of rain
#new
#
#input parameters
rf <- 0.1 # recharge factor, fraction
sy <- 0.005 #specific yield, dimensionless
hini <- 150 # initial condition, metres
hmin.under.flow <- 100
pd.under.flow<-0.967
pars=c("hini"=hini,"rf"=rf,"sy"=sy)##this is the subset of parameters to be SA'd, calibrated, and to inform UA
##
pdf("output.pdf")
output=gw_model(pars,hmin.under.flow,pd.under.flow,rain)
plot(output$h,type='l',xlab='time',ylab='simulated gw level(m)')
dev.off()
##########################
#July 15, 2016
#trial following FME tutorial
#########
library(FME)
#Create dummy observations
#n=240 length of simulated output
#output$h contains the simulated head
#below, observation dataframe is called DataT. It has to(?) contain a "time" and "h" variable like the simulation output
DataT=data.frame(cbind("time"=seq(1,length(n)),"h"=output$h+rnorm(sd=0.5,mean=5,n=240),"sd"=4.5))
DataT2=DataT[,c(1,2)]
##
#Create a function to calculate fit
gwcost=function(pars){
out=gw_model(pars,hmin.under.flow, pd.under.flow, rain)
#cost=modCost(model=out,obs=DataT,err="sd")
cost=modCost(model=out,obs=DataT2)
return(cost=cost)
}
###Above seems to be working(maybe not quite?)
gwcost(pars)
gwcost(pars)$model
#[1] 1075 with no weighting
plot(gwcost(pars)$residuals$res)
#
##Sensitivity using sensFun (here onwards results are not following the JSS article)
##local sensitivity does not seem to be working. does it work only with differential equations?
Sfun=sensFun(gwcost,parms=pars,sensvar=c("h"))
#parameter tinyy to add perturbation? not making any difference
summary(Sfun)
plot(Sfun)
#########