Diffusion SPClibs (feature/fft-split) Browse

fft-split

difequ.f90
View Options

			!>
			!> @file difequ.f90
			!>
			!> @brief
			!>
			!> @copyright
			!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
			!> SPC (Swiss Plasma Center)
			!>
			!> SPClibs is free software: you can redistribute it and/or modify it under
			!> the terms of the GNU Lesser General Public License as published by the Free
			!> Software Foundation, either version 3 of the License, or (at your option)
			!> any later version.
			!>
			!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
			!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
			!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
			!>
			!> You should have received a copy of the GNU Lesser General Public License
			!> along with this program. If not, see <https://www.gnu.org/licenses/>.
			!>
			!> @authors
			!> (in alphabetical order)
			!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
			!>
			subroutine difequ ( mode, xx, v )

			!*************************************************************************
			!
			!! DIFEQU returns information about a differential equation.
			!
			! Reference:
			!
			! Carl DeBoor,
			! A Practical Guide to Splines,
			! Springer Verlag.
			!
			! Parameters:
			!
			! Input, integer MODE, an integer indicating the task to be performed.
			! 1, initialization
			! 2, evaluate de at xx
			! 3, specify the next side condition
			! 4, analyze the approximation
			!
			! Input, real ( kind = 8 ) XX, a point at which information is wanted
			!
			! Output, real ( kind = 8 ) V, depends on the MODE. See comments below
			!
			implicit none

			integer, parameter :: npiece = 100
			integer, parameter :: ncoef = 2000

			real ( kind = 8 ) break
			real ( kind = 8 ) coef
			real ( kind = 8 ) eps
			real ( kind = 8 ) ep1
			real ( kind = 8 ) ep2
			real ( kind = 8 ) error
			real ( kind = 8 ) factor
			integer i
			integer iside
			integer itermx
			integer k
			integer kpm
			integer l
			integer m
			integer mode
			real ( kind = 8 ) rho
			real ( kind = 8 ) s2ovep
			real ( kind = 8 ) solutn
			real ( kind = 8 ) un
			real ( kind = 8 ) v(20)
			real ( kind = 8 ) value
			real ( kind = 8 ) x
			real ( kind = 8 ) xside
			real ( kind = 8 ) xx

			common /approx/ break(npiece),coef(ncoef),l,kpm
			common /side/ m,iside,xside(10)
			common /other/ itermx,k,rho(19)
			!
			! This sample of DIFEQU is for the example in chapter xv. It is a
			! nonlinear second order two point boundary value problem.
			!
			go to (10,50,60,110), mode
			!
			! Initialize everything, Set the order M of the differential equation,
			! the nondecreasing sequence xside(i),i=1,...,m, of points at which side
			! conditions are given and anything else necessary.
			!
			10 continue

			m = 2
			xside(1) = 0.0D+00
			xside(2) = 1.0D+00
			!
			! Print out heading.
			!
			write ( , ) ' '
			write ( , ) ' Carrier''s nonlinear perturb. problem'
			write ( , ) ' '

			eps = 0.005D+00
			write(,)'EPS = ',eps
			!
			! Set constants used in formula for solution below.
			!
			factor = ( sqrt ( 2.0D+00 ) + sqrt ( 3.0D+00 ) )**2
			s2ovep = sqrt ( 2.0D+00 / eps )
			!
			! Initial guess for Newton iteration. un(x)=x*x-1.
			!
			l = 1
			break(1) = 0.0D+00
			do i = 1, kpm
			coef(i) = 0.0D+00
			end do
			coef(1) = -1.0D+00
			coef(3) = 2.0D+00
			itermx = 10
			return
			!
			! Provide value of left side coefficients and right side at xx .
			! specifically, at xx the dif.equ. reads:
			!
			! v(m+1)dm+v(m)d(m-1) + ... + v(1)d**0 = v(m+2)
			!
			! in terms of the quantities v(i),i=1,...,m+2, to be computed here.
			!
			50 continue

			v(3) = eps
			v(2) = 0.0D+00
			call ppvalu(break,coef,l,kpm,xx,0,un)
			v(1) = 2.0D+00 * un
			v(4) = un**2 + 1.0D+00
			return
			!
			! provide the M side conditions. these conditions are of the form
			! v(m+1)dm+v(m)d(m-1) + ... + v(1)d**0 = v(m+2)
			! in terms of the quantities v(i),i=1,...,m+2, to be specified here.
			! note that v(m+1)=0 for customary side conditions.
			!
			60 continue

			v(m+1) = 0.0D+00
			if ( iside == 1 ) then
			v(2) = 1.0D+00
			v(1) = 0.0D+00
			v(4) = 0.0D+00
			iside = iside+1
			else if ( iside == 2 ) then
			v(2) = 0.0D+00
			v(1) = 1.0D+00
			v(4) = 0.0D+00
			iside = iside + 1
			end if

			return
			!
			! Calculate the error near the boundary layer at 1.
			!
			110 continue

			write(,)' '
			write(,)' X, G(X) and G(X)-F(X) at selected points:'
			write(,)' '

			x = 0.75D+00

			do i = 1, 9
			ep1 = exp ( s2ovep * ( 1.0D+00 - x ) ) * factor
			ep2 = exp ( s2ovep * ( 1.0D+00 + x ) ) * factor
			solutn = 12.0D+00 / ( 1.0D+00 + ep1 )*2 ep1 &
			+12.0D+00 / ( 1.0D+00 + ep2 )*2 ep2 - 1.0D+00
			call ppvalu(break,coef,l,kpm,x,0,value)
			error = solutn-value
			write ( *, '(1x,3g14.6)' ) x, solutn, error
			x = x+0.03125
			end do

			return
			end

spclibs/pppack/difequ.f9093a1e9712d04feature/fft-split

difequ.f90View Options

spclibs/pppack/difequ.f90
93a1e9712d04feature/fft-split

difequ.f90
View Options