Demo entry 5661002

model

   

Submitted by anonymous on Jul 07, 2016 at 14:15
Language: Python 3. Code size: 8.4 kB.

#!/usr/bin/python3.5
#-*- coding: utf-8 -*-

'''
###############
# Modèle WASP #
###############
Dernière mise à jour: 01-07-2016.\n
Réalisé par Julien Diot, julien.diot@supagro.fr\n
Dans le cadre d'un stage ingénieur agronome 2eme année (Montpellier SupAgro).
'''

## Importation des modules:
import os #fonctions pour le changement du répertoire courant.
import matplotlib.pyplot as plt  # fonctions graphiques
import time # fonctions de calcul temporel


#modules personnels:
os.chdir("C:/Users\Julien\OneDrive\Documents\SupaAgro\Cours\Stage 2A\Scripts\FINAL\Modèle")# définition du répertoir contenant les modules
import fonctions as fun #importation des fonctions du modèle
import parametres as p  #importation des paramètres du modèle


## Initialisation des variables:
timeT= 30   # durée de la modélisation (en jours)
dt= 60      # pas de temps de la méthode Runge-Kutta (en s)
start="02/05/2001" # date de début de modélisation (à accorder avec les données d'entrées)
date=time.localtime()# récupération de la date et de l'heure pour l'écriture dans le fichier

#initialisation des variables d'états. (Valeurs provenant de la station 2: moyennes des valeures de avril et mai 2001)
LC_do=  [9.975]         #(mgO2/L)
LC_p=   [0.011815*12.5] #(mgC/L)
LC_cbod=[4.52]          #(mgC/L)
LC_NH4= [0.0935]        #(mgN/L)
LC_NO3= [0.4405]        #(mgN/L)
LC_ip=  [0.0015]        #(mgP/L)
LC_on=  [1.781]         #(mgN/L)
LC_op=  [0.068]         #(mgP/L)
LC_chl= [0.011815]      #(mgC/L)
Ltemps=[0]

for i in range(1,int(timeT*24*3600/dt),1):
    Ltemps.append(None)
    LC_p.append(None)
    LC_do.append(None)
    LC_cbod.append(None)
    LC_NH4.append(None)
    LC_NO3.append(None)
    LC_ip.append(None)
    LC_on.append(None)
    LC_op.append(None)
    LC_chl.append(None)


## Importation des entrée:
#chemin des différents fichiers contenant les entrées.
fileT_w="C:/Users\Julien\OneDrive\Documents\SupaAgro\Cours\Stage 2A\Scripts\FINAL\Modèle\T_water_2001_2002_station2.txt"
fileT_a="C:/Users\Julien\OneDrive\Documents\SupaAgro\Cours\Stage 2A\Scripts\FINAL\Modèle\Wuxi 2001-2002_Ta.txt"
fileW="C:/Users\Julien\OneDrive\Documents\SupaAgro\Cours\Stage 2A\Scripts\FINAL\Modèle\Wuxi 2001-2002_wind.txt"
filesun="C:/Users\Julien\OneDrive\Documents\SupaAgro\Cours\Stage 2A\Scripts\FINAL\Modèle\Wuxi 2001-2002_sun.txt"
fileI="C:/Users\Julien\OneDrive\Documents\SupaAgro\Cours\Stage 2A\Scripts\FINAL\Modèle\I_t2001-2003.txt"
#fileS= ???

#Interpolation des fichiers d'entrées:
fT_w=fun.entryImport(fileT_w,start,timeT,sep=";",header=True)
fT_a=fun.entryImport(fileT_a,start,timeT,sep=";",header=True)
fW=fun.entryImport(fileW,start,timeT,sep=";",header=True)
fsun=fun.entryImport(filesun,start,timeT,sep=";",header=True)
fI=fun.entryImport_CtDay(fileI,start,timeT,sep=";",header=True)
#fS=fun.entryImport(fileS,start,timeT,sep=";",header=True)

