Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ode_solvers: simplified macros (dropped hidden index ci) #8

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

ode_solvers: simplified macros (dropped hidden index ci) #8

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
167 changes: 85 additions & 82 deletions src/util/ode_solvers_template.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@

! Non-spatial system. ODE solvers operate on vector with state only.
#define _SIZE_ numc
#define _INCC_(m,i) (m)
#define _INPP_(m,n,i) (m,n)
#define _INCCMAP_(m,ci) (m)
#define _INPPMAP_(m,n,ci) (m,n)
#define _ODE_DECLARE_ITERATOR_
#define _ODE_LOOP_BEGIN_
#define _ODE_LOOP_END_
Expand All @@ -15,17 +15,20 @@

! One spatial dimension. ODE solvers operate on matrix with state and space dimensions.
#define _SIZE_ numc,ni
#define _INCC_(m,i) (i,m)
#define _INPP_(m,n,i) (i,m,n)
#define _INCCMAP_(m,ci) (ci,m)
#define _INPPMAP_(m,n,ci) (ci,m,n)
#define _ODE_DECLARE_ITERATOR_ integer :: ci
#define _ODE_LOOP_BEGIN_ do ci=1,ni
#define _ODE_LOOP_END_ end do

#endif

#define _INCC_(m) _INCCMAP_(m,ci)
#define _INPP_(m,n) _INPPMAP_(m,n,ci)

! Dimension specifications for procedure arguments
#define _DIMCC_ _INCC_(1:numc,_LOWI_:ni)
#define _DIMPP_ _INPP_(1:numc,1:numc,_LOWI_:ni)
#define _DIMCC_ _INCCMAP_(1:numc,_LOWI_:ni)
#define _DIMPP_ _INPPMAP_(1:numc,1:numc,_LOWI_:ni)

! Dimension specifications for automatic (temporary) arrays
#define _LOCDIMCC_ _DIMCC_
Expand Down Expand Up @@ -344,10 +347,10 @@ subroutine patankar()
ppsum=_ZERO_
ddsum=_ZERO_
do j=1,numc
ppsum=ppsum+pp _INPP_(i,j,ci)
ddsum=ddsum+dd _INPP_(i,j,ci)
ppsum=ppsum+pp _INPP_(i,j)
ddsum=ddsum+dd _INPP_(i,j)
end do
cc _INCC_(i,ci)=(cc _INCC_(i,ci)+dt*ppsum)/(_ONE_+dt*ddsum/cc _INCC_(i,ci))
cc _INCC_(i)=(cc _INCC_(i)+dt*ppsum)/(_ONE_+dt*ddsum/cc _INCC_(i))
end do
_ODE_LOOP_END_

Expand Down Expand Up @@ -404,13 +407,13 @@ subroutine patankar_runge_kutta_2()

_ODE_LOOP_BEGIN_
do i=1,numc
ppsum _INCC_(i,ci)=_ZERO_
ddsum _INCC_(i,ci)=_ZERO_
ppsum _INCC_(i)=_ZERO_
ddsum _INCC_(i)=_ZERO_
do j=1,numc
ppsum _INCC_(i,ci)=ppsum _INCC_(i,ci)+pp _INPP_(i,j,ci)
ddsum _INCC_(i,ci)=ddsum _INCC_(i,ci)+dd _INPP_(i,j,ci)
ppsum _INCC_(i)=ppsum _INCC_(i)+pp _INPP_(i,j)
ddsum _INCC_(i)=ddsum _INCC_(i)+dd _INPP_(i,j)
end do
cc1 _INCC_(i,ci)=(cc _INCC_(i,ci)+dt*ppsum _INCC_(i,ci))/(_ONE_+dt*ddsum _INCC_(i,ci)/cc _INCC_(i,ci))
cc1 _INCC_(i)=(cc _INCC_(i)+dt*ppsum _INCC_(i))/(_ONE_+dt*ddsum _INCC_(i)/cc _INCC_(i))
end do
_ODE_LOOP_END_

