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