femwind/fortran/module_ndt_w_assembly.f90

   1 module module_w_assembly
   2        use module_hexa
   3        use module_lin_alg
   4
   5 contains
   6
   7 subroutine w_assembly(                            &
   8     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
   9     ifms, ifme, kfms,kfme, jfms, jfme,            &
  10     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  11     ifts, ifte, kfts, kfte, jfts,jfte,            &
  12     lambda,u0, v0, w0, A,  X, Y, Z,               & !Input from femwind, u0, v0, w0, Spatial Grid Data
  13     U,V,W)                                          !U,V,W
  14 !Purpose: Create Arrays of Wind Vector Component Values at Center Points of Spatial Grid
  15 !In:
  16 !A Coefficient Matrix size 3X3, symmetric positive definite
  17 !u0, v0, w0  Initial wind speed values in x,y,z direction at grid cell centers
  18 !iflags1 iflags1 = 1 returns Kloc and Jg from hexa, iflags2 = 2 returns Floc and Jg from hexa
  19 !iflags2 iflags2 =1  indicates add initial wind to calculated wind
  20 !out:
  21 ! U,V,W Computed wind values in x,y,z direction
  22
  23 implicit none
  24
  25
  26 !*** arguments
  27
  28 integer, intent(in)::                             &
  29     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
  30     ifms, ifme, kfms,kfme, jfms, jfme,            &
  31     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  32     ifts, ifte, kfts, kfte, jfts,jfte
  33
  34
  35
  36
  37
  38 real, intent(in), dimension(ifts:ifte+1, kfts:kfte+1, jfts: jfte+1):: lambda
  39 !Input for hexa
  40
  41
  42 real, intent(in) :: A(3,3)
  43 real, intent(in), dimension(ifms:ifme, kfms:kfme, jfms:jfme)::X,Y,Z,u0, v0, w0
  44
  45 real, intent(out), dimension(ifms:ifme, kfms:kfme, jfms:jfme)::U,V,W
  46 !*** local
  47
  48 integer:: ie1, ie2, ie3, ic1, ic2, ic3, iloc, i, &
  49           k1, k2, k3
  50 real ::  Kloc(8,8), Floc(8), Jg(8,3)
  51 real ::  Xloc(3,8), u0loc(3)
  52 real :: grad(3)
  53 real :: lambda_loc(8)
  54 real :: A_inv(3,3)
  55 lambda_loc = 0.
  56 Xloc = 99999.
  57 Jg = 0.
  58 Kloc = 0.
  59 Floc = 0.
  60 grad = 0.
  61 u0loc =0.
  62
  63
  64 !*** u0loc is an input for module_hexa, but is not used to construct K. Do I need to define this?
  65 !** executable
  66
  67 !print *, 'u0 vector is', u0
  68 do ie2=jfts,jfte
  69     do ie3=kfts, kfte
  70         do ie1=ifts, ifte
  71             do ic2=0,1
  72                 do ic3=0,1
  73                     do ic1=0,1
  74                         iloc=1+ic1+2*(ic2+2*ic3);  !local index of the node in the element
  75                         Xloc(1,iloc)=X(ie1 + ic1, ie3 + ic3, ie2 + ic2)
  76                         Xloc(2,iloc)=Y(ie1 + ic1, ie3 + ic3, ie2 + ic2)
  77                         Xloc(3,iloc)=Z(ie1 + ic1, ie3 + ic3, ie2 + ic2)
  78                         k1 = ie1 + ic1
  79                         k2 = ie2 + ic2
  80                         k3 = ie3 + ic3
  81                         lambda_loc(iloc) = lambda(k1,k3,k2)
  82                     enddo
  83                 enddo
  84             enddo
  85             u0loc(1) = u0(ie1,ie3,ie2)
  86             u0loc(2) = v0(ie1, ie3,ie2)
  87             u0loc(3) = w0(ie1, ie3,ie2)
  88             !fine print *, 'local lambda is', lambda_loc
  89             !print* , 'Xloc is', Xloc
  90             !print* , 'u0loc is', u0loc(1)
  91             call hexa(A,Xloc,u0loc,Kloc,Floc,Jg,3)
  92             !print*, 'Jg is', Jg
  93             !print*, shape(jg)
  94             grad = matmul(transpose(Jg),lambda_loc)
  95             !not ok print *,'Grad before multiplication by A_inv is', grad
  96
  97             call Inv3(A, A_inv)
  98             grad = matmul(transpose(A_inv),grad)
  99             ! Not ok print *,'Grad after multiplication by A_inv is', grad
 100
 101             U(ie1, ie3, ie2)=grad(1)+ u0(ie1, ie3, ie2)
 102             V(ie1, ie3, ie2)=grad(2)+ v0(ie1, ie3, ie2)
 103             W(ie1, ie3, ie2)=grad(3)+ w0(ie1, ie3, ie2)
 104
 105         enddo
 106        !print *, 'lambda array', lambda
 107     enddo
 108 enddo
 109 !print *,'Shape of U', shape(U)
 110
 111 end subroutine w_assembly
 112 end module  module_w_assembly
 113
 114
 115
 116
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132