Skip to content

Commit 6a06a1b

Browse files
pmoczpmocz
authoredMar 4, 2025··
modernizing num tests (#792)
* modernizing num tests * change compile rule * elim .f rules in rates test makefile, format bari* files * updating fortitude rules
1 parent 10bb99b commit 6a06a1b

18 files changed

+1091
-1134
lines changed
 

‎.github/workflows/static-analysis.yml

+1
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ on:
55
workflow_dispatch:
66
release:
77
push:
8+
branches: ["master"]
89
pull_request:
910
types: [opened, synchronize, reopened, closed]
1011

‎atm/test/src/test_atm.f90

+2
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,8 @@ program test_atm
44
use test_atm_support, only: do_test_atm, &
55
test_verbosely, cgrav, eos_handle, kap_handle
66

7+
implicit none
8+
79
logical :: test_verbosely_in
810
real(dp) :: cgrav_in
911
integer :: eos_handle_in, kap_handle_in

‎atm/test/src/test_atm_quietly.f90

+3-2
Original file line numberDiff line numberDiff line change
@@ -4,10 +4,11 @@ program test_atm_quietly
44
use test_atm_support, only: do_test_atm, &
55
test_verbosely, cgrav, eos_handle, kap_handle
66

7+
implicit none
8+
79
call setup
810
test_verbosely = .false.
911

10-
call do_test_atm( &
11-
test_verbosely, cgrav, eos_handle, kap_handle)
12+
call do_test_atm(test_verbosely, cgrav, eos_handle, kap_handle)
1213

1314
end program test_atm_quietly

‎fortitude.toml

+1-1
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,3 @@
11
[check]
2-
ignore = ["E001","OB001","OB011","P001","P011","P021","T001","T002","T011","T031","T041","T042","S001","S101","S102","M001","M011"]
2+
ignore = ["E001","C001","C003","C061","C071","C072","C092","C121","C131","S061","OB011","PORT011"]
33
line-length = 200

‎ionization/test/src/test_ionization.f90

+1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
program test_ionization
22
use test_ionization_support
3+
implicit none
34
call do_test(.false.)
45
end program test_ionization

‎ionization/test/src/test_ionization_quietly.f90

+1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
program test_ionization_quietly
22
use test_ionization_support
3+
implicit none
34
call do_test(.true.)
45
end program test_ionization_quietly

‎num/private/mod_random.f90

+1
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626
module mod_random
2727
use const_def, only: dp
2828

29+
implicit none
2930

3031
contains
3132

‎num/test/make/makefile_base

+14-21
Original file line numberDiff line numberDiff line change
@@ -12,13 +12,13 @@ include $(MESA_DIR)/utils/makefile_header
1212
#
1313
# SOURCES
1414

15-
COMMON_SRCS = test_support.o test_int_support.o \
16-
bari_beam.o test_beam.o \
17-
bari_chemakzo.o test_chemakzo.o \
18-
bari_medakzo.o test_medakzo.o \
19-
bari_vdpol.o test_vdpol.o \
20-
test_diffusion.o test_brent.o test_bobyqa.o test_newuoa.o \
21-
test_simplex.o test_newton.o test_integrate.f
15+
COMMON_SRCS = test_support.f90 test_int_support.f90 \
16+
bari_beam.f90 test_beam.f90 \
17+
bari_chemakzo.f90 test_chemakzo.f90 \
18+
bari_medakzo.f90 test_medakzo.f90 \
19+
bari_vdpol.f90 test_vdpol.f90 \
20+
test_diffusion.f90 test_brent.f90 test_bobyqa.f90 test_newuoa.f90 \
21+
test_simplex.f90 test_newton.f90 test_integrate.f90
2222

2323
TEST_SRCS = $(COMMON_SRCS) test_num.o
2424

@@ -46,9 +46,9 @@ TEST = $(TEST_DIR)/tester
4646
TEST_MPI = $(TEST_DIR)/test_mpi
4747
SAMPLE = $(TEST_DIR)/sample
4848

49-
TEST_OBJS = $(patsubst %.f,%.o,$(patsubst %.f90,%.o,$(TEST_SRCS)))
50-
TEST_MPI_OBJS = $(patsubst %.f,%.o,$(patsubst %.f90,%.o,$(TEST_MPI_SRCS)))
51-
SAMPLE_OBJS = $(patsubst %.f,%.o,$(patsubst %.f90,%.o,$(SAMPLE_SRCS)))
49+
TEST_OBJS = $(patsubst %.f90,%.o,$(TEST_SRCS))
50+
TEST_MPI_OBJS = $(patsubst %.f90,%.o,$(TEST_MPI_SRCS))
51+
SAMPLE_OBJS = $(patsubst %.f90,%.o,$(SAMPLE_SRCS))
5252

5353
all : $(TEST) $(SAMPLE)
5454

@@ -77,20 +77,14 @@ nodeps : $(.DEFAULT_GOAL)
7777
#
7878
# COMPILATION RULES
7979

80-
%.o : %.f
81-
ifneq ($(QUIET),)
82-
@echo TEST_COMPILE fixed $<
83-
@$(TEST_COMPILE) $(FCfixed) -fimplicit-none $<
84-
else
85-
$(TEST_COMPILE) $(FCfixed)-fimplicit-none $<
86-
endif
80+
COMPILE = $(FCbasic) $(FCopenmp) $(FCopt) $(FCfree)
8781

8882
%.o: %.f90
8983
ifneq ($(QUIET),)
90-
@echo TEST_COMPILE free $<
91-
@$(TEST_COMPILE) $(FCfree) -fimplicit-none $<
84+
@echo TEST_COMPILE $<
85+
@$(TEST_COMPILE) $(COMPILE) $<
9286
else
93-
$(TEST_COMPILE) $(FCfree) -fimplicit-none $<
87+
$(TEST_COMPILE) $(COMPILE) $<
9488
endif
9589

9690
#################################################################
@@ -99,7 +93,6 @@ endif
9993

10094
SRC_PATH = $(TEST_SRC_DIR)
10195

102-
vpath %.f $(SRC_PATH)
10396
vpath %.f90 $(SRC_PATH)
10497

10598
vpath %.mod $(LOCAL_LIB_DIR):$(MESA_DIR)/include

‎num/test/src/bari_beam.f

-233
This file was deleted.

‎num/test/src/bari_beam.f90

+228
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,228 @@
1+
! ----------------------------------------------------------------------
2+
! ----------------------------------------------------------------------
3+
!
4+
! This file is part of the Test Set for IVP solvers
5+
! http://www.dm.uniba.it/~testset/
6+
!
7+
! Beam (ODE case)
8+
! ODE of dimension 80
9+
!
10+
! DISCLAIMER: see
11+
! http://www.dm.uniba.it/~testset/disclaimer.php
12+
!
13+
! The most recent version of this source file can be found at
14+
! http://www.dm.uniba.it/~testset/src/problems/beam.f
15+
!
16+
! This is revision
17+
! $Id: beam.F,v 1.2 2006/10/02 10:29:13 testset Exp $
18+
!
19+
! ----------------------------------------------------------------------
20+
subroutine beam_init(neqn, y, yprime, consis)
21+
use const_def, only: dp
22+
integer :: neqn
23+
real(dp) :: y(neqn), yprime(neqn)
24+
logical :: consis
25+
26+
integer :: i
27+
28+
do i = 1, neqn
29+
y(i) = 0d0
30+
end do
31+
32+
end subroutine beam_init
33+
! ----------------------------------------------------------------------
34+
subroutine beam_feval(nvar, t, th, df, ierr, rpar, ipar)
35+
use const_def, only: dp
36+
use math_lib
37+
implicit real(dp) (A - H, O - Z)
38+
integer :: ierr, nvar, i, ipar(*)
39+
integer, parameter :: N = 40, NN = 2*N, NSQ = N*N, NQUATR = NSQ*NSQ
40+
real(dp) :: rpar(*), an, deltas
41+
dimension DF(NN), TH(150), U(150), V(150), W(150)
42+
dimension ALPHA(150), BETA(150), STH(150), CTH(150)
43+
! --- SET DEFAULT VALUES
44+
if (nvar /= nn) stop 'bad nvar for beam_feval'
45+
AN = N
46+
DELTAS = 1.0D+0/AN
47+
! ----- CALCUL DES TH(I) ET DES SIN ET COS -------------
48+
do I = 2, N
49+
THDIFF = TH(I) - TH(I - 1)
50+
STH(I) = sin(THDIFF)
51+
CTH(I) = cos(THDIFF)
52+
end do
53+
! -------- CALCUL DU COTE DROIT DU SYSTEME LINEAIRE -----
54+
IF (T > 3.14159265358979324D0) THEN
55+
! --------- CASE T GREATER PI ------------
56+
! ---------- I=1 ------------
57+
TERM1 = (-3.D0*TH(1) + TH(2))*NQUATR
58+
V(1) = TERM1
59+
! -------- I=2,..,N-1 -----------
60+
DO I = 2, N - 1
61+
TERM1 = (TH(I - 1) - 2.D0*TH(I) + TH(I + 1))*NQUATR
62+
V(I) = TERM1
63+
end do
64+
! ----------- I=N -------------
65+
TERM1 = (TH(N - 1) - TH(N))*NQUATR
66+
V(N) = TERM1
67+
else
68+
! --------- CASE T LESS EQUAL PI ------------
69+
FABS = 1.5D0*sin(T)*sin(T)
70+
FX = -FABS
71+
FY = FABS
72+
! ---------- I=1 ------------
73+
TERM1 = (-3.D0*TH(1) + TH(2))*NQUATR
74+
TERM2 = NSQ*(FY*cos(TH(1)) - FX*sin(TH(1)))
75+
V(1) = TERM1 + TERM2
76+
! -------- I=2,..,N-1 -----------
77+
DO I = 2, N - 1
78+
TERM1 = (TH(I - 1) - 2.D0*TH(I) + TH(I + 1))*NQUATR
79+
TERM2 = NSQ*(FY*cos(TH(I)) - FX*sin(TH(I)))
80+
V(I) = TERM1 + TERM2
81+
end do
82+
! ----------- I=N -------------
83+
TERM1 = (TH(N - 1) - TH(N))*NQUATR
84+
TERM2 = NSQ*(FY*cos(TH(N)) - FX*sin(TH(N)))
85+
V(N) = TERM1 + TERM2
86+
end if
87+
! -------- COMPUTE PRODUCT DV=W ------------
88+
W(1) = STH(2)*V(2)
89+
DO I = 2, N - 1
90+
W(I) = -STH(I)*V(I - 1) + STH(I + 1)*V(I + 1)
91+
end do
92+
W(N) = -STH(N)*V(N - 1)
93+
! -------- TERME 3 -----------------
94+
DO I = 1, N
95+
W(I) = W(I) + TH(N + I)*TH(N + I)
96+
end do
97+
! ------- SOLVE SYSTEM CW=W ---------
98+
ALPHA(1) = 1.D0
99+
DO I = 2, N
100+
ALPHA(I) = 2.D0
101+
BETA(I - 1) = -CTH(I)
102+
end do
103+
ALPHA(N) = 3.D0
104+
DO I = N - 1, 1, -1
105+
Q = BETA(I)/ALPHA(I + 1)
106+
W(I) = W(I) - W(I + 1)*Q
107+
ALPHA(I) = ALPHA(I) - BETA(I)*Q
108+
end do
109+
W(1) = W(1)/ALPHA(1)
110+
DO I = 2, N
111+
W(I) = (W(I) - BETA(I - 1)*W(I - 1))/ALPHA(I)
112+
end do
113+
! -------- COMPUTE U=CV+DW ---------
114+
U(1) = V(1) - CTH(2)*V(2) + STH(2)*W(2)
115+
DO I = 2, N - 1
116+
U(I) = 2.D0*V(I) - CTH(I)*V(I - 1) - CTH(I + 1)*V(I + 1) - STH(I)*W(I - 1) + STH(I + 1)*W(I + 1)
117+
end do
118+
U(N) = 3.D0*V(N) - CTH(N)*V(N - 1) - STH(N)*W(N - 1)
119+
! -------- PUT DERIVATIVES IN RIGHT PLACE -------------
120+
DO I = 1, N
121+
DF(I) = TH(N + I)
122+
DF(N + I) = U(I)
123+
end do
124+
end subroutine beam_feval
125+
! ----------------------------------------------------------------------
126+
subroutine beam_jeval(ldim, neqn, t, y, yprime, dfdy, ierr, rpar, ipar)
127+
use const_def, only: dp
128+
integer :: ldim, neqn, ierr, ipar(*)
129+
real(dp) :: t, y(neqn), yprime(neqn), dfdy(ldim, neqn), rpar(*)
130+
!
131+
! dummy subroutine
132+
!
133+
end subroutine beam_jeval
134+
! ----------------------------------------------------------------------
135+
subroutine beam_solut(neqn, t, y)
136+
use const_def, only: dp
137+
integer :: neqn
138+
real(dp), intent(in) :: t
139+
real(dp), intent(out) :: y(neqn)
140+
141+
! computed using real(dp) RADAU on an
142+
! Alphaserver DS20E, with a 667 MHz EV67 processor.
143+
!
144+
! uround = 1.01d-19
145+
! rtol = atol = h0 = 1.1d-18
146+
147+
y(1) = -0.5792366591285007D-002
148+
y(2) = -0.1695298550721735D-001
149+
y(3) = -0.2769103312973094D-001
150+
y(4) = -0.3800815655879158D-001
151+
y(5) = -0.4790616859743763D-001
152+
y(6) = -0.5738710435274594D-001
153+
y(7) = -0.6645327313454617D-001
154+
y(8) = -0.7510730581979037D-001
155+
y(9) = -0.8335219765414992D-001
156+
y(10) = -0.9119134654647947D-001
157+
y(11) = -0.9862858700132091D-001
158+
y(12) = -0.1056682200378002D+000
159+
y(13) = -0.1123150395409595D+000
160+
y(14) = -0.1185743552727245D+000
161+
y(15) = -0.1244520128755561D+000
162+
y(16) = -0.1299544113264161D+000
163+
y(17) = -0.1350885180610398D+000
164+
y(18) = -0.1398618819194422D+000
165+
y(19) = -0.1442826441015242D+000
166+
y(20) = -0.1483595472463012D+000
167+
y(21) = -0.1521019429001447D+000
168+
y(22) = -0.1555197978061129D+000
169+
y(23) = -0.1586236993420229D+000
170+
y(24) = -0.1614248603702127D+000
171+
y(25) = -0.1639351238193223D+000
172+
y(26) = -0.1661669673440852D+000
173+
y(27) = -0.1681335081778558D+000
174+
y(28) = -0.1698485080602243D+000
175+
y(29) = -0.1713263782440888D+000
176+
y(30) = -0.1725821847462537D+000
177+
y(31) = -0.1736316537975654D+000
178+
y(32) = -0.1744911773840049D+000
179+
y(33) = -0.1751778187863392D+000
180+
y(34) = -0.1757093178712902D+000
181+
y(35) = -0.1761040960228576D+000
182+
y(36) = -0.1763812607175549D+000
183+
y(37) = -0.1765606097564671D+000
184+
y(38) = -0.1766626352260565D+000
185+
y(39) = -0.1767085270807460D+000
186+
y(40) = -0.1767201761075488D+000
187+
y(41) = 0.3747362681329794D-001
188+
y(42) = 0.1099117880217593D+000
189+
y(43) = 0.1798360474312799D+000
190+
y(44) = 0.2472427305442391D+000
191+
y(45) = 0.3121293820596567D+000
192+
y(46) = 0.3744947377019500D+000
193+
y(47) = 0.4343386073492798D+000
194+
y(48) = 0.4916620354601748D+000
195+
y(49) = 0.5464677854586807D+000
196+
y(50) = 0.5987609702624270D+000
197+
y(51) = 0.6485493611110740D+000
198+
y(52) = 0.6958435169132503D+000
199+
y(53) = 0.7406572668520808D+000
200+
y(54) = 0.7830081747813241D+000
201+
y(55) = 0.8229176659201515D+000
202+
y(56) = 0.8604110305190560D+000
203+
y(57) = 0.8955175502377805D+000
204+
y(58) = 0.9282708263127953D+000
205+
y(59) = 0.9587089334522034D+000
206+
y(60) = 0.9868747821728363D+000
207+
y(61) = 0.1012816579961883D+001
208+
y(62) = 0.1036587736679858D+001
209+
y(63) = 0.1058246826481355D+001
210+
y(64) = 0.1077857811433353D+001
211+
y(65) = 0.1095490222005369D+001
212+
y(66) = 0.1111219164294898D+001
213+
y(67) = 0.1125125269286501D+001
214+
y(68) = 0.1137294526609229D+001
215+
y(69) = 0.1147818025153607D+001
216+
y(70) = 0.1156792132004482D+001
217+
y(71) = 0.1164318845130183D+001
218+
y(72) = 0.1170505992596124D+001
219+
y(73) = 0.1175467424299550D+001
220+
y(74) = 0.1179323003228859D+001
221+
y(75) = 0.1182198586299667D+001
222+
y(76) = 0.1184226111223146D+001
223+
y(77) = 0.1185543909805575D+001
224+
y(78) = 0.1186297084203716D+001
225+
y(79) = 0.1186637618908127D+001
226+
y(80) = 0.1186724615113034D+001
227+
228+
end subroutine beam_solut

‎num/test/src/bari_chemakzo.f

-242
This file was deleted.

‎num/test/src/bari_chemakzo.f90

+212
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,212 @@
1+
!-----------------------------------------------------------------------
2+
!-----------------------------------------------------------------------
3+
!
4+
! This file is part of the Test Set for IVP solvers
5+
! http://www.dm.uniba.it/~testset/
6+
!
7+
! Chemical Akzo Nobel problem
8+
! ODE of dimension 6
9+
!
10+
! DISCLAIMER: see
11+
! http://www.dm.uniba.it/~testset/disclaimer.php
12+
!
13+
! The most recent version of this source file can be found at
14+
! http://www.dm.uniba.it/~testset/src/problems/chemakzo.f
15+
!
16+
! This is revision
17+
! $Id: chemakzo.F,v 1.2 2006/10/02 10:29:14 testset Exp $
18+
!
19+
!-----------------------------------------------------------------------
20+
subroutine chemakzo_prob(fullnm, problm, type, neqn, ndisc, t, numjac, mljac, mujac, nummas, mlmas, mumas, ind)
21+
character(len=*) :: fullnm, problm, type
22+
integer :: neqn, ndisc, mljac, mujac, mlmas, mumas, ind(*)
23+
double precision :: t(0:*)
24+
logical :: numjac, nummas
25+
26+
integer :: i
27+
28+
fullnm = 'Chemical Akzo Nobel problem'
29+
problm = 'chemakzo'
30+
type = 'DAE'
31+
neqn = 6
32+
ndisc = 0
33+
t(0) = 0d0
34+
t(1) = 180d0
35+
numjac = .false.
36+
mljac = neqn
37+
mujac = neqn
38+
mlmas = 0
39+
mumas = 0
40+
do i = 1, neqn
41+
ind(i) = 0
42+
end do
43+
44+
end subroutine chemakzo_prob
45+
!-----------------------------------------------------------------------
46+
subroutine chemakzo_init(neqn, y, yprime, consis)
47+
integer :: neqn
48+
double precision :: y(neqn), yprime(neqn)
49+
logical :: consis
50+
51+
integer :: ierr, ipar(0)
52+
double precision :: rpar(0)
53+
54+
double precision :: ks
55+
parameter(ks=115.83d0)
56+
57+
y(1) = 0.444d0
58+
y(2) = 0.00123d0
59+
y(3) = 0d0
60+
y(4) = 0.007d0
61+
y(5) = 0d0
62+
y(6) = ks*y(1)*y(4)
63+
64+
consis = .true.
65+
66+
call chemakzo_feval(neqn, 0d0, y, y, yprime, ierr, rpar, ipar)
67+
68+
end subroutine chemakzo_init
69+
!-----------------------------------------------------------------------
70+
subroutine chemakzo_feval(neqn, t, y, yprime, f, ierr, rpar, ipar)
71+
integer :: neqn, ierr, ipar(*)
72+
double precision :: t, y(neqn), yprime(neqn), f(neqn), rpar(*)
73+
74+
double precision :: k1, k2, k3, k4, kbig, kla, po2, hen, ks
75+
parameter(k1=18.7d0, k2=0.58d0, k3=0.09d0, k4=0.42d0, kbig=34.4d0, kla=3.3d0, ks=115.83d0, po2=0.9d0, hen=737d0)
76+
double precision :: r1, r2, r3, r4, r5, fin, sqy2
77+
include 'formats'
78+
79+
if (y(2) < 0d0) then
80+
!write(*,*) 'y(2)', y(2)
81+
ierr = -1
82+
return
83+
end if
84+
85+
sqy2 = sqrt(y(2))
86+
r1 = k1*y(1)*y(1)*y(1)*y(1)*sqy2
87+
r2 = k2*y(3)*y(4)
88+
r3 = k2/kbig*y(1)*y(5)
89+
r4 = k3*y(1)*y(4)*y(4)
90+
r5 = k4*(y(6)**2)*sqy2
91+
fin = kla*(po2/hen - y(2))
92+
93+
f(1) = ((-2d0*r1 + r2) - r3) - r4
94+
f(2) = ((-0.5d0*r1 - r4) - 0.5d0*r5) + fin
95+
f(3) = (r1 - r2) + r3
96+
f(4) = (-r2 + r3) - 2d0*r4
97+
f(5) = (r2 - r3) + r5
98+
f(6) = ks*y(1)*y(4) - y(6)
99+
100+
end subroutine chemakzo_feval
101+
!-----------------------------------------------------------------------
102+
subroutine chemakzo_jeval(ldim, neqn, t, y, yprime, dfdy, ierr, rpar, ipar)
103+
integer :: ldim, neqn, ierr, ipar(*)
104+
double precision :: t, y(neqn), yprime(neqn), dfdy(ldim, neqn), rpar(*)
105+
106+
double precision :: k1, k2, k3, k4, kbig, kla, ks
107+
parameter(k1=18.7d0, k2=0.58d0, k3=0.09d0, k4=0.42d0, kbig=34.4d0, kla=3.3d0, ks=115.83d0)
108+
integer :: i, j
109+
double precision :: r11, r12, r23, r24, r31, r35, r41, r44, r52, r56, fin2
110+
double precision :: y13
111+
112+
if (y(2) < 0d0) then
113+
ierr = -1
114+
return
115+
end if
116+
117+
do j = 1, neqn
118+
do i = 1, neqn
119+
dfdy(i, j) = 0d0
120+
end do
121+
end do
122+
123+
y13 = y(1)*y(1)*y(1)
124+
r11 = 4d0*k1*y13*sqrt(y(2))
125+
r12 = 0.5d0*k1*y(1)*y13/sqrt(y(2))
126+
r23 = k2*y(4)
127+
r24 = k2*y(3)
128+
r31 = (k2/kbig)*y(5)
129+
r35 = (k2/kbig)*y(1)
130+
r41 = k3*y(4)**2
131+
r44 = 2d0*k3*y(1)*y(4)
132+
r52 = 0.5d0*k4*(y(6)**2)/sqrt(y(2))
133+
r56 = 2d0*k4*y(6)*sqrt(y(2))
134+
fin2 = -kla
135+
136+
dfdy(1, 1) = (-2d0*r11 - r31) - r41
137+
dfdy(1, 2) = -2d0*r12
138+
dfdy(1, 3) = r23
139+
dfdy(1, 4) = r24 - r44
140+
dfdy(1, 5) = -r35
141+
dfdy(2, 1) = -0.5d0*r11 - r41
142+
dfdy(2, 2) = -0.5d0*r12 - 0.5d0*r52 + fin2
143+
dfdy(2, 4) = -r44
144+
dfdy(2, 6) = -0.5d0*r56
145+
dfdy(3, 1) = r11 + r31
146+
dfdy(3, 2) = r12
147+
dfdy(3, 3) = -r23
148+
dfdy(3, 4) = -r24
149+
dfdy(3, 5) = r35
150+
dfdy(4, 1) = r31 - 2d0*r41
151+
dfdy(4, 3) = -r23
152+
dfdy(4, 4) = -r24 - 2d0*r44
153+
dfdy(4, 5) = r35
154+
dfdy(5, 1) = -r31
155+
dfdy(5, 2) = r52
156+
dfdy(5, 3) = r23
157+
dfdy(5, 4) = r24
158+
dfdy(5, 5) = -r35
159+
dfdy(5, 6) = r56
160+
dfdy(6, 1) = ks*y(4)
161+
dfdy(6, 4) = ks*y(1)
162+
dfdy(6, 6) = -1d0
163+
164+
end subroutine chemakzo_jeval
165+
!-----------------------------------------------------------------------
166+
subroutine chemakzo_meval(ldim, neqn, t, yprime, dfddy, ierr, rpar, ipar)
167+
integer :: ldim, neqn, ierr, ipar(*)
168+
double precision :: t, yprime(neqn), dfddy(ldim, neqn), rpar(*)
169+
170+
integer :: i
171+
ierr = 0
172+
173+
do i = 1, neqn - 1
174+
dfddy(1, i) = 1d0
175+
end do
176+
177+
dfddy(1, neqn) = 0d0
178+
179+
end subroutine chemakzo_meval
180+
!-----------------------------------------------------------------------
181+
subroutine chemakzo_solut(neqn, t, y)
182+
integer :: neqn
183+
double precision :: t, y(neqn)
184+
185+
! computed using true double precision on a Cray C90
186+
187+
! Solving Chemical Akzo Nobel problem using PSIDE
188+
189+
! User input:
190+
191+
! give relative error tolerance: 1d-19
192+
! give absolute error tolerance: 1d-19
193+
194+
! Integration characteristics:
195+
196+
! number of integration steps 5535
197+
! number of accepted steps 5534
198+
! number of f evaluations 100411
199+
! number of Jacobian evaluations 7
200+
! number of LU decompositions 368
201+
202+
! CPU-time used: 30.77 sec
203+
!
204+
205+
y(1) = 0.1150794920661702d0
206+
y(2) = 0.1203831471567715d-2
207+
y(3) = 0.1611562887407974d0
208+
y(4) = 0.3656156421249283d-3
209+
y(5) = 0.1708010885264404d-1
210+
y(6) = 0.4873531310307455d-2
211+
212+
end subroutine chemakzo_solut

‎num/test/src/bari_medakzo.f

-556
This file was deleted.

‎num/test/src/bari_medakzo.f90

+550
Large diffs are not rendered by default.

‎num/test/src/bari_vdpol.f

-58
This file was deleted.

‎num/test/src/bari_vdpol.f90

+58
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
!-----------------------------------------------------------------------
2+
! This file is part of the Test Set for IVP solvers
3+
! http://www.dm.uniba.it/~testset/
4+
!
5+
! Problem VAN DER POL
6+
! ODE of dimension 2
7+
!
8+
! DISCLAIMER: see
9+
! http://www.dm.uniba.it/~testset/disclaimer.php
10+
!
11+
! The most recent version of this source file can be found at
12+
! http://www.dm.uniba.it/~testset/src/problems/vdpol.f
13+
!
14+
! This is revision
15+
! $Id: vdpol.F,v 1.2 2006/10/02 10:29:14 testset Exp $
16+
!
17+
!-----------------------------------------------------------------------
18+
19+
subroutine vdpol_init(neqn, y, yprime, consis)
20+
integer :: neqn
21+
double precision :: y(neqn), yprime(neqn)
22+
logical :: consis
23+
24+
y(1) = 2d0
25+
y(2) = 0d0
26+
27+
end subroutine vdpol_init
28+
29+
subroutine vdpol_feval(neqn, t, y, yprime, f, ierr, rpar, ipar)
30+
integer :: neqn, ierr, ipar(*)
31+
double precision :: t, y(neqn), yprime(neqn), f(neqn), rpar(*)
32+
33+
f(1) = y(2)
34+
f(2) = ((1 - y(1)*y(1))*y(2) - y(1))/1.0d-3
35+
36+
end subroutine vdpol_feval
37+
38+
subroutine vdpol_jeval(ldim, neqn, t, y, yprime, dfdy, ierr, rpar, ipar)
39+
integer :: ldim, neqn, ierr, ipar(*)
40+
double precision :: t, y(neqn), yprime(neqn), dfdy(ldim, neqn), rpar(*)
41+
42+
dfdy(1, 1) = 0d0
43+
dfdy(1, 2) = 1d0
44+
dfdy(2, 1) = (-2.0d0*y(1)*y(2) - 1d0)/1.0d-3
45+
dfdy(2, 2) = (1d0 - y(1)*y(1))/1.0d-3
46+
47+
end subroutine vdpol_jeval
48+
49+
subroutine vdpol_solut(neqn, t, y)
50+
integer :: neqn
51+
double precision, intent(in) :: t
52+
double precision, intent(out) :: y(neqn)
53+
54+
! note -- this is for stiffness param = 1d-3
55+
y(1) = 1.7632345401889102d+00
56+
y(2) = -8.3568868191466206d-01
57+
58+
end subroutine vdpol_solut

‎rates/test/make/makefile_base

+2-2
Original file line numberDiff line numberDiff line change
@@ -36,8 +36,8 @@ TEST = $(TEST_DIR)/tester
3636
TESTQ = $(TEST_DIR)/test_quietly
3737
SHOW_RATES = $(TEST_DIR)/show_rates
3838

39-
TEST_OBJS = $(patsubst %.f,%.o,$(patsubst %.f90,%.o,$(TEST_SRCS)))
40-
SHOW_RATES_OBJS = $(patsubst %.f,%.o,$(patsubst %.f90,%.o,$(SHOW_RATES_SRCS)))
39+
TEST_OBJS = $(patsubst %.f90,%.o,$(patsubst %.c,%.o,$(TEST_SRCS)))
40+
SHOW_RATES_OBJS = $(patsubst %.f90,%.o,$(patsubst %.c,%.o,$(SHOW_RATES_SRCS)))
4141

4242
all : $(TEST) $(TESTQ) $(SHOW_RATES)
4343

‎rates/test/src/test_weak.f90

+17-19
Original file line numberDiff line numberDiff line change
@@ -77,30 +77,30 @@ subroutine do_test_weak
7777
d_eta_dlnRho = 0
7878

7979
if (.false.) then ! TESTING
80-
logT = 7.5904236599874348D+00
81-
logRho = 1.0657946486820271D+00
82-
ye = 8.2724691280321605D-01
83-
eta = -5.3262903257381922D+00
84-
d_eta_dlnT = -1.5299344982339016D+00
85-
d_eta_dlnRho = 9.9482489248846617D-01
80+
logT = 7.5904236599874348D+00
81+
logRho = 1.0657946486820271D+00
82+
ye = 8.2724691280321605D-01
83+
eta = -5.3262903257381922D+00
84+
d_eta_dlnT = -1.5299344982339016D+00
85+
d_eta_dlnRho = 9.9482489248846617D-01
8686
else if (.true.) then ! TESTING
8787
logT = log10(9.0d8)
8888
ye = 0.5d0
8989
logRho = log10(4.5d5)
9090
eta = 10d0
9191

9292
! call for cell 565
93-
! logT = 8.3534130765231005
94-
! logRho = 2.4507395003327828
95-
! T = 225638433.79026267
96-
! Rho = 282.31860559981624
97-
! abar = 4.0424973056746829
98-
! zbar = 2.0238390731055702
99-
! z2bar = 4.2987387813744071
100-
! ye = 0.50064079703023978
101-
! eta = -5.3287260155378711
102-
! d_eta_dlnT= -1.5713400060794886
103-
! d_eta_dlnRho = 1.0016532307086357
93+
! logT = 8.3534130765231005
94+
! logRho = 2.4507395003327828
95+
! T = 225638433.79026267
96+
! Rho = 282.31860559981624
97+
! abar = 4.0424973056746829
98+
! zbar = 2.0238390731055702
99+
! z2bar = 4.2987387813744071
100+
! ye = 0.50064079703023978
101+
! eta = -5.3287260155378711
102+
! d_eta_dlnT= -1.5713400060794886
103+
! d_eta_dlnRho = 1.0016532307086357
104104

105105
else
106106
logT = 7.5d0
@@ -306,6 +306,4 @@ end function dfridr_weak_reaction_info
306306

307307
end subroutine do_test_weak
308308

309-
310309
end module test_weak
311-

0 commit comments

Comments
 (0)
Please sign in to comment.