Expand All @@ -419,10 +422,10 @@ subroutine patankar_runge_kutta_2()
_ODE_LOOP_BEGIN_
do i=1,numc
do j=1,numc
ppsum _INCC_(i,ci)=ppsum _INCC_(i,ci)+pp _INPP_(i,j,ci)
ddsum _INCC_(i,ci)=ddsum _INCC_(i,ci)+dd _INPP_(i,j,ci)
ppsum _INCC_(i)=ppsum _INCC_(i)+pp _INPP_(i,j)
ddsum _INCC_(i)=ddsum _INCC_(i)+dd _INPP_(i,j)
end do
cc _INCC_(i,ci)=(cc _INCC_(i,ci)+dt/2*ppsum _INCC_(i,ci))/(_ONE_+dt/2*ddsum _INCC_(i,ci)/cc1 _INCC_(i,ci))
cc _INCC_(i)=(cc _INCC_(i)+dt/2*ppsum _INCC_(i))/(_ONE_+dt/2*ddsum _INCC_(i)/cc1 _INCC_(i))
end do
_ODE_LOOP_END_

Expand Down Expand Up @@ -456,57 +459,57 @@ subroutine patankar_runge_kutta_4()

_ODE_LOOP_BEGIN_
do i=1,numc
ppsum _INCC_(i,ci)=_ZERO_
ddsum _INCC_(i,ci)=_ZERO_
ppsum _INCC_(i)=_ZERO_
ddsum _INCC_(i)=_ZERO_
do j=1,numc
ppsum _INCC_(i,ci)=ppsum _INCC_(i,ci)+pp _INPP_(i,j,ci)
ddsum _INCC_(i,ci)=ddsum _INCC_(i,ci)+dd _INPP_(i,j,ci)
ppsum _INCC_(i)=ppsum _INCC_(i)+pp _INPP_(i,j)
ddsum _INCC_(i)=ddsum _INCC_(i)+dd _INPP_(i,j)
end do
cc1 _INCC_(i,ci)=(cc _INCC_(i,ci)+dt*ppsum _INCC_(i,ci))/(_ONE_+dt*ddsum _INCC_(i,ci)/cc _INCC_(i,ci))
cc1 _INCC_(i)=(cc _INCC_(i)+dt*ppsum _INCC_(i))/(_ONE_+dt*ddsum _INCC_(i)/cc _INCC_(i))
end do
_ODE_LOOP_END_

call get_ppdd(.false.,_SIZE_,cc1,pp,dd)

_ODE_LOOP_BEGIN_
do i=1,numc
ppsum1 _INCC_(i,ci)=_ZERO_
ddsum1 _INCC_(i,ci)=_ZERO_
ppsum1 _INCC_(i)=_ZERO_
ddsum1 _INCC_(i)=_ZERO_
do j=1,numc
ppsum1 _INCC_(i,ci)=ppsum1 _INCC_(i,ci)+pp _INPP_(i,j,ci)
ddsum1 _INCC_(i,ci)=ddsum1 _INCC_(i,ci)+dd _INPP_(i,j,ci)
ppsum1 _INCC_(i)=ppsum1 _INCC_(i)+pp _INPP_(i,j)
ddsum1 _INCC_(i)=ddsum1 _INCC_(i)+dd _INPP_(i,j)
end do
cc1 _INCC_(i,ci)=(cc _INCC_(i,ci)+dt*ppsum1 _INCC_(i,ci))/(_ONE_+dt*ddsum1 _INCC_(i,ci)/cc1 _INCC_(i,ci))
cc1 _INCC_(i)=(cc _INCC_(i)+dt*ppsum1 _INCC_(i))/(_ONE_+dt*ddsum1 _INCC_(i)/cc1 _INCC_(i))
end do
_ODE_LOOP_END_

call get_ppdd(.false.,_SIZE_,cc1,pp,dd)

_ODE_LOOP_BEGIN_
do i=1,numc
ppsum2 _INCC_(i,ci)=_ZERO_
ddsum2 _INCC_(i,ci)=_ZERO_
ppsum2 _INCC_(i)=_ZERO_
ddsum2 _INCC_(i)=_ZERO_
do j=1,numc
ppsum2 _INCC_(i,ci)=ppsum2 _INCC_(i,ci)+pp _INPP_(i,j,ci)
ddsum2 _INCC_(i,ci)=ddsum2 _INCC_(i,ci)+dd _INPP_(i,j,ci)
ppsum2 _INCC_(i)=ppsum2 _INCC_(i)+pp _INPP_(i,j)
ddsum2 _INCC_(i)=ddsum2 _INCC_(i)+dd _INPP_(i,j)
end do
cc1 _INCC_(i,ci)=(cc _INCC_(i,ci)+dt*ppsum2 _INCC_(i,ci))/(_ONE_+dt*ddsum2 _INCC_(i,ci)/cc1 _INCC_(i,ci))
cc1 _INCC_(i)=(cc _INCC_(i)+dt*ppsum2 _INCC_(i))/(_ONE_+dt*ddsum2 _INCC_(i)/cc1 _INCC_(i))
end do
_ODE_LOOP_END_

call get_ppdd(.false.,_SIZE_,cc1,pp,dd)

_ODE_LOOP_BEGIN_
do i=1,numc
ppsum3 _INCC_(i,ci)=_ZERO_
ddsum3 _INCC_(i,ci)=_ZERO_
ppsum3 _INCC_(i)=_ZERO_
ddsum3 _INCC_(i)=_ZERO_
do j=1,numc
ppsum3 _INCC_(i,ci)=ppsum3 _INCC_(i,ci)+pp _INPP_(i,j,ci)
ddsum3 _INCC_(i,ci)=ddsum3 _INCC_(i,ci)+dd _INPP_(i,j,ci)
ppsum3 _INCC_(i)=ppsum3 _INCC_(i)+pp _INPP_(i,j)
ddsum3 _INCC_(i)=ddsum3 _INCC_(i)+dd _INPP_(i,j)
end do
ppsum _INCC_(i,ci)=(ppsum _INCC_(i,ci)/2+ppsum1 _INCC_(i,ci)+ppsum2 _INCC_(i,ci)+ppsum3 _INCC_(i,ci)/2)/3
ddsum _INCC_(i,ci)=(ddsum _INCC_(i,ci)/2+ddsum1 _INCC_(i,ci)+ddsum2 _INCC_(i,ci)+ddsum3 _INCC_(i,ci)/2)/3
cc _INCC_(i,ci)=(cc _INCC_(i,ci)+dt*ppsum _INCC_(i,ci))/(_ONE_+dt*ddsum _INCC_(i,ci)/cc1 _INCC_(i,ci))
ppsum _INCC_(i)=(ppsum _INCC_(i)/2+ppsum1 _INCC_(i)+ppsum2 _INCC_(i)+ppsum3 _INCC_(i)/2)/3
ddsum _INCC_(i)=(ddsum _INCC_(i)/2+ddsum1 _INCC_(i)+ddsum2 _INCC_(i)+ddsum3 _INCC_(i)/2)/3
cc _INCC_(i)=(cc _INCC_(i)+dt*ppsum _INCC_(i))/(_ONE_+dt*ddsum _INCC_(i)/cc1 _INCC_(i))
end do
_ODE_LOOP_END_

Expand Down Expand Up @@ -554,14 +557,14 @@ subroutine modified_patankar()
do i=1,numc
a(i,i)=_ZERO_
do j=1,numc
a(i,i)=a(i,i)+dd _INPP_(i,j,ci)
if (i.ne.j) a(i,j)=-dt*pp _INPP_(i,j,ci)/cc _INCC_(j,ci)
a(i,i)=a(i,i)+dd _INPP_(i,j)
if (i.ne.j) a(i,j)=-dt*pp _INPP_(i,j)/cc _INCC_(j)
end do
a(i,i)=dt*a(i,i)/cc _INCC_(i,ci)
a(i,i)=dt*a(i,i)/cc _INCC_(i)
a(i,i)=_ONE_+a(i,i)
r(i)=cc _INCC_(i,ci)+dt*pp _INPP_(i,i,ci)
r(i)=cc _INCC_(i)+dt*pp _INPP_(i,i)
end do
call matrix(numc,a,r,cc _INCC_(:,ci))
call matrix(numc,a,r,cc _INCC_(:))
_ODE_LOOP_END_

end subroutine modified_patankar
Expand Down Expand Up @@ -633,14 +636,14 @@ subroutine modified_patankar_2()
do i=1,numc
a(i,i)=_ZERO_
do j=1,numc
a(i,i)=a(i,i)+dd _INPP_(i,j,ci)
if (i.ne.j) a(i,j)=-dt*pp _INPP_(i,j,ci)/cc _INCC_(j,ci)
a(i,i)=a(i,i)+dd _INPP_(i,j)
if (i.ne.j) a(i,j)=-dt*pp _INPP_(i,j)/cc _INCC_(j)
end do
a(i,i)=dt*a(i,i)/cc _INCC_(i,ci)
a(i,i)=dt*a(i,i)/cc _INCC_(i)
a(i,i)=_ONE_+a(i,i)
r(i)=cc _INCC_(i,ci)+dt*pp _INPP_(i,i,ci)
r(i)=cc _INCC_(i)+dt*pp _INPP_(i,i)
end do
call matrix(numc,a,r,cc1 _INCC_(:,ci))
call matrix(numc,a,r,cc1 _INCC_(:))
_ODE_LOOP_END_

