Giter Site home page Giter Site logo

ivan-pi / pdecheb Goto Github PK

View Code? Open in Web Editor NEW
2.0 1.0 0.0 176 KB

Chebyshev Polynomial Software for Elliptic-Parabolic Systems of PDEs

Makefile 0.47% Fortran 98.62% CMake 0.91%
method-of-lines parabolic-pde partial-differential-equations elliptic-pde spatial-discretization chebyshev-polynomials fortran pde-solver

pdecheb's Introduction

---
title: An Overview of using DASSL and PDECHEB
---
flowchart BT

    subgraph L4[User provided PDE/ODE]
        l4_1[UVINIT\nPDE/ODE\nInitial\nconditions]
        l4_2[SPDEFN\nBody of PDE]
        l4_3[SBNDR\nBoundary\nconditions]
        l4_4[SODEFN\nResidual of\ncoupled\nODEs]
    end

    subgraph L3["Software to semidiscretize the PDE/ODE System"]
        l3_1[INICHB\nInitialize workspace\nand the solution vector]
        l3_2[PDECHB\nForm ODE residual]
        l3_3[INTERC\nSpatial\ninterpolation]
    end
    l4_1-->l3_1
    l4_2-->l3_2
    l4_3-->l3_2
    l4_4-->l3_2

    subgraph L2[Time integration]
        l2_1[DASSL\nIntegrate the solution\nfrom TSTART to TOUT]
    end
    l3_2-->l2_1

    L1[User's calling program]
    
    l3_1-->L1
    l3_3-->L1
    l2_1-->L1
Loading
---
title: PDECHEB Internal Architecture
---
graph BT

    subgraph A[Public interface]
        A2[INICHB]
        A1[PDECHB]
        A3[INTERC]
    end

    B1[CHINTR]-->A1
    B2[CRES and DRES]-->A1

    subgraph U["User-provided routines"]
        U1[SPDEFN]
        U2[SBNDR]
        U3[SODEFN]
        U4[UVINIT]
    end

    U1-->B2
    U2-->B2

    U1-->B1
    U3-->A1

    C1[INTRCH]-->A3
    U4-->C2[CSET]-->A2
Loading

pdecheb's People

Contributors

ivan-pi avatar

Stargazers

 avatar

Watchers

 avatar

pdecheb's Issues

Aliasing in UVINIT

The V argument is passed as an array of size 1.

pdecheb/pdecheb/cset.f

Lines 1 to 11 in f429435

SUBROUTINE CSET(NPDE,NPTS,U,X,OMEGA,DU,XBK,NEL,NPTL,XC,CCR,XBH,
* IBK,DUTEM,V,NV)
C***********************************************************************
C FORTRAN FUNCTIONS USED: SIN COS .
C***********************************************************************
C .. Scalar Arguments ..
INTEGER IBK, NEL, NPDE, NPTL, NPTS, NV
C .. Array Arguments ..
DOUBLE PRECISIONCCR(NPTL), DU(NPTL,NPTL), DUTEM(NPTL,NPTL),
* OMEGA(NPTL,NPTL), U(NPDE,NPTS), V(1), X(NPTS),
* XBH(IBK), XBK(IBK), XC(NPTL)

...

CALL UVINIT(NPDE,NPTS,X,U,NV,V)

This can have some surprising effects, if NV = 0, but the solution array would be declared as Y(NPDE*NPTS + 1) in the calling program. Here's the initialization routine that exposes the issue:

SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV, V)
IMPLICIT NONE
INTEGER :: NPDE, NPTS, NV
DOUBLE PRECISION :: X(NPTS), U(NPDE,NPTS), TIME, V(NV)
U(1,:) = 0
U(1,NPTS) = 1
V(1) = 0       ! <-- eliminates effect of previous line
END SUBROUTINE

In practice this should be flagged as an out-of-bounds violation (or produce a segfault).

Examples including integration

