Skip to content
Snippets Groups Projects
Commit 03710fa1 authored by Sumbel Shangareeva's avatar Sumbel Shangareeva
Browse files

rothc-friendly version

parent 2f729192
Branches adjoint
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
plot.py
+ 63
− 4
View file @ 03710fa1
import os
from typing import List
from enum import Enum
import f90nml
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import f90nml
'''
Script for plotting.
......@@ -25,8 +24,8 @@ station = nml['config_namelist']['station_name']
model_data_paths = [
f'results/{carbon_model_type}/{station}_1.txt',
f'results/{carbon_model_type}/{station}_2.txt',
f'results/{carbon_model_type}/{station}_3.txt'
# f'results/{carbon_model_type}/{station}_2.txt',
# f'results/{carbon_model_type}/{station}_3.txt'
]
if (station == 'DAO3') or (station == 'DAO4'):
......@@ -153,6 +152,7 @@ def pools_sums_plots(file_names: List[str], obs_data: pd.DataFrame, out_path: st
color_index = 0
for i, file_name in enumerate(file_names):
df, experiment_name, opt, _ = read_data(file_name)
df = df.drop_duplicates(subset='date', keep="last")
pools_sum = df[df.columns[1:]].sum(axis=1)
ax.plot(df['date'], pools_sum, label=f'Результаты расчетов {experiment_name} {int(opt)}', linewidth=2,
color=colors[color_index])
......@@ -188,21 +188,80 @@ def pools_sums_plots(file_names: List[str], obs_data: pd.DataFrame, out_path: st
plt.close()
def optim_results_plot(data_path: str, obs_data: pd.DataFrame, out_path: str) -> None:
plt.rcParams['font.size'] = '12'
fig, ax = plt.subplots(figsize=(15, 6))
df, experiment_name, opt, _ = read_data(data_path)
df['date'] = df['date'].dt.normalize().apply(lambda d: d.replace(day=1))
df_optim = df.copy().drop_duplicates(subset='date', keep="last")
df_prior = df.copy().drop_duplicates(subset='date', keep="first")
df_optim['model'] = df_optim.iloc[:, 1:].sum(axis=1)
df_prior['model'] = df_prior.iloc[:, 1:].sum(axis=1)
obs_data['date'] = obs_data['date'].apply(lambda d: d.replace(day=1))
obs_data.drop(0, inplace=True)
mp = pd.merge(obs_data, df_prior[['date', 'model']], on='date', how='inner')
mo = pd.merge(obs_data, df_optim[['date', 'model']], on='date', how='inner')
rmse_prior = np.sqrt(((mp[int(opt)] - mp['model']) ** 2).mean())
rmse_opt = np.sqrt(((mo[int(opt)] - mo['model']) ** 2).mean())
ax.plot(df_optim['date'], df_optim['model'],
label=f'After 4D-Var optimization, rmse={rmse_opt:.3f}', linewidth=2,
color='r')
ax.plot(df_prior['date'], df_prior['model'],
label=f'Prior, rmse={rmse_prior:.3f}', linewidth=2,
color='blue')
ax.plot(obs_data['date'], obs_data[int(opt)], 'o', label='Observations')
date_min = df_optim['date'].min() - pd.DateOffset(years=1)
date_max = df_optim['date'].max() + pd.DateOffset(years=1)
ax.xaxis.set_major_locator(mdates.YearLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
ax.set_xlim([date_min, date_max])
ax.set_xlabel('Время, годы')
ax.set_ylabel(f'Содержание углерода в почве, кг/м^2')
ax.grid(True)
plt.tight_layout(pad=4.0)
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), fancybox=True, shadow=True, ncol=2)
plt.xticks(rotation=45)
plt.title(f'Сравнение априорных параметров и калиброванных, {carbon_model_type}, {experiment_name}_{opt}')
plt.savefig(out_path, format='png', bbox_inches='tight')
plt.close()
def main():
obs_data = pd.read_csv(obs_data_path)
obs_data['date'] = pd.to_datetime(obs_data['date'])
obs_data.columns = ['date'] + [i + 1 for i in range(len(obs_data.columns) - 1)]
out_dir = ''
experiment_name = ''
for file_name in model_data_paths:
print(f'processing {file_name} started')
data, experiment_name, opt, out_dir = read_data(file_name)
data = data.drop_duplicates(subset='date', keep="last")
pools_plots(data, f'{out_dir}/{experiment_name}_{opt}')
print(f'pools plots for {file_name} saved at {out_dir}/{experiment_name}_{opt}')
pools_sum_plots(data, experiment_name, opt, f'{out_dir}/{experiment_name}_{opt}_sum.png')
print(f'pools sum plot for {file_name} saved at {out_dir}/')
out_path = f'{out_dir}/{experiment_name}_{opt}_optim.png'
optim_results_plot(file_name, obs_data, out_path)
out_path = f'{out_dir}/{experiment_name}_summary_plot.png'
pools_sums_plots(model_data_paths, obs_data, out_path)
print(f'processing all data is done')
......
source/carbon/carbon_models/carbon_model_ROTHC.f90
+ 17
− 17
View file @ 03710fa1
......@@ -20,11 +20,11 @@ module carbon_model_rothc
call set_pool(pid = n_Catm, name = 'catm')
call set_pool(pid = n_Cveg, name = 'cveg')
call set_pool(pid = n_CDPM, name = 'CDPM', initial_value = 0.01944, alias = CDPM)
call set_pool(pid = n_CRPM, name = 'CRPM', initial_value = 1.73544, alias = CRPM)
call set_pool(pid = n_CBIO, name = 'CBIO', initial_value = 0.06188, alias = CBIO)
call set_pool(pid = n_CHUM, name = 'CHUM', initial_value = 1.0416, alias = CHUM)
call set_pool(pid = n_CIOM, name = 'CIOM', initial_value = 0.245, alias = CIOM)
call set_pool(pid = n_CDPM, name = 'CDPM', initial_value = 0.04936, alias = CDPM)
call set_pool(pid = n_CRPM, name = 'CRPM', initial_value = 0.799, alias = CRPM)
call set_pool(pid = n_CBIO, name = 'CBIO', initial_value = 0.11873, alias = CBIO)
call set_pool(pid = n_CHUM, name = 'CHUM', initial_value = 7.07025, alias = CHUM)
call set_pool(pid = n_CIOM, name = 'CIOM', initial_value = 0.814, alias = CIOM)
! station: opt: CDPM: CRPM: CBIO: CHUM: CIOM:
! DAO4 1 0.01944 1.73544 0.06188 1.0416 0.245
......@@ -74,7 +74,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(1))
call set_mult(n_F, 'lin', x_ijn = CDPM)
call set_mult(n_F, 'lin', x_ijn = CDPM, pool_num=n_CDPM)
call set_flux(fid = n_F, pid_out = n_CRPM, pid_in = n_Catm, name = 'RRPM')
......@@ -89,7 +89,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(2))
call set_mult(n_F, 'lin', x_ijn = CRPM)
call set_mult(n_F, 'lin', x_ijn = CRPM, pool_num=n_CRPM)
call set_flux(fid = n_F, pid_out = n_CBIO, pid_in = n_Catm, name = 'RBIO')
......@@ -104,7 +104,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(3))
call set_mult(n_F, 'lin', x_ijn = CBIO)
call set_mult(n_F, 'lin', x_ijn = CBIO, pool_num=n_CBIO)
call set_flux(fid = n_F, pid_out = n_CHUM, pid_in = n_Catm, name = 'RHUM')
......@@ -119,7 +119,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(4))
call set_mult(n_F, 'lin', x_ijn = CHUM)
call set_mult(n_F, 'lin', x_ijn = CHUM, pool_num=n_CHUM)
! atmosphere fall
......@@ -137,7 +137,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(1))
call set_mult(n_F, 'lin', x_ijn = CDPM)
call set_mult(n_F, 'lin', x_ijn = CDPM, pool_num=n_CDPM)
call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CBIO, name = 'BIO_prop*RRPM')
......@@ -154,7 +154,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(2))
call set_mult(n_F, 'lin', x_ijn = CRPM)
call set_mult(n_F, 'lin', x_ijn = CRPM, pool_num=n_CRPM)
call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CBIO, name = 'BIO_prop*RBIO')
......@@ -171,7 +171,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(3))
call set_mult(n_F, 'lin', x_ijn = CBIO)
call set_mult(n_F, 'lin', x_ijn = CBIO, pool_num=n_CBIO)
call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CBIO, name = 'BIO_prop*RHUM')
......@@ -188,7 +188,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(4))
call set_mult(n_F, 'lin', x_ijn = CHUM)
call set_mult(n_F, 'lin', x_ijn = CHUM, pool_num=n_CHUM)
call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CHUM, name = 'HUM_prop*RDPM')
......@@ -205,7 +205,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(1))
call set_mult(n_F, 'lin', x_ijn = CDPM)
call set_mult(n_F, 'lin', x_ijn = CDPM, pool_num=n_CDPM)
call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CHUM, name = 'HUM_prop*RRPM')
......@@ -222,7 +222,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(2))
call set_mult(n_F, 'lin', x_ijn = CRPM)
call set_mult(n_F, 'lin', x_ijn = CRPM, pool_num=n_CRPM)
call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CHUM, name = 'HUM_prop*RBIO')
......@@ -239,7 +239,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(3))
call set_mult(n_F, 'lin', x_ijn = CBIO)
call set_mult(n_F, 'lin', x_ijn = CBIO, pool_num=n_CBIO)
call set_flux(fid = n_F, pid_out = n_Catm, pid_in = n_CHUM, name = 'HUM_prop*RHUM')
......@@ -256,7 +256,7 @@ module carbon_model_rothc
call set_mult(n_F, 'lin', y0 = 0.6, k = -0.4, x_ij = veg, x0 = 1.)
call set_mult(n_F, 'lin', x = ks(4))
call set_mult(n_F, 'lin', x_ijn = CHUM)
call set_mult(n_F, 'lin', x_ijn = CHUM, pool_num=n_CHUM)
end subroutine carbon_model_assembly
......
source/carbon/carbon_models/carbon_model_ROTHC_aux.f90
+ 2
− 0
View file @ 03710fa1
......@@ -192,6 +192,8 @@ contains
end subroutine
subroutine carbon_model_deallocate()
deallocate(s0, smin, SA, CO2_prop, BIO_prop, HUM_prop)
fdpm = 0.
end subroutine carbon_model_deallocate
......
source/main_4d_var.f90
+ 20
− 7
View file @ 03710fa1
......@@ -15,16 +15,29 @@ program main
implicit none
real, allocatable :: cost_hist(:)
real, allocatable :: lr_param(:)
real, allocatable :: B_inv_diag(:)
integer :: idx
call config_init()
nparam = 2
allocate(opt_pool_num(nparam), theta_b(nparam), theta_opt(nparam), B_inv(nparam, nparam))
opt_pool_num = [3, 4] !< id of pools to optimize
theta_b = [ 0., 6. ] !< a priori parameters
B_inv = reshape([1.0,0.0,0.0,0.1],[nparam,nparam]) !< error covariance of a priori
nparam = 4
allocate(opt_pool_num(nparam), theta_b(nparam), theta_opt(nparam), B_inv(nparam, nparam), lr_param(nparam), B_inv_diag(nparam))
! [3, 4, 5, 6], [3,4]
opt_pool_num = [3, 4, 5, 6] !< id of pools to optimize
! [.01, 0.1, .05, 5.], [.0, 6.]
theta_b = [.01, 0.1, .05, 5.] !< a priori parameters
! [1., 1., 5., 5.], [5., 5.]
lr_param = [1., 1., 5., 5.]
! [1., 1., 1.,.1], [1., .1]
B_inv_diag = [1., 1., 1.,.1] !< error covariance of a priori
B_inv = 0.
do idx = 1, nparam
B_inv(idx, idx) = B_inv_diag(idx)
end do
R_inv = 1.0/0.1
call grid_init()
......@@ -53,7 +66,7 @@ program main
! 4D-Var
! ---------------------------------------------------------------------------------
call fourd_var(n_iter = 20, lr = 2.5, tol = 1.0e-4, cost_hist = cost_hist)
call fourd_var(n_iter=20, lr_param=lr_param, tol=1.0e-6, cost_hist=cost_hist)
write(*,'("Оптимальные IC: ",2F12.6)') theta_opt
......
source/optimization/carbon_obs_operator.f90
+ 10
− 6
View file @ 03710fa1
......@@ -56,9 +56,11 @@ contains
residual = 0.
do m = 1, nobs
if (obs_idx(m) == t-1) then
do k = 1, size(opt_pool_num)
residual = residual + state_ts(opt_pool_num(k),ii,jj,nn,t-1)
end do
residual = residual + state_ts(3,ii,jj,nn,t-1) + state_ts(4,ii,jj,nn,t-1)
residual = residual + state_ts(5,ii,jj,nn,t-1) + state_ts(6,ii,jj,nn,t-1)
residual = residual + state_ts(7,ii,jj,nn,t-1)
residual = (residual - obs_val(m))/yrs
do k = 1, size(opt_pool_num)
forcing(opt_pool_num(k)) = residual * R_inv
......@@ -79,9 +81,11 @@ contains
do k = 1, nobs
model_obs_diff = 0.
do m = 1, size(opt_pool_num)
model_obs_diff = model_obs_diff + state_ts(opt_pool_num(m),i0,j0,1,obs_idx(k))
end do
model_obs_diff = model_obs_diff + state_ts(3,i0,j0,1,obs_idx(k)) + state_ts(4,i0,j0,1,obs_idx(k))
model_obs_diff = model_obs_diff + state_ts(5,i0,j0,1,obs_idx(k)) + state_ts(6,i0,j0,1,obs_idx(k))
model_obs_diff = model_obs_diff + state_ts(7,i0,j0,1,obs_idx(k))
model_obs_diff = model_obs_diff - obs_val(k)
cost_o = cost_o + 0.5*(model_obs_diff)**2 * R_inv
end do
......
source/optimization/carbon_optimize.f90
+ 4
− 3
View file @ 03710fa1
......@@ -17,12 +17,13 @@ module carbon_optimize
contains
subroutine fourd_var(n_iter, lr, tol, cost_hist)
subroutine fourd_var(n_iter, lr_param, tol, cost_hist)
use grid, only : i0,i1,j0,j1, ntime
integer, intent(in) :: n_iter !< number of iterations of opt algo
real, intent(in) :: lr, tol !< learning rate for Adagrad
real, intent(in) :: tol
real, intent(in) :: lr_param(nparam) ! ← learning rate for Adagrad
real, allocatable, intent(out) :: cost_hist(:)
integer :: k, it
......@@ -74,7 +75,7 @@ subroutine fourd_var(n_iter, lr, tol, cost_hist)
! ---------------------------------------------------------------------------------
do k = 1, nparam
cache(k) = cache(k) + grad(k)**2
theta_new(k) = theta(k) - lr*grad(k)/sqrt(cache(k)+eps)
theta_new(k) = theta(k) - lr_param(k)*grad(k)/sqrt(cache(k)+eps)
end do
print*, 'Niter: ', it, 'J: ', cost, 'θ: ', theta, 'θ_new: ', theta_new, 'grad: ', grad
......
ui1_config.nml
+ 1
− 1
View file @ 03710fa1
......@@ -12,7 +12,7 @@
! -----------------------------------------------------------------
!> Углеродная модель:
carbon_model_type = 'socs'
carbon_model_type = 'rothc'
!
! 'inmcm' - модель inmc [1,2]
! 'socs' - модель SOCS [3]
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment