From 905858a00b1a31dd960ac57320e58468295da07a Mon Sep 17 00:00:00 2001 From: rois1995 <43813346+rois1995@users.noreply.github.com> Date: Sun, 9 Dec 2018 16:16:46 +0100 Subject: [PATCH] Add files via upload Ancora non funzionano lo script Boundary_conditions.f90 --- Boundary_conditions.f90 | 116 +++ LOGO.f90 | 63 ++ Makefile | 141 +++ RK.f90 | 207 ++++ SW_completo.f90 | 171 ++++ conversione.f90 | 72 ++ fix_BC.f90 | 121 +++ fluxes.f90 | 98 ++ gnufor_mio.f90 | 2029 +++++++++++++++++++++++++++++++++++++++ inversa.f90 | 94 ++ problema_vero.f90 | 192 ++++ slope_calc.f90 | 131 +++ slope_det.f90 | 49 + 13 files changed, 3484 insertions(+) create mode 100644 Boundary_conditions.f90 create mode 100644 LOGO.f90 create mode 100644 Makefile create mode 100644 RK.f90 create mode 100644 SW_completo.f90 create mode 100644 conversione.f90 create mode 100644 fix_BC.f90 create mode 100644 fluxes.f90 create mode 100644 gnufor_mio.f90 create mode 100644 inversa.f90 create mode 100644 problema_vero.f90 create mode 100644 slope_calc.f90 create mode 100644 slope_det.f90 diff --git a/Boundary_conditions.f90 b/Boundary_conditions.f90 new file mode 100644 index 0000000..9cf40be --- /dev/null +++ b/Boundary_conditions.f90 @@ -0,0 +1,116 @@ +module Boundary_conditions + + ! Verifica delle condizioni al contorno da applicare + + use problema_vero + + use conversione + + implicit none + +contains + + subroutine BC_dx(stato_dx, BC) + + real(kind = 8), dimension(:), intent(in) :: stato_dx + + real(kind = 8), dimension(size(stato_dx), size(stato_dx)), intent(inout) :: BC + + real(kind = 8), dimension(size(stato_dx)) :: lambda + + real(kind = 8), dimension(size(stato_dx), size(stato_dx)) :: R, L + + real(kind = 8), dimension(size(stato_dx)) :: stato_ext,& ! stato_ext = [rho_ext, m_ext] + delta_w,& ! incremento delle variabili conservative + delta_v,& ! incremento variabili caratteristiche + BC_cons ! BC in forma conservativa + + + call eigenstructure(stato_dx, lambda, R, L) + + !write(*,*) " BC_fis(:, 2) prima di trasformazioni stato dx = ", BC(:, 2) + + call fis2cons(BC(:, 2), BC_cons) ! Calcolo la m all' uscita dell'estremo destro + ! BC = [ altezza , velocità ] + + delta_w = BC_cons - stato_dx + + delta_v = matmul(L, delta_w) ! passo alle variabili caratteristiche + + delta_v = delta_v * 0.5d0 * (1 - sign(1.d0, lambda)) ! controllo sul segno delle velocità di propagazione + ! applico condizioni al contorno solo se + ! è minore o uguale a 0 + + delta_w = matmul(R, delta_v) ! torno alle variabili conservative + + BC_cons = stato_dx + delta_w + + call cons2fis(BC_cons, BC(:, 2)) ! Torno alle variabili fisiche + ! per le condizioni al contorno + ! BC(1, 2) = h + + !write(*,*) " BC_fis(:, 2) dopo trasformazioni stato dx = ", BC(:, 2) + + + endsubroutine BC_dx + + subroutine BC_sx(stato_sx, BC) + + real(kind = 8), dimension(:), intent(in) :: stato_sx + + real(kind = 8), dimension(size(stato_sx), size(stato_sx)), intent(inout) :: BC + + real(kind = 8), dimension(size(stato_sx)) :: lambda + + real(kind = 8), dimension(size(stato_sx), size(stato_sx)) :: R, L + + real(kind = 8), dimension(size(stato_sx)) :: stato_ext,& ! stato_ext = [rho_ext, m_ext] + delta_w,& ! incremento delle variabili conservative + delta_v,& ! incremento variabili caratteristiche + BC_cons ! BC in forma conservativa + + + call eigenstructure(stato_sx, lambda, R, L) + + write(*,*) " BC_fis(:, 1) prima di trasformazioni stato sx = ", BC(:, 1) + + + call fis2cons(BC(:, 1), BC_cons) ! Calcolo la m all' uscita dell'estremo sinistro + ! BC = [ velocità, pressione ] + + delta_w = BC_cons - stato_sx + + !write(*,*) "R = ", R + + !write(*,*) "prima del checking lambda delta_w = ", delta_w + + delta_v = matmul(L, delta_w) ! passo alle variabili caratteristiche + + !write(*,*) "L = ", L + + !write(*,*) "stato_sx = ", stato_sx + + !write(*,*) "lambda_sx = ", lambda + + where (lambda <= 0) ! controllo sul segno delle velocità di propagazione + + delta_v = 0 ! applico condizioni al contorno solo se è minore o uguale a 0 + + endwhere + + delta_w = matmul(R, delta_v) ! torno alle variabili conservative + + BC_cons = stato_sx + delta_w + + !write(*,*) "delta_w = ", delta_w + + call cons2fis(BC_cons, BC(:, 1)) ! Torno alle variabili fisiche + ! per le condizioni al contorno + ! BC(1, 1) = u + + !write(*,*) " BC_fis(:, 1) dopo trasformazioni stato sx = ", BC(:, 1) + + + endsubroutine BC_sx + +endmodule Boundary_conditions diff --git a/LOGO.f90 b/LOGO.f90 new file mode 100644 index 0000000..b25527d --- /dev/null +++ b/LOGO.f90 @@ -0,0 +1,63 @@ +module LOGO + +implicit none + +contains + subroutine disegna_logo() + + +write(*,*) "------------------------------------BENVENUTI NEL FOTTUTISSIMO RI",& + "SOLUTORE DELLE SW-EQUATIONS-----------------------------------" + +write(*,*) " ", & + " ------------------------***---------------------------" +write(*,*) " ", & + " -------------------*************----------------------" +write(*,*) " __ __ ___ ___ ___ _ _ _ _ __ __ ", & + " ----------------********************------------------" +write(*,*) "|__ | | | | | | | | |_ |_ | | | | | || | ", & + " ------------**************************----------------" +write(*,*) "| |__| | | |__| | | _| _| | | | |__||__| ", & + " ---------*******************************--------------" +write(*,*) " ", & + " -- ________******************________*****------------" +write(*,*) " ", & + " --/ \****************/ \*****-----------" +write(*,*) " ", & + " -/ \**************/ \*****----------" +write(*,*) " ", & + " -| . |--------------| . |*****----------" +write(*,*) " ", & + " -\ /******/\******\ /*****----------" +write(*,*) " ", & + " --\________/******/++\******\________/******----------" +write(*,*) " ", & + " --***************/++++\*********************----------" +write(*,*) " ", & + " --**************/++++++|*******************-----------" +write(*,*) " ", & + " --*************/+++++++|*******************-----------" +write(*,*) " ", & + " --************( . )( . )******************------------" +write(*,*) " ", & + " --***********_____**********************--------------" +write(*,*) " ", & + " --**********/ \********************---------------" +write(*,*) " ", & + " ---********/ \******************----------------" +write(*,*) " ", & + " ---********| O |*****************-----------------" +write(*,*) " ", & + " ---********\ /***************-------------------" +write(*,*) " ", & + " ---*********\_____/**************---------------------" +write(*,*) " ", & + " ----***************************-----------------------" +write(*,*) " ", & + " ------********************----------------------------" +write(*,*) " ", & + " ----------*******-------------------------------------" + +endsubroutine disegna_logo + +endmodule LOGO diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..55e983d --- /dev/null +++ b/Makefile @@ -0,0 +1,141 @@ +############################################################ +# MAKEFILE FOR F95 # +############################################################ + +# Zero axiom: the file name of a makefile must be mandatorily +# +# makefile OR Makefile + +# To start the execution, type in terminal: make + +# The names (identifictars of variables) for aliases +# are completely free! They are usually in CAPITALS. + +# FUNDAMENTAL WARNINGS +# +# The actions, typically $(CMP) and $(LNK), must be +# NECESSARILY preceded by a DAMNED! TAB character +# + +############################################################ +# +# Compiler and linker +# +# The following definitions are for Fortran 95 only, +# but they can extended if another language is used, +# or modified when other compilers are considered. + +# LNK alias to be used as $(LNK) to compile and link main files .f90 +# CMP alias to be used as $(CMP) to compile the source files .f90 +# OPT alias to be used as $(OPT) to select the compilation options + +# NAG Compiler +# ============ +# +CMP = f95 -c +LNK = f95 +OPT = -g -mcmodel=large +# OPT = -gline # Line number is indicated for errors + # occurring in runtime (for debugging) +# OPT = -O3 # Execution optimization + + +# Intel Compiler +# ============== +# +# CMP = ifort -c -heap-arrays 100 +# LNK = ifort +# OPT = +# OPT = -g -DD -debug -traceback -fpe0 -check all +# OPT = -fast +# OPT = -O3 + + +# GNU Fortran Compiler +# =================== +# +# CMP = gfortran -c +# LNK = gfortran +# OPT = -fno-automatic -fbounds-check -ffpe-trap=invalid,zero,overflow -ggdb3 +# OPT = -fno-automatic -O3 + +############################################################ +# +# Objects: list of all objects *.o + +OBJS = gnufor_mio.o LOGO.o conversione.o inversa.o problema_vero.o fluxes.o\ + Boundary_conditions.o #fix_BC.o +# OBJS alias to be used as $(OBJS) + +# Alternative way for line feed +#OBJS = Riemann_solvers.o \ +# polytropic_ideal_gas.o + +############################################################ +# Executable generation by the linker / the effect of command +# / make starts from here + +finale.exe: SW_completo.o $(OBJS) + $(LNK) $(OPT) SW_completo.o $(OBJS) \ + -o finale.exe # Name chosen for .exe + +############################################################ +# Objects generation + +SW_completo.o: SW_completo.f90 $(OBJS) + $(CMP) $(OPT) SW_completo.f90 + +gnufor_mio.o: gnufor_mio.f90 + $(CMP) $(OPT) gnufor_mio.f90 + +LOGO.o: LOGO.f90 + $(CMP) $(OPT) LOGO.f90 + +conversione.o: conversione.f90 + $(CMP) $(OPT) conversione.f90 + +inversa.o: inversa.f90 + $(CMP) $(OPT) inversa.f90 + +problema_vero.o: problema_vero.f90 inversa.o + $(CMP) $(OPT) problema_vero.f90 + +Boundary_conditions.o: Boundary_conditions.f90 problema_vero.o conversione.o + $(CMP) $(OPT) Boundary_conditions.f90 + +#fix_BC.o: fix_BC.f90 conversione.o +# $(CMP) $(OPT) fix_BC.f90 + +fluxes.o: fluxes.f90 problema_vero.o conversione.o Boundary_conditions.o + $(CMP) $(OPT) fluxes.f90 + + + +#plot.o: plot.f90 checking.o compl_Newton.o gnufor_mio.o +# $(CMP) $(OPT) plot.f90 + +# +# il carattere TAB c'e' ma non appare a causa di # +# provare a fare return appena dopo # per credere +# il TAB (con il suo fottutissimo doppio spazio) resuscita + + +############################################################ +# Cleaning command to remove all objects *.o, files *.mod and +# the executables *.exe + +clean: + @echo cleaning objects, modules, executables + rm *.o *.mod *.exe + +cleanimages: + @echo cleaning images files + rm video/plot_video/* + +cleanvideo: + @echo cleaning video files + rm plot_video/* + +# The command "make clean" deletes all the indicated files + +############################################################ diff --git a/RK.f90 b/RK.f90 new file mode 100644 index 0000000..863355e --- /dev/null +++ b/RK.f90 @@ -0,0 +1,207 @@ +module RK + + ! Integratore nel tempo tramite Runge-Kutta a 2 ed a 3 passi + + use fluxes + + use slope_det + + implicit none + + +contains + + function cose(stato, k, flag, lambda, sslope, BC_fis, V, R, Dx_vec) result(U_next) + + real(kind = 8), dimension(:, :, :), intent(inout) :: R, sslope + + real(kind = 8), dimension(:, :), intent(inout) :: stato, V, lambda + + real(kind = 8), intent(in) :: k + + real(kind = 8), dimension(:), intent(in) :: Dx_vec + + integer, intent(in) :: flag + + real(kind = 8), dimension(:, :), intent(inout) :: BC_fis + + real(kind = 8), dimension(size(BC_fis,1), size(BC_fis,2)) :: BC_cons + + real(kind = 8), dimension(size(stato, 1), size(stato, 2)) :: U_next + + real(kind = 8), dimension(size(stato, 1), size(stato, 2)) :: U_app_1, U_app_2,& + U_app_3, U_app_4,& + U_app_1_2, U_app_2_2,& + U_app_3_2, U_app_4_2,& + d_flux, Dx_mat + + real(kind = 8), dimension(size(stato, 1), size(stato, 2) + 1) :: flusso + + real(kind = 8) :: a = 1.d0 + + integer :: j + + ! Calcolo delle pendenze + + call slope_fin(Dx_vec, stato, BC_fis, sslope, BC_cons, V, R) + + ! Calcolo del vettore differenza dei flussi + + call fluxo(stato, k, Dx_vec, sslope, BC_cons, V, R, flusso) + + d_flux = flusso(:, 2:size(flusso, 2)) - flusso(:, 1:size(flusso, 2) - 1) + + do j = 1, size(Dx_mat, 1) + + Dx_mat(j, :) = Dx_vec + + enddo + + + select case (flag) + + case (1) + + write(*,*) "-----------------------------------------------" + write(*,*) "No Runge-Kutta" + write(*,*) "-----------------------------------------------" + + U_next = stato - k * d_flux / Dx_mat + + + case (2) + + write(*,*) "-----------------------------------------------" + write(*,*) "Esecuzione del metodo di Runge-Kutta a 2 passi" + write(*,*) "-----------------------------------------------" + + ! Primo passo + + U_app_1 = stato - 0.5d0 * k * d_flux / Dx_mat + + ! Calcolo delle pendenze + + call slope_fin(Dx_vec, U_app_1, BC_fis, sslope, BC_cons, V, R) + + ! Ricalcolo i flussi + + call fluxo(U_app_1, k, Dx_vec, sslope, BC_cons, V, R, flusso) + + d_flux = flusso(:, 2:size(flusso, 2)) - flusso(:, 1:size(flusso, 2) - 1) + + ! Secondo passo + + U_next = stato - 0.5d0 * k * d_flux / Dx_mat + + + case (3) + + write(*,*) "-----------------------------------------------" + write(*,*) "Esecuzione del metodo di Runge-Kutta a 3 passi" + write(*,*) "-----------------------------------------------" + + ! Primo passo + + U_app_1 = stato - 0.5d0 * k * d_flux / Dx_mat + + ! Calcolo delle pendenze + + call slope_fin(Dx_vec, U_app_1, BC_fis, sslope, BC_cons, V, R) + + ! Ricalcolo i flussi + + call fluxo(U_app_1, k, Dx_vec, sslope, BC_cons, V, R, flusso) + + d_flux = flusso(:, 2:size(flusso, 2)) - flusso(:, 1:size(flusso, 2) - 1) + + ! Secondo passo + + U_app_2 = (3/4.d0) * stato + 0.25 * U_app_1 - (k / 8.d0) * d_flux / Dx_mat + + + ! Calcolo delle pendenze + + call slope_fin(Dx_vec, U_app_2, BC_fis, sslope, BC_cons, V, R) + + ! Ricalcolo i flussi + + call fluxo(U_app_2, k, Dx_vec, sslope, BC_cons, V, R, flusso) + + d_flux = flusso(:, 2:size(flusso, 2)) - flusso(:, 1:size(flusso, 2) - 1) + + ! Secondo passo + + U_next = (1/3.d0) * stato + (2/3.d0) * U_app_2 - (1/3.d0) * k * d_flux / Dx_mat + + !!!!!!!!!!!!!!!!!!!!!!!!!! Controlloare se è giusto mettere i /Dx_mat anche in RK3 e RK4 + + + case (4) + + write(*,*) "-----------------------------------------------" + write(*,*) "Esecuzione del metodo di Runge-Kutta a 4 passi" + write(*,*) "-----------------------------------------------" + + ! Primo passo + + U_app_1 = - 0.5d0 * k * d_flux / Dx_mat + + U_app_1_2 = U_app_1 / 2.d0 + stato + + ! Calcolo delle pendenze + + call slope_fin(Dx_vec, U_app_1_2, BC_fis, sslope, BC_cons, V, R) + + ! Ricalcolo i flussi + + call fluxo(U_app_1_2, k, Dx_vec, sslope, BC_cons, V, R, flusso) + + d_flux = flusso(:, 2:size(flusso, 2)) - flusso(:, 1:size(flusso, 2) - 1) + + ! Secondo passo + + U_app_2 = - 0.5d0 * k * d_flux / Dx_mat + + U_app_2_2 = U_app_2 / 2.d0 + stato + + ! Calcolo delle pendenze + + call slope_fin(Dx_vec, U_app_2_2, BC_fis, sslope, BC_cons, V, R) + + ! Ricalcolo i flussi + + call fluxo(U_app_2_2, k, Dx_vec, sslope, BC_cons, V, R, flusso) + + d_flux = flusso(:, 2:size(flusso, 2)) - flusso(:, 1:size(flusso, 2) - 1) + + ! Terzo passo + + U_app_3 = - 0.5d0 * k * d_flux / Dx_mat + + U_app_3_2 = U_app_3 + stato + + ! Calcolo delle pendenze + + call slope_fin(Dx_vec, U_app_3_2, BC_fis, sslope, BC_cons, V, R) + + ! Ricalcolo i flussi + + call fluxo(U_app_3_2, k, Dx_vec, sslope, BC_cons, V, R, flusso) + + d_flux = flusso(:, 2:size(flusso, 2)) - flusso(:, 1:size(flusso, 2) - 1) + + ! Quarto passo + + U_app_4 = - 0.5d0 * k * d_flux / Dx_mat + + ! Update della soluzione + + U_next = stato + (1/6.d0) * (U_app_1 + 2 * U_app_2 + 2 * U_app_3 + U_app_4) + + endselect + + + + end function + +endmodule RK diff --git a/SW_completo.f90 b/SW_completo.f90 new file mode 100644 index 0000000..d4166d3 --- /dev/null +++ b/SW_completo.f90 @@ -0,0 +1,171 @@ +PROGRAM SW_completo + + ! Solutore per le Shallow-Water Equations in 1D + + USE gnufor_mio + + USE LOGO + + use problema_vero + + use fluxes + + ! Discretizzazione nello spazio + + INTEGER, parameter :: Np = 100 + + real(kind = 8), dimension(Np) :: x_vec + + real(kind = 8), dimension(Np - 2) :: Dx_vec + + real(kind = 8) :: Dx, t, Dt + + real(kind = 8), parameter :: t0 = 0, Tf = 10.d0, xs = 0.d0, xd = 10.d0 + + real(kind = 8) :: xm = (xd + xs)/2.d0 + + real(kind = 8), parameter :: CFL = 0.9 + + ! Problema fisico + + real(kind = 8), dimension(Np) :: h_i, zb = 0.0d0, u_i = 0.0d0, d_zb = 0.0d0, cf = 0.1 + + real(kind = 8) :: cf_wind = 0.001, W_wind = 10.d0 + + real(kind = 8), dimension(2, Np) :: w, w_old, lambda, s + + real(kind = 8), dimension(2, 2, Np) :: R, L + + real(kind = 8), dimension(2, 2) :: BC + + real(kind = 8), dimension(2, Np + 1) :: LW, UP + + real(kind = 8), dimension(Np, 1000) :: h_mat, u_mat, x_mat + + ! Variabili utili per i plot + + CHARACTER(len = 50) :: title + + CHARACTER(len = 50), dimension(:), allocatable :: legend + + ! Interi utili + + INTEGER :: j + + !---------------------------------------------------------------------------- + ! LOGO FIGO + !---------------------------------------------------------------------------- + + call disegna_logo + + !---------------------------------------------------------------------------- + ! DEFINIZIONE DEL DOMINIO DI CALCOLO + !---------------------------------------------------------------------------- + + ! Costruzione del dominio di calcolo + + Dx = abs(xd - xs) / (Np - 1) + + x_vec(1) = xs; x_vec(Np) = xd + + do j = 2, Np - 1 + + x_vec(j) = x_vec(j - 1) + Dx + + if ( j == 2 .or. j == Np - 1 ) then + + Dx_vec(j - 1) = Dx/2 + + else + + Dx_vec(j - 1) = Dx + + endif + + enddo + + !h_i = 1 + 0.02 * sin(2 * 3.14d0 * (xd - x_vec)/ (xd / 4)) ! Condizione iniziale sull'altezza sinusoidale + h_i = 1 - 0.02 * exp(-(x_vec - xm) ** 2) ! Condizione iniziale sull'altezza gaussiana + + if(allocated(legend)) deallocate(legend); + + allocate(legend(1)) + + legend(1) = 'Altezza pelo libero' + + title = 'boh' + + !call plot_function(x_vec, h_i, legend, title) + + w_old(1, :) = (h_i - zb); w_old(2, :) = (h_i - zb) * u_i + + h_mat(:, 1) = h_i + u_mat(:, 1) = u_i + x_mat(:, 1) = x_vec + + !write(*,*) "w_old(1, :) = ", w_old(1, :) + + j = 2 + + do while (t < 7.6) + + call eigenstructure(w_old, lambda, R, L) + + Dt = Dx * CFL / (max(maxval(lambda(1, :)), maxval(lambda(2, :)))) + + t = t + Dt + + write(*,*) "---------------------------------------------------------------" + write(*,*) "Istante t = ", t + write(*,*) "---------------------------------------------------------------" + + ! Bordi entrambe aperti + + BC(1, 1) = w_old(1, 1); BC(2, 1) = w_old(2, 1) / w_old(1, 1); ! Sono scritte in variabili fisiche + BC(1, 2) = w_old(1, Np); BC(2, 2) = w_old(2, Np) / w_old(1, Np); + + ! Bordi entrambi chiusi, non funziona per ora + + !BC(1, 1) = w_old(1, 1); BC(2, 1) = 0; + !BC(1, 2) = w_old(1, Np); BC(2, 2) = 0; + + call flux_LW(w_old, Dt, Dx_vec, BC, LW) + call sorgente(w_old, cf, d_zb, s) + + w = w_old - (Dt / (2 * Dx)) * (LW(:, 2:size(LW, 2)) - LW(:, 1:size(LW, 2) - 1)) -& + s * Dt + + + write(*,*) "LW(:, 1) = ", LW(:, 1) + write(*,*) "LW(:, 2) = ", LW(:, 2) + write(*,*) "j = ", j + + h_mat(:, j) = w(1, :) - zb + u_mat(:, j) = w(2, :) / w(1, :) + x_mat(:, j) = x_vec + + j = j + 1 + w_old = w + + + enddo + + if(allocated(legend)) deallocate(legend); + + allocate(legend(2)) + + legend(1) = 'altezza' + legend(2) = 'fondale' + + !call plot_2functions(x_vec, h_mat(:, 2), zb, legend, title) + + if(allocated(legend)) deallocate(legend); + + allocate(legend(1)) + + legend(1) = 'altezza' + + !call plot_video(x_mat, h_mat, title, legend) + + +endprogram SW_completo diff --git a/conversione.f90 b/conversione.f90 new file mode 100644 index 0000000..5668622 --- /dev/null +++ b/conversione.f90 @@ -0,0 +1,72 @@ +module conversione + + ! Modulo per la conversione da variabili conservative + ! a variabili fisiche e viceversa + ! w = variabili conservative = [a, au] + ! p = variabili fisiche = [a, u] + + implicit none + + interface cons2fis + + module procedure cons2fis_s, cons2fis_v + + endinterface cons2fis + + interface fis2cons + + module procedure fis2cons_s, fis2cons_v + + endinterface fis2cons + +contains + + subroutine cons2fis_s(w, p) + + real(kind = 8), dimension(:), intent(in) :: w + + real(kind = 8), dimension(size(w)), intent(out) :: p + + p(1) = w(1) + + p(2) = w(2) / w(1) + + endsubroutine cons2fis_s + + subroutine cons2fis_v(ww, pp) + + real(kind = 8), dimension(:, :), intent(in) :: ww + + real(kind = 8), dimension(size(ww, 1), size(ww, 2)), intent(out) :: pp + + pp(1, :) = ww(1, :) + + pp(2, :) = ww(2, :) / ww(1, :) + + endsubroutine cons2fis_v + + subroutine fis2cons_s(p, w) + + real(kind = 8), dimension(:), intent(in) :: p + + real(kind = 8), dimension(size(p)), intent(out) :: w + + w(1) = p(1) + + w(2) = p(2) * p(1) + + endsubroutine fis2cons_s + + subroutine fis2cons_v(pp, ww) + + real(kind = 8), dimension(:, :), intent(in) :: pp + + real(kind = 8), dimension(size(pp, 1), size(pp, 2)), intent(out) :: ww + + ww(1, :) = pp(1, :) + + ww(2, :) = pp(2, :) * pp(1, :) + + endsubroutine fis2cons_v + +endmodule conversione diff --git a/fix_BC.f90 b/fix_BC.f90 new file mode 100644 index 0000000..ce57911 --- /dev/null +++ b/fix_BC.f90 @@ -0,0 +1,121 @@ +module fix_BC + + ! modulo per fissare le BC (periodiche, stazionarie, muro, parete etc..) + ! Le condizioni al contorno vengono assegnate in variabili fisiche + + use conversione + + implicit none + +contains + + subroutine BCs(flag, stato_sx, stato_dx, BC, Bound_speed, t) + + real(kind = 8), dimension(:), intent(in) :: stato_sx, stato_dx + + integer, dimension(:), intent(in) :: flag + + logical, dimension(:), intent(in) :: Bound_speed + + real(kind = 8), dimension(:), optional :: t ! t = [Tf t_att] + + real(kind = 8), dimension(size(stato_dx), size(stato_dx)), intent(inout) :: BC + + + + select case (flag(1)) + + case (1) + + write(*,*) "A sinistra ho condizioni al contorno Wall" + + ! Si suppone che gli stato abbiano m = 0 e che non ci sia la formazione + ! di vuoto + + call cons2fis(stato_sx, BC(:, 1)) + + if(Bound_speed(1)) then + + if(present(t)) then + + BC(2, 1) = 1 - sin(0.1 * 2 * 3.14d0 * (t(1) - t(2)) / t(1)) + + else + + write(*,*) "Manca il tempo come input" + + return + + endif + + else + + BC(2, 1) = 0 + + endif + + write(*,*) "BC_s = ", BC(:, 1) + + case(2) + + write(*,*) "A sinistra ho condizioni al contorno Open" + + ! Si suppone che gli stato abbiano m = 0 e che non ci sia la formazione + ! di vuoto + + call cons2fis(BC(:, 1), BC(:, 1)) + + endselect + + !!!!! FIX IL CASO DI WALL PERCHé NON DEVE RICALCOLARLE NELLE ALTRE FUNZIONI. + !!!! LA U è FISSATA!! + + select case (flag(2)) + + case (1) + + write(*,*) "A destra ho condizioni al contorno Wall" + + ! Si suppone che gli stato abbiano m = 0 e che non ci sia la formazione + ! di vuoto + + call cons2fis(stato_dx, BC(:, 2)) + + if(Bound_speed(2)) then + + if(present(t)) then + + BC(2, 2) = 1 - sin(0.1 * 2 * 3.14d0 * (t(1) - t(2)) / t(1)) + + else + + write(*,*) "Manca il tempo come input" + + return + + endif + + else + + BC(2, 2) = 0 + + endif + + write(*,*) "BC_d = ", BC(:, 2) + + case(2) + + write(*,*) "A destra ho condizioni al contorno Open" + + ! Si suppone che gli stato abbiano m = 0 e che non ci sia la formazione + ! di vuoto + + call cons2fis(BC(:, 2), BC(:, 2)) + + endselect + + endsubroutine + + + +endmodule fix_BC diff --git a/fluxes.f90 b/fluxes.f90 new file mode 100644 index 0000000..7858418 --- /dev/null +++ b/fluxes.f90 @@ -0,0 +1,98 @@ +module fluxes + + ! Modulo per il calcolo dei flussi + + use problema_vero + + use Boundary_conditions + + use conversione + + implicit none + +contains + + subroutine flux_LW(w, Dt, Dx_Vec, BC, LW) + + real(kind = 8), dimension(:, :), intent(inout) :: w, BC + + real(kind = 8), intent(in) :: Dt + + real(kind = 8), dimension(:), intent(in) :: Dx_vec + + real(kind = 8), dimension(size(w, 1), size(w, 2) + 1), intent(out) :: LW + + real(kind = 8), dimension(size(w, 1), size(w, 1), size(w, 2) - 1) :: A + + real(kind = 8), dimension(size(w, 1), size(w, 2)) :: flusso + + real(kind = 8), dimension(size(w, 1), size(w, 2) - 1) :: w_medio, flusso_medio,& + flusso_meno_medio!, A_appoggio + + !real(kind = 8), dimension(size(w, 1)) :: A_appoggio + real(kind = 8), dimension(size(w, 1), size(w, 1)) :: A_appoggio + + real(kind = 8), dimension(size(w, 1), 2) :: BC_fis + + real(kind = 8), dimension(size(Dx_vec)) :: C + integer :: j + + C = Dt / (2 * Dx_vec) + + call flusso_vero(w, flusso) + + ! BC brutte + + w_medio = (w(:, 1:size(w, 2) - 1) + w(:, 2:size(w, 2))) / 2 + + call jacobian(w_medio, A) + + flusso_medio = (flusso(:, 1:size(w, 2) - 1) + flusso(:, 2:size(w, 2))) / 2 + flusso_meno_medio = (flusso(:, 2:size(w, 2)) - flusso(:, 1:size(w, 2) - 1)) / 2 + + do j = 1, size(A, 3) + + A_appoggio = A(:, :, j) + + LW(:, j + 1) = flusso_medio(:, j) -& + C * matmul(A_appoggio, flusso_meno_medio(:, j)) + + enddo + + !call cons2fis(BC, BC_fis) + + !write(*,*) "BC_fis = ", BC_fis + + !call BC_dx(w(:, size(w, 2)), BC_fis(:, 2)) !! ERRORE QUI, NON SO PERCHè + !call BC_sx(w(:, 1), BC_fis(:, 1)) + + !write(*,*) "BC_fis(:, 1) = ", BC_fis(:, 1) + + !call fis2cons(BC_fis, BC) + + !write(*,*) "BC = ", BC + + call cons2fis(BC, BC) + + !write(*,*) "BC_fis = ", BC + + call BC_dx(w(:, size(w, 2)), BC(:, 2)) + call BC_sx(w(:, 1), BC(:, 1)) + + !write(*,*) "BC_fis(:, 1) = ", BC(:, 1) + + call fis2cons(BC, BC) + + !write(*,*) "BC = ", BC + + call flusso_vero(BC(:, 1), LW(:, 1)) + !print *, "BC(:, 1) = ", BC(:, 1) + !print *, "BC(:, 2) = ", BC(:, 2) + call flusso_vero(BC(:, 2), LW(:, size(LW, 2))) + + + endsubroutine flux_LW + + + +endmodule fluxes diff --git a/gnufor_mio.f90 b/gnufor_mio.f90 new file mode 100644 index 0000000..0aa7da1 --- /dev/null +++ b/gnufor_mio.f90 @@ -0,0 +1,2029 @@ +MODULE gnufor_mio + + IMPLICIT NONE + + PRIVATE get_unit, run_gnuplot, timestamp, write_vector_data, & + write_vector_plot, write_xy_data, write_xy_data_video, & + write_xy_plot, write_xy_plot_video + + ! useful functions + PUBLIC plot_function, plot_video + +CONTAINS + +! LQ AG + +subroutine plot_2functions (x, y1, y2, legend, title, N) + + implicit none + + real(kind=8), dimension(:), intent(in) :: x, y1, y2 + + real(kind=8), dimension(size(x), 3) :: xx + + character(len = *), dimension(:), intent(in) :: legend + + integer, dimension(:,:), optional, intent(inout) :: N + + character(len = *), intent(in) :: title + + character ( len = 100 ) :: command_file_name = 'plot_function_commands.txt' + character ( len = 100 ) :: data_file_name = 'plot_function_data.txt' + + integer :: ierror + + xx(:,1) = x + xx(:,2) = y1 + xx(:,3) = y2 + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a)' ) ' To plot a simple set of (X,Y) data,' + write ( *, '(a)' ) ' WRITE_XY_DATA writes the data file,' + write ( *, '(a)' ) ' WRITE_XY_PLOT writes the plot command file.' + + + call write_xyy_data ( data_file_name, size(x), size(x), 3, xx, ierror) + ! write_xyy_data ( data_file_name, lda, nrow, ncol, x, ierror) + + + if ( ierror /= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a,i6)' ) ' WRITE_XY_DATA returned IERROR = ', ierror + end if + +if(present(N)) then + + call write_xyy_plots ( command_file_name, data_file_name, 3, legend, title, ierror, N ) + ! write_xyy_plots ( command_file_name, data_file_name, ncol, legend, title, ierror, N ) + +else + + call write_xyy_plots ( command_file_name, data_file_name, 3, legend, title, ierror) + ! write_xyy_plots ( command_file_name, data_file_name, ncol, legend, title, ierror) + +endif + + if ( ierror /= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a,i6)' ) ' WRITE_XY_PLOT returned IERROR = ', ierror + end if + + call run_gnuplot ( command_file_name ) + +end subroutine plot_2functions + +subroutine plot_3functions (x, y1, y2, y3, legend, title, N) + + implicit none + + real(kind=8), dimension(:), intent(in) :: x, y1, y2, y3 + + real(kind=8), dimension(size(x), 4) :: xx + + character( len = *), dimension(:), intent(in) :: legend + + integer, dimension(:,:), optional, intent(inout) :: N + + character(len = *), intent(in) :: title + + character ( len = 100 ) :: command_file_name = 'plot_function_commands.txt' + character ( len = 100 ) :: data_file_name = 'plot_function_data.txt' + + integer :: ierror + + xx(:,1) = x + xx(:,2) = y1 + xx(:,3) = y2 + xx(:,4) = y3 + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a)' ) ' To plot a simple set of (X,Y) data,' + write ( *, '(a)' ) ' WRITE_XY_DATA writes the data file,' + write ( *, '(a)' ) ' WRITE_XY_PLOT writes the plot command file.' + + + call write_xyy_data ( data_file_name, size(x), size(x), 4, xx, ierror) + ! write_xyy_data ( data_file_name, lda, nrow, ncol, x, ierror) + + if ( ierror /= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a,i6)' ) ' WRITE_XY_DATA returned IERROR = ', ierror + end if + +if(present(N)) then + + call write_xyy_plots ( command_file_name, data_file_name, 4, legend, title, ierror, N ) + ! write_xyy_plots ( command_file_name, data_file_name, ncol, legend, title, ierror, N ) + +else + + call write_xyy_plots ( command_file_name, data_file_name, 4, legend, title, ierror) + ! write_xyy_plots ( command_file_name, data_file_name, ncol, legend, title, ierror) + +endif + if ( ierror /= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a,i6)' ) ' WRITE_XY_PLOT returned IERROR = ', ierror + end if + + call run_gnuplot ( command_file_name ) + +end subroutine plot_3functions + +subroutine plot_4functions (x, y1, y2, y3, y4, legend, title, N) + + implicit none + + real(kind=8), dimension(:), intent(in) :: x, y1, y2, y3, y4 + + real(kind=8), dimension(size(x), 5) :: xx + + character( len = *), dimension(:), intent(in) :: legend + + integer, dimension(:,:), optional, intent(inout) :: N + + character(len = *), intent(in) :: title + + character ( len = 100 ) :: command_file_name = 'plot_function_commands.txt' + character ( len = 100 ) :: data_file_name = 'plot_function_data.txt' + + integer :: ierror + + xx(:,1) = x + xx(:,2) = y1 + xx(:,3) = y2 + xx(:,4) = y3 + xx(:,5) = y4 + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a)' ) ' To plot a simple set of (X,Y) data,' + write ( *, '(a)' ) ' WRITE_XY_DATA writes the data file,' + write ( *, '(a)' ) ' WRITE_XY_PLOT writes the plot command file.' + + + call write_xyy_data ( data_file_name, size(x), size(x), 5, xx, ierror) + ! write_xyy_data ( data_file_name, lda, nrow, ncol, x, ierror) + + + if ( ierror /= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a,i6)' ) ' WRITE_XY_DATA returned IERROR = ', ierror + end if + +if(present(N)) then + + call write_xyy_plots ( command_file_name, data_file_name, 5, legend, title, ierror, N ) + ! write_xyy_plots ( command_file_name, data_file_name, ncol, legend, title, ierror, N ) + +else + + call write_xyy_plots ( command_file_name, data_file_name, 5, legend, title, ierror) + ! write_xyy_plots ( command_file_name, data_file_name, ncol, legend, title, ierror) + + +endif + if ( ierror /= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a,i6)' ) ' WRITE_XY_PLOT returned IERROR = ', ierror + end if + + call run_gnuplot ( command_file_name ) + +end subroutine plot_4functions + +subroutine plot_muchmorethan2functions (x, yy, N, legend, title, M) + + ! x = vettore delle x ( Capitan Ovvio) + ! yy = vettore lungo N * size(x) contenente incolonnate una sotto l'altra + ! le y delle funzioni da plottare + ! N = numero di funzioni da plottare + + implicit none + + real(kind=8), dimension(:), intent(in) :: x + + real(kind = 8), dimension(:), intent(in) :: yy + + character( len = *), dimension(:), intent(in) :: legend + + integer, dimension(:,:), optional, intent(inout) :: M + + character(len = *), intent(in) :: title + + integer, intent(in) :: N + + real(kind=8), dimension(size(x), N + 1) :: xx + + integer :: j, k, l + + + character ( len = 100 ) :: command_file_name = 'plot_function_commands.txt' + character ( len = 100 ) :: data_file_name = 'plot_function_data.txt' + + integer :: ierror + + xx(:, 1) = x + + + do j = 2, (N + 1) + + xx(:, j) = yy( 1 + (j - 2) * size(x) : (j - 1) * size(x) ) + + enddo + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a)' ) ' To plot a simple set of (X,Y) data,' + write ( *, '(a)' ) ' WRITE_XY_DATA writes the data file,' + write ( *, '(a)' ) ' WRITE_XY_PLOT writes the plot command file.' + + + call write_xyy_data ( data_file_name, size(x), size(x), N + 1, xx, ierror) + ! write_xyy_data ( data_file_name, lda, nrow, ncol, x, ierror) + + + if ( ierror /= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a,i6)' ) ' WRITE_XY_DATA returned IERROR = ', ierror + end if + +if(present(M)) then + + call write_xyy_plots ( command_file_name, data_file_name, N + 1, legend, title, ierror, M ) + ! write_xyy_plots ( command_file_name, data_file_name, ncol, legend, title, ierror, M ) + +else + + call write_xyy_plots ( command_file_name, data_file_name, N + 1, legend, title, ierror) + ! write_xyy_plots ( command_file_name, data_file_name, ncol, legend, title, ierror) + +endif + if ( ierror /= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a,i6)' ) ' WRITE_XY_PLOT returned IERROR = ', ierror + end if + + call run_gnuplot ( command_file_name ) + +end subroutine plot_muchmorethan2functions + + +subroutine plot_function (x, y, legend, title) + + implicit none + + real(kind=8), dimension(:), intent(in) :: x, y + + character( len = *), dimension(:), intent(in) :: legend + + character(len = *), intent(in) :: title + + character ( len = 100 ) :: command_file_name = 'plot_function_commands.txt' + character ( len = 100 ) :: data_file_name = 'plot_function_data.txt' + + integer :: ierror + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a)' ) ' To plot a simple set of (X,Y) data,' + write ( *, '(a)' ) ' WRITE_XY_DATA writes the data file,' + write ( *, '(a)' ) ' WRITE_XY_PLOT writes the plot command file.' + + call write_xy_data ( data_file_name, size(x), x, y, ierror) + + if ( ierror /= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a,i6)' ) ' WRITE_XY_DATA returned IERROR = ', ierror + end if + + call write_xy_plot ( command_file_name, data_file_name, legend, title, ierror ) + + if ( ierror /= 0 ) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_function' + write ( *, '(a,i6)' ) ' WRITE_XY_PLOT returned IERROR = ', ierror + end if + + call run_gnuplot ( command_file_name ) + +end subroutine plot_function +!******************************************************************************* +! RA, JB Riccardo Albi e Jacopo Banchetti 7 gennaio 2015 + +subroutine plot_video (XX, YY, title, legend, fn, t, save_flag ) + + ! INPUT: + ! X, matrice (eventualmente di una sola colonna) + ! Y, matrice + ! + ! fn, radice del filename, opzionale + ! t, tempo di pausa, opzionale + ! save_flag, tipo logico, se vero salva un file contenente + ! il video senza pero' mostrare subito l'animazione. + ! + ! Mostra un video, dove ogni fotogramma e' un grafico 2D + ! composto da una colonna della matrice XX e dalla corrispondente + ! colonna della matrice YY. + ! + ! Se XX contiene UNA sola colonna, i grafici saranno disegnati + ! utilizzando la stessa colonna XX(:,1) per tutte le colonne + ! della matrice YY. + ! + ! La funzione salva in ogni caso i dati relativi ai grafici + ! per un successivo eventuale uso. + ! Il parametro opzionale t permette di impostare la pausa + ! tra un fotogramma e il successivo. + ! Il parametro save_flag permette di salvare un video ottenuto + ! dall'unione dei vari grafici in formato .avi. + ! Ogni fotogramma viene salvato in formato .png prima di essere + ! assemblato con l'aiuto del programma 'ffmpeg'. + ! + ! NB: 'ffmpeg' deve essere installato per salvare il video. + + implicit none + + real(kind=8), intent(in), dimension(:,:) :: XX ! puo' essere (:,1) + real(kind=8), intent(in), dimension(:,:) :: YY + character(len=*), intent(in), optional :: fn + character(len = *), dimension(:), intent(in) :: legend + character(len = *), intent(in) :: title + real, intent(in), optional :: t + logical, intent(in), optional :: save_flag + + + character(len=64) :: fn_ + real :: t_ + + character (len = 100) :: command_file_name = 'video/plot_function_commands.txt' + character (len = 100) :: data_file_name = 'video/plot_function_data.txt' + + integer :: ierror, n_active_columns, n + + ! Se trova una colonna di 0 sulle YY, non plotta da quella colonna in poi + ! LQ + + n_active_columns = SIZE(YY, 2) + + do n = 1, SIZE(YY, 2) + + if (COUNT(YY(:, n) == 0) == SIZE(YY,1)) then + n_active_columns = n - 1 + exit + end if + + end do +! RA +! do n_active_columns = 0, SIZE(YY, 2) +! if (COUNT(YY(:, n_active_columns + 1) == 0) == SIZE(YY,1)) exit +! end do + + + if (PRESENT(fn)) THEN + fn_ = fn + else + fn_ = "__ddd" ! nome file di default + end if + + call write_xy_data_video (size(YY,1), n_active_columns, XX, YY, ierror, fn_) + + if (ierror /= 0) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_video' + write ( *, '(a,i6)' ) ' WRITE_XY_DATA_VIDEO returned IERROR = ', ierror + end if + + if (PRESENT(t)) THEN + if (t < 0 .OR. t > 10) THEN + write (*, *) '[gnufor.f90] WARNING: pause time must be in the 0-10 interval' + t_ = 0.1 + else + t_= t + end if + else + t_ = 0.1 ! tempo di default tra i fotogrammi + end if + + call write_xy_plot_video (XX, YY, command_file_name, data_file_name, & + legend, title, ierror, n_active_columns, fn_,& + t_, save_flag) + + if (ierror /= 0) then + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'plot_video' + write ( *, '(a,i6)' ) ' WRITE_XY_PLOT_VIDEO returned IERROR = ', ierror + end if + + call run_gnuplot ( command_file_name ) + + if (PRESENT(save_flag) .AND. save_flag) then + + write(*, *) 'ffmpeg -i video/'// fn_ // '%05d.png' // fn_ // '.avi' + call system ('ffmpeg -i video/'// TRIM(fn_) // '%05d.png ' // TRIM(fn_) // '.avi') + + end if + +end subroutine plot_video + +!******************************************************************************* +! +subroutine get_unit ( iunit ) +! +!******************************************************************************* +! +!! GET_UNIT returns a free FORTRAN unit number. +! +! +! Discussion: +! +! A "free" FORTRAN unit number is an integer between 1 and 99 which +! is not currently associated with an I/O device. A free FORTRAN unit +! number is needed in order to open a file with the OPEN command. +! +! Modified: +! +! 02 March 1999 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Output, integer IUNIT. +! +! If IUNIT = 0, then no free FORTRAN unit could be found, although +! all 99 units were checked (except for units 5 and 6). +! +! Otherwise, IUNIT is an integer between 1 and 99, representing a +! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 +! are special, and will never return those values. +! + implicit none +! + integer i + integer ios + integer iunit + logical lopen + + iunit = 0 + + do i = 1, 99 + + if ( i /= 5 .and. i /= 6 ) then + + inquire ( unit = i, opened = lopen, iostat = ios ) + + if ( ios == 0 ) then + if ( .not. lopen ) then + iunit = i + return + end if + end if + + end if + + end do + +end subroutine get_unit + +function pi ( ) +! +!******************************************************************************* +! +!! PI returns the value of pi. +! +! +! Modified: +! +! 04 December 1998 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Output, real PI, the value of pi. +! + implicit none +! + real(kind=8) :: pi +! + pi = 3.14159265358979323846264338327950288419716939937510E+00 + +end function pi + +subroutine run_gnuplot ( command_file_name ) + +! LQ AG +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +! NECESSARY ONLY FOR NAGWARE FORTRAN +! +! Use added for system call using function SYSTEM +! ====== + !!!!!!!!!!!!!!!1 USE f90_unix_proc + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! LQ AG + +! +!******************************************************************************* +! +!! RUN_GNUPLOT runs GNUPLOT with a given command file. +! +! +! Discussion: +! +! The GNUPLOT program, version 3.7, must be available. To check whether +! this is so, try typing +! +! which gnuplot +! +! If the response is +! +! gnuplot: command not found +! +! then you're going to have to make GNUPLOT available. +! +! At ISU, this may require that you issue the command +! +! add gnu +! +! You may need to set the environment variable GNUTERM: +! +! setenv GNUTERM x11 +! +! so that GNUPLOT automatically displays to your X window terminal. +! +! +! This routine expects that there is a text file containing the appropriate +! commands to GNUPLOT to display your picture. There are a number of +! routines in this package that will do this for simple plotting tasks. +! Most of them require that you also set up a file of data to be plotted. +! +! Once this routine invokes GNUPLOT, a graphics window should open +! up, and the FORTRAN program will pause. Hitting RETURN should advance +! to the next picture, or terminate the window at the end, allowing the +! FORTRAN routine to proceed. +! +! +! You can look at the data and command files created by the routines. +! Moreover, you can easily modify the command file to change the options +! used in GNUPLOT, and then run GNUPLOT interactively, as in: +! +! gnuplot commands +! +! In particular, if you want a PostScript version of your graphics files, +! insert the command "set term postscript" at the beginning of the command +! file and run gnuplot as follows: +! +! gnuplot commands > mypicture.ps +! +! You will also have to hit RETURN once for each plot that is made. +! +! Modified: +! +! 21 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) COMMAND_FILE_NAME, the name of the +! command file. +! + implicit none +! + character ( len = 100 ) command + character ( len = * ) command_file_name + +! LQ AG integer status +! LQ AG integer system + +! +! Issue a command to the system that will startup GNUPLOT, using +! the file we just wrote as input. +! + write ( command, * ) 'gnuplot ' // trim ( command_file_name ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'Issuing the command:' // trim ( command ) + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) ' Press RETURN to proceed.' + +! LQ AG status = system ( trim ( command ) ) + + call system ( trim ( command ) ) + +! LQ AG +! +! if ( status == 0 ) then +! write ( *, '(a)' ) ' ' +! write ( *, '(a)' ) 'RUN_GNUPLOT:' +! write ( *, '(a)' ) ' Normal end of execution of the GNUPLOT program.' +! else +! write ( *, '(a)' ) ' ' +! write ( *, '(a)' ) 'RUN_GNUPLOT - Fatal error!' +! write ( *, '(a)' ) ' An error code was returned when the GNUPLOT command' +! write ( *, '(a)' ) ' was issued. Perhaps GNUPLOT is not in your path.' +! write ( *, '(a)' ) ' Type "which gnuplot" to check this.' +! stop +! end if +! LQ AG + +end subroutine run_gnuplot + +subroutine timestamp ( ) +! +!******************************************************************************* +! +!! TIMESTAMP prints the current YMDHMS date as a time stamp. +! +! +! Example: +! +! May 31 2001 9:45:54.872 AM +! +! Modified: +! +! 31 May 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! None +! + implicit none +! + character ( len = 8 ) ampm + integer d + character ( len = 8 ) date + integer h + integer m + integer mm + character ( len = 9 ), parameter, dimension(12) :: month = (/ & + 'January ', 'February ', 'March ', 'April ', & + 'May ', 'June ', 'July ', 'August ', & + 'September', 'October ', 'November ', 'December ' /) + integer n + integer s + character ( len = 10 ) time + integer values(8) + integer y + character ( len = 5 ) zone +! + call date_and_time ( date, time, zone, values ) + + y = values(1) + m = values(2) + d = values(3) + h = values(5) + n = values(6) + s = values(7) + mm = values(8) + + if ( h < 12 ) then + ampm = 'AM' + else if ( h == 12 ) then + if ( n == 0 .and. s == 0 ) then + ampm = 'Noon' + else + ampm = 'PM' + end if + else + h = h - 12 + if ( h < 12 ) then + ampm = 'PM' + else if ( h == 12 ) then + if ( n == 0 .and. s == 0 ) then + ampm = 'Midnight' + else + ampm = 'AM' + end if + end if + end if + + write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & + trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) + +end subroutine timestamp + +subroutine write_vector_data ( data_file_name, n, x, y, dx, dy, ierror ) +! +!******************************************************************************* +! +!! WRITE_VECTOR_DATA writes vector data to a file, for plotting by GNUPLOT. +! +! +! Discussion: +! +! Each vector is described by 4 values, X, Y, dX, dY, indicating that +! a vector is to be drawn from (X,Y) to (X+dX,Y+dY). +! +! Modified: +! +! 22 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Input, integer N, the number of vectors. +! +! Input, real X(N), Y(N), DX(N), DY(N), the vector data. +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + integer n +! + character ( len = * ) data_file_name + real(kind=8) :: dx(n) + real(kind=8) :: dy(n) + integer file_unit + integer i + integer ierror + integer ios + real(kind=8) :: x(n) + real(kind=8) :: y(n) +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_VECTOR_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = data_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_VECTOR_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + do i = 1, n + write ( file_unit, * ) x(i), y(i), dx(i), dy(i) + end do + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_VECTOR_DATA:' + write ( *, '(a)' ) ' Wrote the GNUPLOT vector data file "' // & + trim ( data_file_name ) // '"' + +end subroutine write_vector_data + +subroutine write_vector_plot ( command_file_name, data_file_name, & + ierror ) +! +!******************************************************************************* +! +!! WRITE_VECTOR_PLOT writes GNUPLOT commands to plot vectors. +! +! +! Modified: +! +! 22 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) COMMAND_FILE_NAME, the name of the +! command file. +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + character ( len = * ) command_file_name + character ( len = * ) data_file_name + integer file_unit + integer ierror + integer ios +! +! Write the data file. +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_VECTOR_PLOT - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = command_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_VECTOR_PLOT - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + write ( file_unit, '(a)' ) 'set title "GNUFOR plot"' + write ( file_unit, '(a)' ) 'set xlabel "x"' + write ( file_unit, '(a)' ) 'set ylabel "y"' + write ( file_unit, '(a)' ) 'set data style vector' + write ( file_unit, '(a,i2,a)' ) 'plot "' // trim ( data_file_name ) + write ( file_unit, '(a)' ) 'pause -1' + write ( file_unit, '(a)' ) 'q' + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_VECTOR_PLOT:' + write ( *, '(a)' ) ' Wrote the GNUPLOT table plots command file "' // & + trim ( command_file_name ) // '"' + +end subroutine write_vector_plot + +subroutine write_xy_data ( data_file_name, n, x, y, ierror ) +! +!******************************************************************************* +! +!! WRITE_XY_DATA writes X(1:N), Y(1:N) data to a file. +! +! +! Modified: +! +! 23 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Input, integer N, the number of data items. +! +! Input, real X(N), Y(N), the X and Y data +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + integer n +! + character ( len = * ) data_file_name + integer file_unit + integer i + integer ierror + integer ios + real(kind=8) :: x(n) + real(kind=8) :: y(n) +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = data_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + do i = 1, n + write ( file_unit, * ) x(i), y(i) + end do + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_DATA:' + write ( *, '(a)' ) ' Wrote the GNUPLOT XY data file "' // & + trim ( data_file_name ) // '"' + +end subroutine write_xy_data + +! RA, JB + +subroutine write_xy_data_video ( nr, nc, x, y, ierror , data_file_name_base) + + implicit none + + integer, intent(in) :: nr, nc + real(kind=8), dimension(:, :), intent(in) :: x + real(kind=8), dimension(:, :), intent(in) :: y + real(kind=8), dimension(size(x, 1), size(y, 1)/size(x, 1) + 1) :: xx + integer, intent(out) :: ierror + + character ( len = * ), intent(in) :: data_file_name_base + + character ( len = 128 ) data_file_name_current + character ( len = * ), parameter :: base_dir = "video/" + integer file_unit + integer i, j, k, s + integer ios + + ierror = 0 + + write(*,*) "Creating directory for frames storage: "//base_dir + call system ( "mkdir -p "//base_dir ) + + call get_unit ( file_unit ) + + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + write(*,*) "size(y, 1) = ", size(y, 1) + write(*,*) "size(x, 1) = ", size(x, 1) + + if(SIZE(X, 2) == SIZE(Y, 2)) then + + ! many X versus many Y + + do j = 1, nc + + WRITE(data_file_name_current, '(a, a, i5.5)') TRIM(base_dir), TRIM(data_file_name_base), j + + open ( unit = file_unit, file = data_file_name_current, status = 'replace', & + iostat = ios ) + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + +!AR + + xx(:, 1) = x(:, j) + + do s = 1, size(y, 1) / size(x, 1) + + xx(:, s + 1) = y( 1 + (s - 1) * size(x, 1) : s * size(x, 1), j) + + enddo + + !write ( file_unit, * ) xx + +! LQ + do i = 1, size(x, 1) + !k = ((i-1)/size(x(:,j))) * size(x(:,j)) ! stile anti Auteri + write ( file_unit, * ) xx(i, :) + ! write ( file_unit, * ) x(i-((i-1)/size(x(:,j))) * size(x(:,j)), 1), y(i, j) + end do + + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_DATA_VIDEO:' + write ( *, '(a)' ) ' writing "' // trim ( data_file_name_current ) // '"' + + end do + + else if (size(X, 2) == 1) then + + ! print many Y versus one X only + + do j = 1, nc + + WRITE(data_file_name_current, '(a, a, i5.5)') TRIM(base_dir), TRIM(data_file_name_base), j + + open ( unit = file_unit, file = data_file_name_current, status = 'replace', & + iostat = ios ) + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_DATA_VIDEO - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + +!AR + + xx(:, 1) = x(:, 1) + + do s = 1, size(y, 1) / size(x, 1) + + xx(:, s + 1) = y( 1 + (s - 1) * size(x, 1) : s * size(x, 1), j) + + enddo + + !write ( file_unit, * ) xx + +! LQ + do i = 1, size(x, 1) + !k = ((i-1)/size(x(:,1))) * size(x(:,1)) ! stile anti Auteri + write ( file_unit, * ) xx(i, :) + end do + + close ( unit = file_unit ) +! +! write ( *, '(a)' ) ' ' +! write ( *, '(a)' ) 'WRITE_XY_DATA_VIDEO:' +! write ( *, '(a)' ) ' writing "' // trim ( data_file_name_current ) // '"' +! + end do + + else + + ierror = 3 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_DATA_VIDEO - Fatal error!' + write ( *, '(a)' ) ' dimension mismatch' + return + + end if + +end subroutine write_xy_data_video + +subroutine write_xy_plot ( command_file_name, data_file_name, legend, title, ierror ) +! +!******************************************************************************* +! +!! WRITE_XY_PLOT writes GNUPLOT commands to plot X(1:N), Y(1:N) data. +! +! +! Modified: +! +! 23 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) COMMAND_FILE_NAME, the name of the +! command file. +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + character ( len = * ) command_file_name + character ( len = * ) data_file_name + character ( len = * ), dimension(:) :: legend + character ( len = * ) title + integer file_unit +! LQ integer i + integer ierror + integer ios +! +! Write the data file. +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_PLOT - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = command_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_PLOT - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + write ( file_unit, '(a)' ) 'set title "'// trim(title) // '"' + write ( file_unit, '(a)' ) 'set xlabel "x"' + write ( file_unit, '(a)' ) 'set ylabel "y"' + write ( file_unit, '(a, i2, a)' ) 'plot "' // trim ( data_file_name ) // & + '" using 1:2 title "' // trim(legend(1)) // '" with lines lc 1' + write ( file_unit, '(a)' ) 'pause mouse' + write ( file_unit, '(a)' ) 'q' + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_PLOT:' + write ( *, '(a)' ) ' Wrote the GNUPLOT XY plot command file "' // & + trim ( command_file_name ) // '"' + +end subroutine write_xy_plot + +! RA, JB + +subroutine write_xy_plot_video (x, y, command_file_name, data_file_name, legend, title, ierror, nc, fn, t, save_flag) + + implicit none + + real(kind=8), dimension(:, :) :: x, y + character ( len = * ) command_file_name + character ( len = * ) data_file_name + character ( len = * ) title + character(len = *), dimension(:) :: legend + integer ierror + + character ( len = * ) fn + + integer file_unit + integer ios + integer nc + integer j + + real a1, b1, a2, b2, c1, c2 + + integer int_inf, int_sup, o + logical, optional :: save_flag + + real t + + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_PLOT_VIDEO - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = command_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XY_PLOT_VIDEO - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + write ( file_unit, '(a)' ) 'cd "video/"' + write ( file_unit, '(a)' ) 'set title "' // trim(title) // '"' + write ( file_unit, '(a)' ) 'set xlabel "x"' + write ( file_unit, '(a)' ) 'set ylabel "y"' + + if(size(x, 2) == 1) then + a1 = minval(x(:, 1)) * (1 - sign(0.05d0, minval(x(:, 1)))) + b1 = maxval(x(:, 1)) * (1 + sign(0.05d0, maxval(x(:, 1)))) + write ( file_unit, '(a, f10.2, a, f10.2,a)' ) 'set xrange [', a1, ':', b1 , ']' + else + a2 = minval(x(:, 1:nc)) * (1 - sign(0.05d0, minval(x(:, 1:nc)))) + b2 = maxval(x(:, 1:nc)) * (1 + sign(0.05d0, maxval(x(:, 1:nc)))) + write ( file_unit, '(a, f10.2, a, f10.2,a)' ) 'set xrange [', a2,':', b2,']' + end if + + c1 = minval(y(:, 1:nc)) * (1 - sign(0.05d0,minval(y(:, 1:nc)))) + c2 = maxval(y(:, 1:nc)) * (1 + sign(0.05d0, maxval(y(:, 1:nc)))) + + write ( file_unit, '(a, f10.2, a, f10.2,a)' ) 'set yrange [', c1,':',c2,']' + + + if (PRESENT(save_flag)) THEN + if (save_flag) THEN + write ( file_unit, '(a)' ) 'set terminal png size 1920,1080' ! 'set term png' + end if + end if + + do j=1, nc + + write (data_file_name, '(a, i5.5)') TRIM(fn), j + + + if(PRESENT(save_flag) .AND. save_flag) THEN + write ( file_unit, '(a)' ) 'set output ''' // trim ( data_file_name ) // '.png''' + end if + + if (size(x, 1) /= size(y, 1)) then + + write ( file_unit, '(a)' ) 'set multiplot ' + + do o = 1, size(y, 1)/size(x, 1) + + int_inf = (o - 1)*size(x, 1) + int_sup = o*size(x, 1) - 1 + !int_sup = int_inf + 1 + !write ( file_unit, '(a, i4, a, i4, a, f6.5, a, i1, a)' ) 'plot "' // trim ( data_file_name ) // & + ! '" every ::',int_inf, '::',int_sup,' using 1:2 title "' // trim( legend(o) ) // & + ! '" offset 0,',-(o-1)/5.d0,'with lines lc ',(o+1),' lw 2' + + if( o == 1) then + + write ( file_unit, '(a, i4, a, i4, a, i1, a)' ) 'plot "' // trim ( data_file_name ) // & + '" using 1:',(o+1),' title "' // trim(legend(o)) // & + '" with lines lc ',(o+1),' lw 2,\' + + elseif( o /= size(y, 1)/size(x, 1) .and. o /= 1) then + + write ( file_unit, '(a, i4, a, i4, a, i1, a)' ) '"' // trim ( data_file_name ) // & + '" using 1:',(o+1),' title "' // trim(legend(o)) // & + '" with lines lc ',(o+1),' lw 2,\' + + else + + write ( file_unit, '(a, i4, a, i4, a, i1, a)' ) '"' // trim ( data_file_name ) // & + '" using 1:',(o+1),' title "' // trim(legend(o)) // & + '" with lines lc ',(o+1),' lw 2' + + endif + + end do + + if (PRESENT(save_flag) .AND. save_flag) THEN + write ( file_unit, '(a)' ) 'unset multiplot ' + end if + + else + write ( file_unit, '(a, i2, a)' ) 'plot "' // trim ( data_file_name ) // & + '" using 1:2 title "' // trim( legend(1) ) // & + '" with lines lc 1' + end if + + if(.NOT. PRESENT(save_flag) .OR. (.NOT. save_flag)) THEN + write ( file_unit, '(a, F4.2)' ) 'pause ', t + end if + + end do + + if(.NOT.PRESENT(save_flag).OR. (.NOT. save_flag)) THEN + write ( file_unit, '(a)' ) 'pause -1' + end if + + close ( unit = file_unit ) + +end subroutine write_xy_plot_video + +subroutine write_xyy_data ( data_file_name, lda, nrow, ncol, x, ierror) +! +!******************************************************************************* +! +!! WRITE_XYY_DATA writes a table of data to a file, for plotting by GNUPLOT. +! +! +! Discussion: +! +! The first column of data is assumed to be the independent variable, X. +! Separate plots are made of X versus all the other columns of data. +! +! Modified: +! +! 21 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Input, integer LDA, the leading dimension of X. +! +! Input, integer NROW, NCOL, the dimensions of X. +! +! Input, real X(LDA,NCOL), the NROW by NCOL data to be written. +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + integer lda + integer ncol +! + character ( len = * ) data_file_name + integer file_unit + integer i + integer ierror + integer ios + + integer nrow + real(kind=8) :: x(lda,ncol) +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYY_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = data_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYY_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + +do i = 1, nrow + write ( file_unit, * ) x(i,1:ncol) +end do + +close ( unit = file_unit ) + + + + + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYY_DATA:' + write ( *, '(a)' ) ' Wrote the GNUPLOT XYY data file "' // & + trim ( data_file_name ) // '"' + +end subroutine write_xyy_data + +subroutine write_xyy_plots ( command_file_name, data_file_name, ncol, legend, title, ierror, N ) +! +!******************************************************************************* +! +!! WRITE_XYY_PLOTS writes GNUPLOT commands to make multiple (X,Y) plots. +! +! +! Discussion: +! +! The first column of data is assumed to be the independent variable, X. +! Separate plots are made of X versus all the other columns of data. +! +! Modified: +! +! 23 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) COMMAND_FILE_NAME, the name of the +! command file. +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Input, integer NCOL, the number of columns of data. +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + character ( len = * ) command_file_name + character ( len = * ) data_file_name + character ( len = * ), dimension(:) :: legend + character ( len = * ) title + character ( len = 10) :: start + character ( len = 10) :: end + integer, dimension(:,:), optional, intent(inout) :: N + integer file_unit + integer i + integer ierror + integer ios + integer ncol +! +! Write the data file. +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYY_PLOTS - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = command_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYY_PLOTS - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + write ( file_unit, '(a)' ) 'set title "'// trim(title) // '"' + write ( file_unit, '(a)' ) 'set xlabel "x"' + write ( file_unit, '(a)' ) 'set ylabel "y"' + !write ( file_unit, '(a)' ) 'set label "stato sinistro" at "50, 0" offset 0, 0.5 ' + !write ( file_unit, '(a)' ) 'set label "stato destro" at "30, 0" offset 0, 0.5 ' + + + ! modifications by Jacopo Banchetti 26 settembre 2018 + !==================================================== + do i = 2, ncol + + if (i == 2) then + + write(*,*) "Sono in i == 2" + + write ( file_unit, '(a,i2,a)' ) 'plot "' // trim ( data_file_name ) // & + '" using 1:', i,' title "' // trim(legend(1)) // '" with lines ,\' + ! to split line + else + + if (i /= ncol) then + + if (.not.present(N)) then + + write(*,*) "Sono in i /= ncol e .not.present(N)" + + write ( file_unit, '(a,i2,a)' ) '"' //& + trim ( data_file_name ) // & + '" using 1:', i,'title "' // trim(legend( i - 1 )) // '" with lines ,\' + + else + + if( i <= ncol - N(1, 1) ) then ! Sono prima delle funzioni non plottate del tutto + + write(*,*) "Sono in i /= ncol, present(N) e prima delle funzioni plottate in parte" + + write ( file_unit, '(a,i2,a)' ) '"' // trim ( data_file_name ) // & + '" using 1:', i,'title "' //& + trim(legend( i - 1 )) // '" with lines ,\' + + else ! Sono nelle funzioni da plottare solo in parte + + write(*,*) "Sono in i /= ncol, present(N) e nelle delle funzioni plottate in parte" + + write(start, "(I0)") N(i - N(1,1), 1) + write(end, "(I0)") N(i - N(1,1), 2) + + write ( file_unit, '(a,i2,a)' ) '"' // trim ( 'plot_function_data.txt' ) // & + '" every :: ' // trim(start) // '::'// trim(end) //& + ' using 1:', i,'title "' //& + trim(legend( i - 1 )) // '" with lines ,\' + + + endif + + endif + + else + + if(.not.present(N)) then + + write(*,*) "Sono in i == ncol, .not.present(N) " + + ! only when i == ncol + write ( file_unit, '(a,i2,a)' ) '"' // trim ( data_file_name ) // & + '" using 1:', i,'title "' // trim(legend(size(legend))) // '" with lines' + + else + + write(*,*) "Sono in i == ncol, present(N) " + + write(*,*) "i = ", i + write(*,*) "N(1, 1) = ", N(1, 1) + write(*,*) "N(i - N(1,1), 1) = ", N(i - N(1,1), 1) + + write(start, "(I0)") N(i - N(1,1), 1) + write(end, "(I0)") N(i - N(1,1), 2) + + write(*,*) "start = ", start + + write ( file_unit, '(a,i2,a)' ) '"' // trim ( 'plot_function_data.txt' ) // & + '" every :: ' // trim(start) // '::'// trim(end) //& + ' using 1:', i,'title "' //& + trim(legend(size(legend))) // '" with lines' + + + endif + + endif + + endif + + ! JB LQ write ( file_unit, '(a)' ) 'pause -1' + + end do + + ! JB LQ + write ( file_unit, '(a)' ) 'pause -1' + ! JB LQ + + ! end of modifications by Jacopo Banchetti 26 settembre 2018 + !=========================================================== + + write ( file_unit, '(a)' ) 'q' + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYY_PLOTS:' + write ( *, '(a)' ) ' Wrote the GNUPLOT XYY plots command file "' // & + trim ( command_file_name ) // '"' + +end subroutine write_xyy_plots + +subroutine write_xyz_data ( data_file_name, n, x, y, z, ierror ) +! +!******************************************************************************* +! +!! WRITE_XYZ_DATA writes X(1:N), Y(1:N), Z(1:N) data to a file. +! +! +! Modified: +! +! 23 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Input, integer N, the number of data items. +! +! Input, real X(N), Y(N), Z(N), the X, Y, Z data +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + integer n +! + character ( len = * ) data_file_name + integer file_unit + integer i + integer ierror + integer ios + real(kind=8) :: x(n) + real(kind=8) :: y(n) + real(kind=8) :: z(n) +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZ_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = data_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZ_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + do i = 1, n + write ( file_unit, * ) x(i), y(i), z(i) + end do + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZ_DATA:' + write ( *, '(a)' ) ' Wrote the GNUPLOT XYZ data file "' // & + trim ( data_file_name ) // '"' + +end subroutine write_xyz_data + +subroutine write_xyz_plot ( command_file_name, data_file_name, ierror ) +! +!******************************************************************************* +! +!! WRITE_XYZ_PLOT writes commands to plot parametric (X,Y,Z) data. +! +! +! Discussion: +! +! This routine tries to write a command file suitable for displaying +! a 3D arc specified by points (X,Y,Z). A grid data file, containing +! values of X, Y and Z, will also be needed. +! +! Modified: +! +! 22 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) COMMAND_FILE_NAME, the name of the +! command file. +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + character ( len = * ) command_file_name + character ( len = * ) data_file_name + integer file_unit +! LQ integer i + integer ierror + integer ios +! LQ integer ncol +! +! Write the data file. +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZ_PLOT - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = command_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZ_PLOT - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + write ( file_unit, '(a)' ) 'set title "GNUFOR plot"' + write ( file_unit, '(a)' ) 'set xlabel "x"' + write ( file_unit, '(a)' ) 'set ylabel "y"' + write ( file_unit, '(a)' ) 'set parametric' + write ( file_unit, '(a)' ) 'splot "' // trim ( data_file_name ) // & + '" using 1:2:3 with lines 1' + write ( file_unit, '(a)' ) 'pause -1' + write ( file_unit, '(a)' ) 'q' + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZ_PLOT:' + write ( *, '(a)' ) ' Wrote the GNUPLOT SPLOT command file "' // & + trim ( command_file_name ) // '"' + +end subroutine write_xyz_plot + +subroutine write_xyzgrid_contour ( command_file_name, data_file_name, ierror ) +! +!******************************************************************************* +! +!! WRITE_XYZGRID_CONTOUR writes commands to plot contours of Z(X,Y). +! +! +! Discussion: +! +! This routine tries to write a command file suitable for displaying +! contours of Z(X,Y) gridded data. A grid data file, containing values +! of X, Y and Z, will also be needed. +! +! Modified: +! +! 22 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) COMMAND_FILE_NAME, the name of the +! command file. +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + character ( len = * ) command_file_name + character ( len = * ) data_file_name + integer file_unit +! LQ integer i + integer ierror + integer ios +! LQ integer ncol +! +! Write the data file. +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZGRID_CONTOUR - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = command_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZGRID_CONTOUR - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + write ( file_unit, '(a)' ) 'set title "GNUFOR plot"' + write ( file_unit, '(a)' ) 'set xlabel "x"' + write ( file_unit, '(a)' ) 'set ylabel "y"' + write ( file_unit, '(a)' ) 'set parametric' + write ( file_unit, '(a)' ) 'set nosurface' + write ( file_unit, '(a)' ) 'set contour' + write ( file_unit, '(a)' ) 'set cntrparam levels 10' + write ( file_unit, '(a)' ) 'set term table' + write ( file_unit, '(a)' ) 'set out "table.txt"' + write ( file_unit, '(a)' ) 'splot "' // trim ( data_file_name ) // & + '" using 1:2:3 with lines' + write ( file_unit, '(a)' ) 'set term x11' + write ( file_unit, '(a)' ) 'plot "table.txt" with lines' + write ( file_unit, '(a)' ) 'pause -1' + write ( file_unit, '(a)' ) 'q' + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZGRID_CONTOUR:' + write ( *, '(a)' ) ' Wrote the GNUPLOT XYZGRID contour plot command file "' // & + trim ( command_file_name ) // '"' + +end subroutine write_xyzgrid_contour + +subroutine write_xyzgrid_data ( data_file_name, nx, ny, xyz, ierror ) +! +!******************************************************************************* +! +!! WRITE_XYZGRID_DATA writes a file of XYZ grid data. +! +! +! Discussion: +! +! It is assumed that values of Z are available on a regular NX by NY grid +! of (X,Y) points. +! +! The form of the data file requires that all the data for a given value +! of Y be listed, followed by a blank line, followed by the data for +! another value of Y. +! +! Example: +! +! Here is a grid data file for a 3 by 3 grid, with Z = X + Y. +! +! 0.0 0.0 0.0 +! 1.0 0.0 1.0 +! 2.0 0.0 2.0 +! +! 0.0 1.0 1.0 +! 1.0 1.0 2.0 +! 2.0 1.0 3.0 +! +! 0.0 2.0 2.0 +! 1.0 2.0 3.0 +! 2.0 2.0 4.0 +! +! Modified: +! +! 23 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Input, integer NX, NY, the dimensions of the grid. +! +! Input, real XYZ(3,NX,NY), the XYZ grid data to be written. +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + integer nx + integer ny +! + character ( len = * ) data_file_name + integer file_unit + integer i + integer ierror + integer ios + integer j + real(kind=8) :: xyz(3,nx,ny) +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZGRID_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = data_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZGRID_DATA - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + do j = 1, ny + do i = 1, nx + write ( file_unit, * ) xyz(1:3,i,j) + end do + write ( file_unit, '(a)' ) + end do + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZGRID_DATA:' + write ( *, '(a)' ) ' Wrote the GNUPLOT XYZ grid data file "' // & + trim ( data_file_name ) // '"' + +end subroutine write_xyzgrid_data + +subroutine write_xyzgrid_surface ( command_file_name, data_file_name, ierror ) +! +!******************************************************************************* +! +!! WRITE_XYZGRID_SURFACE writes a file of GNUPLOT commands to plot a 3D surface. +! +! +! Discussion: +! +! This routine tries to write a command file suitable for displaying +! a surface Z(X,Y). A grid data file, containing values of X, Y and Z, +! will also be needed. The routine WRITE_XYZGRID_DATA can write this file. +! +! Modified: +! +! 22 February 2001 +! +! Author: +! +! John Burkardt +! +! Parameters: +! +! Input, character ( len = * ) COMMAND_FILE_NAME, the name of the +! command file. +! +! Input, character ( len = * ) DATA_FILE_NAME, the name of the data file. +! +! Output, integer IERROR, nonzero if an error occurred. +! + implicit none +! + character ( len = * ) command_file_name + character ( len = * ) data_file_name + integer file_unit +! LQ integer i + integer ierror + integer ios +! LQ integer ncol +! +! Write the data file. +! + ierror = 0 + + call get_unit ( file_unit ) + + if ( file_unit == 0 ) then + ierror = 1 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZGRID_SURFACE - Fatal error!' + write ( *, '(a)' ) ' Could not get a free FORTRAN unit.' + return + end if + + open ( unit = file_unit, file = command_file_name, status = 'replace', & + iostat = ios ) + + if ( ios /= 0 ) then + ierror = 2 + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_XYZGRID_SURFACE - Fatal error!' + write ( *, '(a)' ) ' Could not open the output file.' + return + end if + + write ( file_unit, '(a)' ) 'set title "GNUFOR plot"' + write ( file_unit, '(a)' ) 'set xlabel "x"' + write ( file_unit, '(a)' ) 'set ylabel "y"' + write ( file_unit, '(a)' ) 'set parametric' + write ( file_unit, '(a)' ) 'set hidden3d' + write ( file_unit, '(a)' ) 'set contour' + write ( file_unit, '(a)' ) 'splot "' // trim ( data_file_name ) // & + '" using 1:2:3 with lines' + write ( file_unit, '(a)' ) 'pause -1' + write ( file_unit, '(a)' ) 'q' + + close ( unit = file_unit ) + + write ( *, '(a)' ) ' ' + write ( *, '(a)' ) 'WRITE_SURFACE_COMMANDS:' + write ( *, '(a)' ) ' Wrote the GNUPLOT surface plot command file "' // & + trim ( command_file_name ) // '"' + +end subroutine write_xyzgrid_surface + +! LQ AG + +END MODULE gnufor_mio + +! LQ AG diff --git a/inversa.f90 b/inversa.f90 new file mode 100644 index 0000000..86bea31 --- /dev/null +++ b/inversa.f90 @@ -0,0 +1,94 @@ +module inversa + + implicit none +contains + + subroutine inverse(a,c,n) + !============================================================ + ! Inverse matrix + ! Method: Based on Doolittle LU factorization for Ax=b + ! Alex G. December 2009 + !----------------------------------------------------------- + ! input ... + ! a(n,n) - array of coefficients for matrix A + ! n - dimension + ! output ... + ! c(n,n) - inverse matrix of A + ! comments ... + ! the original matrix a(n,n) will be destroyed + ! during the calculation + !=========================================================== + + integer, intent(in) :: n + + real(kind = 8), dimension(n, n), intent(inout) :: a + + real(kind = 8), dimension(n, n), intent(out) :: c + + real(kind = 8), dimension(n,n) :: L, U + + real(kind = 8), dimension(n) :: b, d, x + + real(kind = 8) :: coeff + + integer :: i, j, k + + ! step 0: initialization for matrices L and U and b + ! Fortran 90/95 aloows such operations on matrices + L=0.0 + U=0.0 + b=0.0 + + ! step 1: forward elimination + do k=1, n-1 + do i=k+1,n + coeff=a(i,k)/a(k,k) + L(i,k) = coeff + do j=k+1,n + a(i,j) = a(i,j)-coeff*a(k,j) + end do + end do + end do + + ! Step 2: prepare L and U matrices + ! L matrix is a matrix of the elimination coefficient + ! + the diagonal elements are 1.0 + do i=1,n + L(i,i) = 1.0 + end do + ! U matrix is the upper triangular part of A + do j=1,n + do i=1,j + U(i,j) = a(i,j) + end do + end do + + ! Step 3: compute columns of the inverse matrix C + do k=1,n + b(k)=1.0 + d(1) = b(1) + ! Step 3a: Solve Ld=b using the forward substitution + do i=2,n + d(i)=b(i) + do j=1,i-1 + d(i) = d(i) - L(i,j)*d(j) + end do + end do + ! Step 3b: Solve Ux=d using the back substitution + x(n)=d(n)/U(n,n) + do i = n-1,1,-1 + x(i) = d(i) + do j=n,i+1,-1 + x(i)=x(i)-U(i,j)*x(j) + end do + x(i) = x(i)/u(i,i) + end do + ! Step 3c: fill the solutions x(n) into column k of C + do i=1,n + c(i,k) = x(i) + end do + b(k)=0.0 + end do + end subroutine inverse + +endmodule inversa diff --git a/problema_vero.f90 b/problema_vero.f90 new file mode 100644 index 0000000..476e8fd --- /dev/null +++ b/problema_vero.f90 @@ -0,0 +1,192 @@ +MODULE problema_vero + + ! Modulo contenente i flussi per Shallow-Water Equations in 1D + ! W = vettore variabili conservative = [a; au] + ! u = velocità orizzontale + ! a = h - zb + ! h = altezza colonna d'acqua dallo 0 + ! zb = altezza fondale + + use inversa + + implicit none + + real(kind = 8), parameter :: g = 9.81d0 + + interface flusso_vero + + module procedure flusso_vero_s, flusso_vero_v + + endinterface flusso_vero + + interface sorgente + + module procedure sorgente_s, sorgente_v + + endinterface sorgente + + interface eigenstructure + + module procedure eigenstructure_s, eigenstructure_v + + endinterface eigenstructure + + interface jacobian + + module procedure jacobian_s, jacobian_v + + endinterface jacobian + + +contains + + subroutine flusso_vero_s(w, flusso) + + real(kind = 8), dimension(:), intent(in) :: w + + real(kind = 8), dimension(size(w)), intent(out) :: flusso + + flusso(1) = w(2) ! f(1) = a u + + flusso(2) = (w(2) ** 2) / w(1) + 0.5d0 * g * (w(1) ** 2) ! f(2) = a u^2 + 1/2 g a^2 + + endsubroutine flusso_vero_s + + subroutine flusso_vero_v(ww, fflusso) + + real(kind = 8), dimension(:, :), intent(in) :: ww + + real(kind = 8), dimension(size(ww, 1), size(ww, 2)), intent(out) :: fflusso + + fflusso(1, :) = ww(2, :) ! f(1) = a u + + fflusso(2, :) = (ww(2, :) ** 2) / ww(1, :) + 0.5d0 * g * (ww(1, :) ** 2) ! f(2) = a u^2 + 1/2 g a^2 + + endsubroutine flusso_vero_v + + subroutine sorgente_s(w, cf, d_zb, S) + + real(kind = 8), dimension(:), intent(in) :: w + + real(kind = 8), intent(in) :: cf, d_zb + + real(kind = 8), dimension(size(w)), intent(out) :: S + + S(1) = 0 + + S(2) = cf * (w(2) / w(1)) * sqrt(w(1)**2 + w(2)**2) + g * w(1) * d_zb ! S(2) = cf u |w| + g a dx(zb) + + endsubroutine sorgente_s + + subroutine sorgente_v(ww, ccf, dd_zb, SS) + + real(kind = 8), dimension(:, :), intent(in) :: ww + + real(kind = 8), dimension(:), intent(in) :: ccf, dd_zb + + real(kind = 8), dimension(size(ww, 1), size(ww, 2)), intent(out) :: SS + + SS(1, :) = 0 + + SS(2, :) = ccf * (ww(2, :) / ww(1, :)) * sqrt(ww(1, :)**2 + ww(2, :)**2) +& + g * ww(1, :) * dd_zb ! S(2) = cf u |w| + g a dx(zb) + + endsubroutine sorgente_v + + subroutine jacobian_s(w, J) + + ! Subroutine per il calcolo dello jacobiano in 1D + + real(kind = 8), dimension(:), intent(in) :: w + + real(kind = 8), dimension(size(w), size(w)), intent(out) :: J + + J(1, 1) = 0; J(1, 2) = 1.d0; ! J(1, :) = [ 0, 1 ] + + J(2, 1) = -(w(2) / w(1))**2 + g * w(1) ! J(2, 1) = - u^2 + ga + + J(2, 2) = 2 * w(2) / w(1) ! J(2, 2) = 2u + + + endsubroutine jacobian_s + + subroutine jacobian_v(ww, JJ) + + ! Subroutine per il calcolo dello jacobiano in 1D + + real(kind = 8), dimension(:, :), intent(in) :: ww + + real(kind = 8), dimension(size(ww, 1), size(ww, 1), size(ww, 2)), intent(out) :: JJ + + JJ(1, 1, :) = 0; JJ(1, 2, :) = 1.d0; ! JJ(1, :, :) = [ 0, 1 ] + + JJ(2, 1, :) = -(ww(2, :) / ww(1, :))**2 + g * ww(1, :) ! JJ(2, 1, :) = - u^2 + ga + + JJ(2, 2, :) = 2 * ww(2, :) / ww(1, :) ! JJ(2, 2, :) = 2u + + + + + endsubroutine jacobian_v + + subroutine eigenstructure_s(w, lambda, R, L) + + real(kind = 8), dimension(:), intent(in) :: w + + real(kind = 8), dimension(size(w)), intent(out) :: lambda + + real(kind = 8), dimension(size(w), size(w)), intent(out) :: R, L + + real(kind = 8), dimension(size(w), size(w)) :: B ! Matrice di appoggio per il calcolo dell'inversa + + lambda(1) = w(2) / w(1) - sqrt(g * w(1)) ! lambda(1) = u - sqrt(ga) + lambda(2) = w(2) / w(1) + sqrt(g * w(1)) ! lambda(2) = u + sqrt(ga) + + R(1, 1) = 1.d0; R(1, 2) = 1.d0; ! R(1, :) = [1, 1] + + R(2, 1) = w(2) / w(1) - sqrt(g * w(1)) ! R(2, 1) = u - sqrt(ga) + R(2, 2) = w(2) / w(1) + sqrt(g * w(1)) ! R(2, 2) = u + sqrt(ga) + + !R(2, 1) = - sqrt(g * w(1)) ! R(2, 1) = - sqrt(ga) + !R(2, 2) = sqrt(g * w(1)) ! R(2, 2) = sqrt(ga) + + B = R ! Mi appoggio su B + + call inverse(B, L, size(w)) + + endsubroutine eigenstructure_s + + subroutine eigenstructure_v(ww, llambda, RR, LL) + + real(kind = 8), dimension(:, :), intent(in) :: ww + + real(kind = 8), dimension(size(ww, 1), size(ww, 2)), intent(out) :: llambda + + real(kind = 8), dimension(size(ww, 1), size(ww, 1), size(ww, 2)), intent(out) :: RR, LL + + real(kind = 8), dimension(size(ww, 1), size(ww, 1)) :: B ! Matrice di appoggio per il calcolo dell'inversa + + integer :: j + + llambda(1, :) = ww(2, :) / ww(1, :) - sqrt(g * ww(1, :)) ! lambda(1) = u - sqrt(ga) + llambda(2, :) = ww(2, :) / ww(1, :) + sqrt(g * ww(1, :)) ! lambda(2) = u + sqrt(ga) + + RR(1, 1, :) = 1.d0; RR(1, 2, :) = 1.d0; ! R(1, :) = [1, 1] + + RR(2, 1, :) = ww(2, :) / ww(1, :) - sqrt(g * ww(1, :)) ! R(2, 1) = u - sqrt(ga) + RR(2, 2, :) = ww(2, :) / ww(1, :) + sqrt(g * ww(1, :)) ! R(2, 2) = u + sqrt(ga) + + !RR(2, 1, :) = - sqrt(g * ww(1, :)) ! R(2, 1) = - sqrt(ga) + !RR(2, 2, :) = sqrt(g * ww(1, :)) ! R(2, 2) = sqrt(ga) + + do j = 1, size(ww, 2) + + B = RR(:, :, j) + + call inverse(B, LL(:, :, j), size(ww, 1)) + + enddo + + endsubroutine eigenstructure_v + +ENDMODULE problema_vero diff --git a/slope_calc.f90 b/slope_calc.f90 new file mode 100644 index 0000000..6941182 --- /dev/null +++ b/slope_calc.f90 @@ -0,0 +1,131 @@ +module slope_calc + + ! Calcolo delle pendenze di Upwind e di Lax-Wendroff + + use Boundary_conditions + + implicit none + +contains + + subroutine Up_slope(Np, order, slope_vec) + + ! Calcolo del vettore delle pendenze col metodo Upwind + + integer, intent(in) :: Np, order + + real(kind = 8), dimension(order, order, Np), intent(inout) :: slope_vec + + slope_vec = 0 + + + endsubroutine Up_slope + + + + + subroutine LW_slope(Dx_vec, stato, BC_fis, BC_cons, slope_vec_for, slope_vec_back, V, R) + + ! Calcolo del vettore delle pendenze col metodo Lax-Wendroff + + real(kind = 8), dimension(:), intent(in) :: Dx_vec + + real(kind = 8), dimension(:, :), intent(inout) :: BC_cons, BC_fis, stato + + real(kind = 8), dimension(size(BC_cons, 1), size(BC_cons, 2),& + size(stato, 2)), intent(inout) :: slope_vec_for, slope_vec_back + + real(kind = 8), dimension(size(BC_cons, 1),& + size(stato, 2) - 1) :: Beta_vec_for, Beta_vec_back + + real(kind = 8), dimension(size(BC_cons, 1), size(stato, 2)), intent(inout) :: V + + real(kind = 8), dimension(size(BC_cons, 1)) :: stato_dx, stato_sx + + real(kind = 8), dimension(size(BC_cons, 1), size(BC_cons, 1), size(stato, 2)), intent(inout) :: R + + real(kind = 8), dimension(size(BC_cons, 1), size(BC_cons, 1), size(stato, 2)) :: L + + real(kind = 8), dimension(size(BC_cons, 1), size(BC_cons, 1)) :: L_app, B + + real(kind = 8), dimension(size(BC_cons, 1)) :: R_app + + real(kind = 8), dimension(size(BC_cons, 1), size(stato, 2) - 1) :: delta_V_for, delta_V_back + + real(kind = 8), dimension(size(BC_cons, 1), size(stato, 2) - 1) :: delta_U + + + integer :: dim, j, k + + + delta_U = stato(:, 2:size(stato, 2)) - stato(:, 1:size(stato, 2) - 1) + + ! Mi porto nelle variabili caratteristiche + + do j= 1, size(stato, 2) + + call R_calc(stato(:, j), R(:, :, j)) + + B = R(:, :, j) + + call L_calc(B, L(:, :, j)) + + L_app = L(:, :, j) + + V(:, j) = matmul(L_app, stato(:, j)) + + enddo + + do j= 1, size(stato, 2) - 1 + + L_app = L(:, :, j + 1) + + delta_V_for(:, j) = matmul(L_app, delta_U(:, j)) + + L_app = L(:, :, j) + + delta_V_back(:, j) = matmul(L_app, delta_U(:, j)) + + enddo + + dim = size(stato, 2) + + do k = 1, size(stato, 1) + + Beta_vec_for(k, :) = delta_V_for(k, :)/Dx_vec(1:dim - 1) + + Beta_vec_back(k, :) = delta_V_back(k, :)/Dx_vec(2:dim) + + ! Torno alle pendenze nelle variabili conservative + + do j = 1, dim - 1 + + R_app = R(:, k, j + 1) + slope_vec_for(:, k, j) = R_app * Beta_vec_for(k, j) + + R_app = R(:, k, j) + slope_vec_back(:, k, j + 1) = R_app * Beta_vec_back(k, j) + + enddo + + enddo + + slope_vec_for(:, :, size(slope_vec_for, 3)) = slope_vec_for(:, :, size(slope_vec_for, 3) - 1) + slope_vec_back(:, :, 1) = slope_vec_back(:, :, 2) + + stato_dx = stato(:, size(stato, 2)) + + call BC_dx(stato_dx, BC_fis) + + stato_sx = stato(:, 1) + + call BC_sx(stato_sx, BC_fis) + + call fis2cons(BC_fis, BC_cons) + + + !write(*,*) "BC_dx = ", BC_cons(:, 2) + + endsubroutine LW_slope + +endmodule slope_calc diff --git a/slope_det.f90 b/slope_det.f90 new file mode 100644 index 0000000..a9bbf22 --- /dev/null +++ b/slope_det.f90 @@ -0,0 +1,49 @@ +module slope_det + + ! Modulo per la determinazione della pendenza + + use slope_calc + + use checking1 + + implicit none + +contains + + subroutine slope_fin(Dx_vec, stato, BC_fis, sslope_def, BC_cons, V, R) + + real(kind = 8), dimension(:), intent(in) :: Dx_vec + + real(kind = 8), dimension(:, :), intent(inout) :: stato, BC_fis + + real(kind = 8), dimension(size(stato, 1), size(stato, 1), size(stato, 2)), intent(inout) :: sslope_def, R + + real(kind = 8), dimension(size(stato, 1), size(stato, 1), size(stato, 2)) :: sslope_LW_for,& + sslope_LW_back,& + sslope_Up + + real(kind = 8), dimension(size(BC_fis, 1), size(BC_fis, 2)), intent(inout) :: BC_cons + + real(kind = 8), dimension(size(stato, 1), size(stato, 2)), intent(inout) :: V + + integer :: j + + !----------------------------------------------------------------------------- + ! Calcolo vettore delle pendenze per rho e per m + !----------------------------------------------------------------------------- + + call Up_slope(size(stato, 2), size(stato, 1), sslope_Up) ! Vettore delle pendenze di Upwind (tutte 0) + + ! Vettore delle pendenze di Lax_Wendroff + + call LW_slope(Dx_vec, stato, BC_fis, BC_cons, sslope_LW_for, sslope_LW_back, V, R) + + do j = 1, size(Dx_vec) + + call slope(sslope_LW_for(:, :, j), sslope_LW_back(:, :, j), sslope_def(:, :, j)) + + enddo + +endsubroutine slope_fin + +endmodule slope_det