#Entrées sans fichiers:
# T_w=
# T_a=
# W=
# t_U=
# t_D=
# f=
# I=
# I_a=
S=0

## Algorithme:
tic=time.time() #lancement du chronomètre.
num=0      #initiation du compteur de calculs (pour l'affichage du temps de modélisation restant)
coefRK=[float('inf'),2,2,1] #liste de coéfficients pour la méthode Runge-Kutta.
listODE=[fun.S_do,fun.S_p,fun.S_cbod,fun.S_NH4,fun.S_NO3,fun.S_ip,fun.S_on,fun.S_op] #liste des ODE des variables d'état.

for i in range(1,int(timeT*24*3600/dt),1): #int(timeT*24*3600/dt) = convertion en nb d'itération
    t=i*dt # calcule de la variable de temps
    num+=1 #incrémentation du compteur
    
    C_do=LC_do[i-1]
    C_p= LC_p[i-1]
    C_cbod=LC_cbod[i-1]
    C_NH4= LC_NH4[i-1]
    C_NO3= LC_NO3[i-1]
    C_ip=  LC_ip[i-1]
    C_on=  LC_on[i-1]
    C_op=  LC_op[i-1]
    C_chl= LC_chl[i-1]
    
    #Méthode de Runge-Kutta d'ordre 4:
    k=[[0 for i in range(8)] for i in range(5)] # initialisation de la liste des variations élémentaires
    
    for order in range(4):
        t2=t+dt/coefRK[order] # positionnement dans le temps.
        
        #récupération des entrées au temps adéquate.
        #il faut commenter ce paragraphe si les entrées sont constantes.
        T_w=fT_w(t2) 
        T_a=fT_a(t2)
        W=fW(t2)
        sun=fsun(t2)
        I=fI(t2)
        
        t_U= 12-(sun/2)
        t_D= 12+(sun/2)
        f=sun/24
        if f!=0:
            I_a= 0.9*I/f
        else: I_a=0
        
        
        #création de la liste d'arguments pour les ODE.
        X=[t2,C_do,C_p,C_cbod,C_NH4,C_NO3,C_ip,C_on,C_op,C_chl,T_w,T_a,W,t_U,t_D,f,I,I_a,S]
        var=1
        for ode in listODE:
            X[var]+=k[order][var-1]/coefRK[order]
            k[order+1][var-1]=(dt*ode(*X)/(24*3600))
            var+=1
    
    #Attribution des nouvelles valeurs pour les variables d'état.
    #(avec test de négativité)
    Ltemps[i]=t/3600/24 #convertion en jour
    
    LC_do[i] = LC_do[i-1]+k[1][0]/6+k[2][0]/3+k[3][0]/3+k[4][0]/6
    if LC_do[i]<0:
        LC_do[i]=0
    LC_p[i] = LC_p[i-1]+k[1][1]/6+k[2][1]/3+k[3][1]/3+k[4][1]/6
    if LC_p[i]<0:
        LC_p[i]=0
    LC_cbod[i] = LC_cbod[i-1]+k[1][2]/6+k[2][2]/3+k[3][2]/3+k[4][2]/6
    if LC_cbod[i]<0:
        LC_cbod[i]=0
    LC_NH4[i] = LC_NH4[i-1]+k[1][3]/6+k[2][3]/3+k[3][3]/3+k[4][3]/6
    if LC_NH4[i]<0:
        LC_NH4[i]=0
    LC_NO3[i] = LC_NO3[i-1]+k[1][4]/6+k[2][4]/3+k[3][4]/3+k[4][4]/6
    if LC_NO3[i]<0:
        LC_NO3[i]=0
    LC_ip[i] = LC_ip[i-1]+k[1][5]/6+k[2][5]/3+k[3][5]/3+k[4][5]/6
    if LC_ip[i]<0:
        LC_ip[i]=0
    LC_on[i] = LC_on[i-1]+k[1][6]/6+k[2][6]/3+k[3][6]/3+k[4][6]/6
    if LC_on[i]<0:
        LC_on[i]=0
    LC_op[i] = LC_op[i-1]+k[1][7]/6+k[2][7]/3+k[3][7]/3+k[4][7]/6
    if LC_op[i]<0:
        LC_op[i]=0
    LC_chl[i] = LC_p[i]/12.5

    #Affichage du temps restant avant la fin de la modélisation (toutes les 5 secondes).
    if time.time()-tic>5 or i==1000:
        vit=num/(time.time()-tic)        #calcul de la vitesse de calcul. à partir du compteur "num"
        num=0       #réinitialisation du compteur.
        tic=time.time()      #réinitialisation du chronomètre.
        tRestant=((timeT*24*3600/dt)-i)/vit     #calcul du temps restant. 
        avancement=i/int(timeT*24*3600/dt)*100  #Calcul de l'avancement de la modélisation (en %)
        #affichage.
        print(int(round(avancement,1)),"% -- Encore environ ",int(tRestant/60)," minutes et ",int(tRestant%60)," secondes")