call get_ppdd(.false.,_SIZE_,cc1,pp1,dd1)
Expand All @@ -652,14 +655,14 @@ subroutine modified_patankar_2()
do i=1,numc
a(i,i)=_ZERO_
do j=1,numc
a(i,i)=a(i,i)+dd _INPP_(i,j,ci)
if (i.ne.j) a(i,j)=-dt*pp _INPP_(i,j,ci)/cc1 _INCC_(j,ci)
a(i,i)=a(i,i)+dd _INPP_(i,j)
if (i.ne.j) a(i,j)=-dt*pp _INPP_(i,j)/cc1 _INCC_(j)
end do
a(i,i)=dt*a(i,i)/cc1 _INCC_(i,ci)
a(i,i)=dt*a(i,i)/cc1 _INCC_(i)
a(i,i)=_ONE_+a(i,i)
r(i)=cc _INCC_(i,ci)+dt*pp _INPP_(i,i,ci)
r(i)=cc _INCC_(i)+dt*pp _INPP_(i,i)
end do
call matrix(numc,a,r,cc _INCC_(:,ci))
call matrix(numc,a,r,cc _INCC_(:))
_ODE_LOOP_END_

end subroutine modified_patankar_2
Expand Down Expand Up @@ -698,14 +701,14 @@ subroutine modified_patankar_4()
do i=1,numc
a(i,i)=_ZERO_
do j=1,numc
a(i,i)=a(i,i)+dd _INPP_(i,j,ci)
if (i.ne.j) a(i,j)=-dt*pp _INPP_(i,j,ci)/cc _INCC_(j,ci)
a(i,i)=a(i,i)+dd _INPP_(i,j)
if (i.ne.j) a(i,j)=-dt*pp _INPP_(i,j)/cc _INCC_(j)
end do
a(i,i)=dt*a(i,i)/cc _INCC_(i,ci)
a(i,i)=dt*a(i,i)/cc _INCC_(i)
a(i,i)=_ONE_+a(i,i)
r(i)=cc _INCC_(i,ci)+dt*pp _INPP_(i,i,ci)
r(i)=cc _INCC_(i)+dt*pp _INPP_(i,i)
end do
call matrix(numc,a,r,cc1 _INCC_(:,ci))
call matrix(numc,a,r,cc1 _INCC_(:))
_ODE_LOOP_END_

call get_ppdd(first,_SIZE_,cc1,pp1,dd1)
Expand All @@ -714,14 +717,14 @@ subroutine modified_patankar_4()
do i=1,numc
a(i,i)=_ZERO_
do j=1,numc
a(i,i)=a(i,i)+dd1 _INPP_(i,j,ci)
if (i.ne.j) a(i,j)=-dt*pp1 _INPP_(i,j,ci)/cc1 _INCC_(j,ci)
a(i,i)=a(i,i)+dd1 _INPP_(i,j)
if (i.ne.j) a(i,j)=-dt*pp1 _INPP_(i,j)/cc1 _INCC_(j)
end do
a(i,i)=dt*a(i,i)/cc1 _INCC_(i,ci)
a(i,i)=dt*a(i,i)/cc1 _INCC_(i)
a(i,i)=_ONE_+a(i,i)
r(i)=cc _INCC_(i,ci)+dt*pp1 _INPP_(i,i,ci)
r(i)=cc _INCC_(i)+dt*pp1 _INPP_(i,i)
end do
call matrix(numc,a,r,cc1 _INCC_(:,ci))
call matrix(numc,a,r,cc1 _INCC_(:))
_ODE_LOOP_END_

call get_ppdd(first,_SIZE_,cc1,pp2,dd2)
Expand All @@ -730,14 +733,14 @@ subroutine modified_patankar_4()
do i=1,numc
a(i,i)=_ZERO_
do j=1,numc
a(i,i)=a(i,i)+dd2 _INPP_(i,j,ci)
if (i.ne.j) a(i,j)=-dt*pp2 _INPP_(i,j,ci)/cc1 _INCC_(j,ci)
a(i,i)=a(i,i)+dd2 _INPP_(i,j)
if (i.ne.j) a(i,j)=-dt*pp2 _INPP_(i,j)/cc1 _INCC_(j)
end do
a(i,i)=dt*a(i,i)/cc1 _INCC_(i,ci)
a(i,i)=dt*a(i,i)/cc1 _INCC_(i)
a(i,i)=_ONE_+a(i,i)
r(i)=cc _INCC_(i,ci)+dt*pp2 _INPP_(i,i,ci)
r(i)=cc _INCC_(i)+dt*pp2 _INPP_(i,i)
end do
call matrix(numc,a,r,cc1 _INCC_(:,ci))
call matrix(numc,a,r,cc1 _INCC_(:))
_ODE_LOOP_END_

