1!
2!     CalculiX - A 3-dimensional finite element program
3!     Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19!     diagonal preconditioning of a matrix
20!
21      subroutine inviscidpre(nk,inomat,ntmat_,shcon,nshcon,physcon,
22     &     xmach2,nactdoh,b,vold,v,mi)
23!
24      implicit none
25!
26      integer i,nk,imat,inomat(*),ntmat_,nshcon(*),nactdoh(0:4,*),
27     &     idof,mi(*)
28!
29      real*8 temp,shcon(0:3,ntmat_,*),cp,physcon(*),r,xkappa,vel2,
30     &     xmach,xmach2(*),constant,b(*),vold(0:mi(2),*),
31     &     v(0:mi(2),*)
32!
33      do i=1,nk
34        imat=inomat(i)
35        if(imat.eq.0) cycle
36!
37        idof=nactdoh(0,i)
38        temp=vold(0,i)
39!
40!     material properties r, cp and kappa
41!
42        call materialdata_cp_sec(imat,ntmat_,temp,shcon,
43     &       nshcon,cp,physcon)
44        r=shcon(3,1,imat)
45        xkappa=cp/(cp-r)
46!
47!     size of the velocity
48!
49        vel2=vold(1,i)**2+vold(2,i)**2+vold(3,i)**2
50!
51        xmach=dsqrt(vel2)/(xkappa*r*(temp-physcon(1)))
52!
53!     determine the relevant Mach number
54!
55        if(xmach.lt.1.d-5) then
56          xmach2(idof)=1.d-10
57        elseif(xmach.lt.1.d0) then
58          xmach2(idof)=xmach**2
59        else
60          xmach2(idof)=1.d0
61        endif
62!
63        constant=1.d0-1.d0/xmach2(idof)
64!
65!     modifying the right hand side
66!
67        b(idof)=b(idof)+vel2*constant/3.d0*v(0,i)-constant*
68     &       -(vold(1,i)*v(1,i)+vold(2,i)*v(2,i)+vold(3,i)*v(3,i))
69!
70      enddo
71!
72      return
73      end
74