##Affichage graphique:
n=len(Ltemps) # changer cette variable pour n'afficher qu'une partie de la modélisation
step=1     # pas de l'affichage 

plt.plot(Ltemps[0:n:step],LC_do[0:n:step],label="DO")
plt.legend(loc='best')
 
plt.plot(Ltemps[0:n:step],LC_p[0:n:step],label="Phyto")
plt.legend(loc='best')

plt.plot(Ltemps[0:n:step],LC_cbod[0:n:step],label="CBOD")
plt.legend(loc='best')

plt.plot(Ltemps[0:n:step],LC_NH4[0:n:step],label="NH4")
plt.legend(loc='best')

plt.plot(Ltemps[0:n:step],LC_NO3[0:n:step],label="NO3")
plt.legend(loc='best')

plt.plot(Ltemps[0:n:step],LC_ip[0:n:step],label="IP")
plt.legend(loc='best')

plt.plot(Ltemps[0:n:step],LC_on[0:n:step],label="ON")
plt.legend(loc='best')

plt.plot(Ltemps[0:n:step],LC_op[0:n:step],label="OP")
plt.legend(loc='best')

plt.plot(Ltemps[0:n:step],LC_chl[0:n:step],label="chla")
plt.legend(loc='best')

plt.show()

##Ecriture des résultats dans un fichier texte:
#Récupération des information sur la modélisation:
strDate=str(date.tm_year)+"-"+str(date.tm_mon)+"-"+str(date.tm_mday)+"_"+str(date.tm_hour)+"-"+str(date.tm_min)
info="durée de modélisation: "+str(timeT)+" jours\nPas de temps utilisé: "+str(dt)+" secondes\nModélisation à partir du "+start+"\n"

#écriture des informations
nameFile="WASP_results_"+strDate+".txt"
file= open(nameFile, "w")
file.write("Modélisation lancée le"+str(date.tm_mday)+"/"+str(date.tm_mon)+"/"+str(date.tm_year)+" à "+str(date.tm_hour)+"h"+str(date.tm_min)+"min\n"+info)

#écriture des résultats
file.write("temps:"+str(Ltemps)+"\n")
file.write("DO:"+str(LC_do)+"\n")
file.write("Phyto:"+str(LC_p)+"\n")
file.write("CBOD:"+str(LC_cbod)+"\n")
file.write("NH4:"+str(LC_NH4)+"\n")
file.write("NO3:"+str(LC_NO3)+"\n")
file.write("IP:"+str(LC_ip)+"\n")
file.write("ON:"+str(LC_on)+"\n")
file.write("OP:"+str(LC_op)+"\n")
file.write("chla:"+str(LC_chl)+"\n")

file.close()

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).