call get_ppdd(first,_SIZE_,cc1,pp3,dd3)
Expand All @@ -749,14 +752,14 @@ subroutine modified_patankar_4()
do i=1,numc
a(i,i)=_ZERO_
do j=1,numc
a(i,i)=a(i,i)+dd _INPP_(i,j,ci)
if (i.ne.j) a(i,j)=-dt*pp _INPP_(i,j,ci)/cc1 _INCC_(j,ci)
a(i,i)=a(i,i)+dd _INPP_(i,j)
if (i.ne.j) a(i,j)=-dt*pp _INPP_(i,j)/cc1 _INCC_(j)
end do
a(i,i)=dt*a(i,i)/cc1 _INCC_(i,ci)
a(i,i)=dt*a(i,i)/cc1 _INCC_(i)
a(i,i)=_ONE_+a(i,i)
r(i)=cc _INCC_(i,ci)+dt*pp _INPP_(i,i,ci)
r(i)=cc _INCC_(i)+dt*pp _INPP_(i,i)
end do
call matrix(numc,a,r,cc _INCC_(:,ci))
call matrix(numc,a,r,cc _INCC_(:))
_ODE_LOOP_END_

end subroutine modified_patankar_4
Expand Down Expand Up @@ -798,8 +801,8 @@ subroutine emp_1()
call get_rhs(.true.,_SIZE_,cc,derivative)

_ODE_LOOP_BEGIN_
call findp_bisection(numc, cc _INCC_(:,ci), derivative _INCC_(:,ci), dt, 1.d-9, pi)
cc _INCC_(:,ci) = cc _INCC_(:,ci) + dt*derivative _INCC_(:,ci)*pi
call findp_bisection(numc, cc _INCC_(:), derivative _INCC_(:), dt, 1.d-9, pi)
cc _INCC_(:) = cc _INCC_(:) + dt*derivative _INCC_(:)*pi
_ODE_LOOP_END_

end subroutine emp_1
Expand Down Expand Up @@ -858,23 +861,23 @@ subroutine emp_2()
call get_rhs(.true.,_SIZE_,cc,rhs)

_ODE_LOOP_BEGIN_
call findp_bisection(numc, cc _INCC_(:,ci), rhs _INCC_(:,ci), dt, 1.d-9, pi)
cc_med _INCC_(:,ci) = cc _INCC_(:,ci) + dt*rhs _INCC_(:,ci)*pi
call findp_bisection(numc, cc _INCC_(:), rhs _INCC_(:), dt, 1.d-9, pi)
cc_med _INCC_(:) = cc _INCC_(:) + dt*rhs _INCC_(:)*pi
_ODE_LOOP_END_

call get_rhs(.false.,_SIZE_,cc_med,rhs_med)

_ODE_LOOP_BEGIN_
rhs _INCC_(:,ci) = (rhs _INCC_(:,ci) + rhs_med _INCC_(:,ci))/2
rhs _INCC_(:) = (rhs _INCC_(:) + rhs_med _INCC_(:))/2

! Correct for the state variables that will be included in 'p'.
do i=1,numc
if (rhs _INCC_(i,ci) .lt. 0.) rhs _INCC_(:,ci) = rhs _INCC_(:,ci) * cc _INCC_(i,ci)/cc_med _INCC_(i,ci)
if (rhs _INCC_(i) .lt. 0.) rhs _INCC_(:) = rhs _INCC_(:) * cc _INCC_(i)/cc_med _INCC_(i)
end do

call findp_bisection(numc, cc _INCC_(:,ci), rhs _INCC_(:,ci), dt, 1.d-9, pi)
call findp_bisection(numc, cc _INCC_(:), rhs _INCC_(:), dt, 1.d-9, pi)

cc _INCC_(:,ci) = cc _INCC_(:,ci) + dt*rhs _INCC_(:,ci)*pi
cc _INCC_(:) = cc _INCC_(:) + dt*rhs _INCC_(:)*pi
_ODE_LOOP_END_ ! ci (z-levels)

end subroutine emp_2
Expand Down