For completeness, we need some examples that include

  • integration of the PDE variables over the entire domain or individual elements/material domains, for $m = 0, 1, 2$
  • calculation of cumulative integrals, at the knots of the Chebyshev basis (see this question: https://scicomp.stackexchange.com/q/23398/37438)

The coefficients in the Chebyshev expansion can be obtained from the work array wkres as follows (TOMS Algorithm 660, pg. 195):

C     
      NPTL = NPOLY + 1
      DO 100 I = 1,NEL
C      ITH ELEMENT
C      IU IS THE COMPONENT OF SOLUTION VECTOR AT LHS
C      OF ELEMENT I
       IU = (I - 1)*(NPOLY)*NPDE + 1
       DO 80 IS = 1,NPTL
        DO 80 JK = 1,NPDE
         COEFF(JK,IS,I) = 0.0D0
   80  CONTINUE
       DO 90 IS = 1,NPTL
        DO 90 JS = 1,NPTL
         DO 90 JK = 2,NPDE
         COEFF(JK,IS,I) = COEFF(JK,IS,I) +
     1                    WKRES(IS + (JS-1)*NPTL)*Y(IU + JS + JK - 2)
   90  CONTINUE
  100 CONTINUE
C

Or using free-form Fortran:

nptl = npoly + 1
do i = 1, nel
   ! ith element
   ! iu is the component of solution vector at LHS of element i
   iu = (i - 1)*npoly*npde + 1
   do is = 1, nptl
      do jk = 1, npde
         coeff(jk,is,i) = 0.0d0
      end do
   end do
   do is = 1, nptl
      do js = 1, nptl
         do jk = 2, npde
            coeff(jk,is,i) = coeff(jk,is,i) + &
               wkres(is + (js - 1)*nptl) * y(iu + js + jk - 2)
         end do
      end do
   end do
end do

The coefficient $a_{j,n}(t)$ of the $k$-th PDE is stored in COEFF(K,J,N), where $j$ runs along elements in the mesh, and $n$ counts the Chebyshev basis polynomials $T_n(x)$.

More test problems

We should add more test problems to exercise the PDECHEB code in as many modes as possible.

A few relevant codes in this area are:

A few more solvers are given in the GAMS Class I2a1a.

Also the following works (and works citing these) may be worth looking into:

Lots of parabolic-elliptic initial-value problems can be found in textbooks related to chemical engineering, electrochemistry, transport phenomena, etc.

DAE Solvers

A few improvements could be made to DASSL:

  • replace LINPACK routines with their LAPACK counterparts
  • link against an optimized BLAS library

It would also be of interest to provide a PDECHB callback compatible with the newer DASPK solver. DASPK can be found on Netlib: https://netlib.org/ode/. The method is described in

Brown, P. N., Hindmarsh, A. C., & Petzold, L. R. (1994). Using Krylov methods in the solution of large-scale differential-algebraic systems. SIAM Journal on Scientific Computing, 15(6), 1467-1488. https://doi.org/10.1137/0915088

A few more DAE solvers can be found in the Netlib folder (e.g. daesolve). Other notable ones include:

Constant for machine precision

The code below

pdecheb/pdecheb/inichb.f

Lines 268 to 279 in f429435

C
C CALCULATE ROUGH ESTIMATE OF UNIT ROUND-OFF ERROR FOR CHECKING
C
TWOU = 0.1D0
40 TEMP = 1.0D0 + TWOU
IF (1.0D0.EQ.TEMP) THEN
TWOU = TWOU*2.0D0
ELSE
TWOU = TWOU*0.5D0
GO TO 40
END IF
C

can be replaced with a call to the intrinsic function epsilon:

twou = epsilon(1.0D0)

This could also become a parameter instead of a dynamic variable in a common block.

Out of bounds access in CHINTR

See run: https://github.com/ivan-pi/pdecheb/actions/runs/4349430829

When runtime checks are turned on with -fcheck=all, we get an out of bounds access at the following line:

IF (XP(IP).LT.(XBK(I)*TEM1-TWOU)) GO TO 20

At line 95 of file /home/runner/work/pdecheb/pdecheb/pdecheb/chintr.f
 TEST PROBLEM 1
 ***********
 POLY OF DEGREE =   2 NO OF ELEMENTS =   20
Fortran runtime error: Index '0' of dimension 1 of array 'xp' below lower bound of 1

Error termination. Backtrace:
#0  0x7faf8482dad0 in ???
#1  0x7faf8482e649 in ???
#2  0x7faf8482ec46 in ???
#3  0x55d0fb8a71ae in chintr_
	at /home/runner/work/pdecheb/pdecheb/pdecheb/chintr.f:95
#4  0x55d0fb8a5297 in pdechb_
	at /home/runner/work/pdecheb/pdecheb/pdecheb/pdechb.f:67
#5  0x55d0fb8adc18 in ddastp_
	at /home/runner/work/pdecheb/pdecheb/ddassl.f:2417
#6  0x55d0fb8b1a67 in ddassl_
	at /home/runner/work/pdecheb/pdecheb/ddassl.f:[13](https://github.com/ivan-pi/pdecheb/actions/runs/4349430829/jobs/7599092007#step:7:14)[22](https://github.com/ivan-pi/pdecheb/actions/runs/4349430829/jobs/7599092007#step:7:23)
#7  0x55d0fb89f676 in MAIN__
	at /home/runner/work/pdecheb/pdecheb/examples/example1.f:55
#8  0x55d0fb89faba in main
	at /home/runner/work/pdecheb/pdecheb/examples/example1.f:78
Error: Process completed with exit code 2.

Replace LINPACK with LAPACK in DASSL

The DDASSL routine uses the following double precision LINPACK routines:

  • DGEFA/DGESL
  • DGBFA/DGBSL

with dependencies on the following BLAS routines:

  • DSCAL
  • DAXPY
  • DDOT
  • IDAMAX

Ideally, we should replace the LINPACK routines with their equivalent LAPACK replacements: DGETRF/DGETRS (general) and DGBTRF/DGBTRS (banded).

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.