Commit 1199140d authored by Adam Powell's avatar Adam Powell

Imported Upstream version 6.1.0.svn.5272.dfsg

parent bd325599
......@@ -227,7 +227,7 @@
<parameter name="Viscosity Model" >Power law</parameter>
<parameter name="Viscosity Exponent" >$ (1.0/3.0)</parameter>
<parameter name="Critical Shear Rate" >1.0E-6</parameter>
<parameter name="Viscosity" >Variable Temperature; Real MATC "(2.0*3.0*1.916E03 * exp( -139.0E03/(8.314 *(tx+273.15)))^(-1.0/3.0)</parameter>
<parameter name="Viscosity" >Variable Temperature; Real MATC "(2.0*3.0*1.916E03 * exp( -139.0E03/(8.314 *(tx+273.15))))^(-1.0/3.0)</parameter>
</material>
......
This diff is collapsed.
......@@ -493,7 +493,7 @@ int LoadElmerInput(struct FemType *data,struct BoundaryType *bound,
{
int noknots,noelements,nosides,maxelemtype;
int sideind[MAXNODESD1],tottypes,elementtype;
int i,j,k,dummyint,cdstat;
int i,j,k,l,dummyint,cdstat;
FILE *in;
char line[MAXLINESIZE],filename[MAXFILESIZE],directoryname[MAXFILESIZE];
......@@ -582,8 +582,13 @@ int LoadElmerInput(struct FemType *data,struct BoundaryType *bound,
bigerror("Cannot continue with invalid elements");
}
data->elementtypes[j] = elementtype;
for(k=0;k< elementtype%100 ;k++)
fscanf(in,"%d",&(data->topology[j][k]));
for(k=0;k< elementtype%100 ;k++) {
fscanf(in,"%d",&l);
data->topology[j][k] = l;
if(l < 0 || l > noknots ) {
printf("node out of range: %d %d %d %d %d\n",i,j,elementtype,k,l);
}
}
}
fclose(in);
......@@ -2033,8 +2038,8 @@ int PartitionSimpleElements(struct FemType *data,int dimpart[],int dimper[],
}
else {
cx = 1.0;
cy = 0.01;
cz = 0.0001;
cy = 0.0001;
cz = cy*cy;
}
z = 0.0;
......@@ -2269,8 +2274,8 @@ int PartitionSimpleNodes(struct FemType *data,int dimpart[],int dimper[],
}
else {
cx = 1.0;
cy = 0.01;
cz = 0.0001;
cy = 0.0001;
cz = cy*cy;
}
z = 0.0;
......@@ -2825,7 +2830,7 @@ int PartitionMetisNodes(struct FemType *data,struct BoundaryType *bound,
static void CheckPartitioning(struct FemType *data,int info)
{
int i,j,partitions,part,part2,noknots,noelements,mini,maxi,sumi,hit,ind,nodesd2,elemtype;
int i,j,k,partitions,part,part2,noknots,noelements,mini,maxi,sumi,hit,ind,nodesd2,elemtype;
int *elempart, *nodepart,*elemsinpart,*nodesinpart,*sharedinpart;
noknots = data->noknots;
......@@ -2882,6 +2887,7 @@ static void CheckPartitioning(struct FemType *data,int info)
}
if(info) {
printf("Information on partition bandwidth\n");
if(partitions <= 4) {
printf("Distribution of elements, nodes and shared nodes\n");
printf(" %-10s %-10s %-10s %-10s\n","partition","elements","nodes","shared");
......@@ -2894,14 +2900,14 @@ static void CheckPartitioning(struct FemType *data,int info)
mini = MIN( elemsinpart[i], mini);
maxi = MAX( elemsinpart[i], maxi);
}
printf("There are in average %d elements with range %d in partition\n",noelements/partitions,maxi-mini);
printf("Average %d elements with range %d in partition\n",noelements/partitions,maxi-mini);
mini = maxi = nodesinpart[1];
for(i=1;i<=partitions;i++) {
mini = MIN( nodesinpart[i], mini);
maxi = MAX( nodesinpart[i], maxi);
}
printf("There are in average %d nodes with range %d in partition\n",noknots/partitions,maxi-mini);
printf("Average %d nodes with range %d in partition\n",noknots/partitions,maxi-mini);
sumi = 0;
mini = maxi = sharedinpart[1];
......@@ -2910,7 +2916,7 @@ static void CheckPartitioning(struct FemType *data,int info)
maxi = MAX( sharedinpart[i], maxi);
sumi += sharedinpart[i];
}
printf("There are in average %d shared nodes with range %d in partition\n",sumi/partitions,maxi-mini);
printf("Average %d shared nodes with range %d in partition\n",sumi/partitions,maxi-mini);
}
}
......@@ -2919,21 +2925,25 @@ static void CheckPartitioning(struct FemType *data,int info)
if(0) printf("Checking that each node in elements belongs to nodes\n");
for(i=1;i<=data->noelements;i++) {
part = elempart[i];
elemtype = data->elementtypes[j];
nodesd2 = elemtype%100;
elemtype = data->elementtypes[i];
nodesd2 = elemtype % 100;
for(j=0;j < nodesd2;j++)
for(j=0;j < nodesd2;j++) {
ind = data->topology[i][j];
hit = FALSE;
for(j=1;j<=data->maxpartitiontable;j++) {
part2 = data->partitiontable[j][ind];
if( part == part2 ) hit = TRUE;
if(hit && !part) break;
}
if(!hit) {
printf("******** Warning *******\n");
printf("Node %d in element %d does not belong to partition %d (%d)\n",ind,i,part,j);
hit = FALSE;
for(k=1;k<=data->maxpartitiontable;k++) {
part2 = data->partitiontable[k][ind];
if( part == part2 ) hit = TRUE;
if(hit && !part) break;
}
if(!hit) {
printf("******** Warning *******\n");
printf("Node %d in element %d does not belong to partition %d (%d)\n",ind,i,part,part2);
printf("elemtype = %d nodesd2 = %d\n",elemtype,nodesd2);
for(k=0;k < nodesd2;k++)
printf("ind[%d] = %d\n",k,data->topology[i][k]);
}
}
}
......@@ -3196,6 +3206,7 @@ static int RenumberCuthillMckee( int nrows, int *rows, int *cols, int *iperm )
for(i=0; i<nrows; i++ )
for( j=rows[i]; j<rows[i+1]; j++ )
bw_aft = MAX( bw_aft, ABS(iperm[cols[j]]-iperm[i])+1 );
printf( "RCM: Bandwidth after: %d\n", bw_aft );
free_Ivector(level,0,nrows-1);
......@@ -3288,8 +3299,6 @@ static void RenumberPartitions(struct FemType *data,int info)
}
xadj[partitions] = totcon;
perm = Ivector(0,partitions-1);
bw_reduced = RenumberCuthillMckee( partitions, xadj, adjncy, perm );
......@@ -3390,30 +3399,33 @@ int OptimizePartitioning(struct FemType *data,struct BoundaryType *bound,int noo
if(partbw) RenumberPartitions(data,info);
/* Check partitioning after table is created for the first time */
printf("Checking partitioning before optimization\n");
CheckPartitioning(data,info);
/* Distribute the shared nodes as evenly as possible.
These store the load balancing information. */
/* Calculate how many nodes is owned by each partition */
neededvector = Ivector(1,partitions);
for(i=1;i<=partitions;i++)
neededvector[i] = 0;
/* Make the initial distribution that points the ownerships. */
for(i=1;i<=noknots;i++)
neededvector[nodepart[i]] += 1;
optimize = 1;
probnodes = Ivector(1,noknots);
for(i=1;i<=noknots;i++)
probnodes[i] = 0;
if(!noopt) printf("Applying aggressive optimization for load balancing\n");
if(!noopt) {
optimize = 1;
probnodes = Ivector(1,noknots);
for(i=1;i<=noknots;i++)
probnodes[i] = 0;
printf("Applying aggressive optimization for load balancing\n");
}
optimizeownership:
dshared = CheckSharedDeviation(neededvector,partitions,info);
if(!noopt) {
/* Distribute the shared nodes as evenly as possible. */
int target, same, nochanges, maxrounds, dtarget;
target = noknots / partitions;
......@@ -3596,7 +3608,7 @@ int OptimizePartitioning(struct FemType *data,struct BoundaryType *bound,int noo
ind = data->topology[i][j];
k = nodepart[ind];
if((k == i1 && e1 < e2) || (k == i2 && e1 >= e2)) {
probnodes[ind] += 1;
if(!noopt) probnodes[ind] += 1;
nodepart[ind] = elempart[i];
neededvector[elempart[i]] += 1;
neededvector[k] -= 1;
......@@ -3618,10 +3630,6 @@ int OptimizePartitioning(struct FemType *data,struct BoundaryType *bound,int noo
}
}
printf("A posteriori checking\n");
CheckPartitioning(data,info);
/* This seems to work also iteratively */
if(!noopt && m+n > 10 && optimize < 50) {
optimize++;
......@@ -3630,9 +3638,15 @@ int OptimizePartitioning(struct FemType *data,struct BoundaryType *bound,int noo
}
free_Ivector(neededvector,1,partitions);
free_Ivector(probnodes,1,noknots);
if(!noopt) free_Ivector(probnodes,1,noknots);
if(info) printf("The partitioning was optimized.\n");
printf("Checking partitioning after optimization\n");
CheckPartitioning(data,info);
return(0);
}
......
......@@ -245,7 +245,7 @@ void InitParameters(struct ElmergridType *eg)
eg->nodes3d = 0;
eg->metis = 0;
eg->partopt = 0;
eg->partoptim = 0;
eg->partoptim = FALSE;
eg->partjoin = 0;
eg->partitionhalo = FALSE;
eg->partitionindirect = FALSE;
......
......@@ -2133,8 +2133,6 @@ int SetDiscontinuousBoundary(struct FemType *data,struct BoundaryType *bound,
int mat1,mat2,par1,par2,mat1old,mat2old,material;
static int hitsexist=FALSE,hitslength,*hits;
printf("a1\n");
if(boundtype < 0) {
newbc = TRUE;
......@@ -3241,7 +3239,6 @@ int CloneMeshes(struct FemType *data,struct BoundaryType *bound,
newy = Rvector(1,noknots);
if(data->dim == 3) newz = Rvector(1,noknots);
printf("a1\n");
for(l=0;l<ncopies[2];l++) {
for(k=0;k<ncopies[1];k++) {
......
......@@ -504,7 +504,6 @@ int main(int argc, char *argv[])
case 19:
printf("a1\n");
boundaries[nofile] = (struct BoundaryType*)
malloc((size_t) (MAXBOUNDARIES)*sizeof(struct BoundaryType));
......@@ -877,6 +876,7 @@ int main(int argc, char *argv[])
printf( "------------------------------\n");
timer_show();
if(eg.periodicdim[0] || eg.periodicdim[1] || eg.periodicdim[2])
FindPeriodicNodes(&data[k],eg.periodicdim,info);
......
......@@ -1366,6 +1366,100 @@ END FUNCTION CRS_RowSum
!------------------------------------------------------------------------------
SUBROUTINE CRS_BlockMatrixPick(A,B,Blocks,Nrow,Ncol)
!------------------------------------------------------------------------------
!******************************************************************************
!
! DESCRIPTION:
! Pics a block from matrix A to build matrix B.
!------------------------------------------------------------------------------
TYPE(Matrix_t) :: A, B
INTEGER :: Blocks, Nrow,Ncol
INTEGER :: i,j,k,l,kb,n,Nrow0,Ncol0,nsub
INTEGER :: lsub,isub,istat,modNcol
LOGICAL :: NewMatrix, Diagonal
IF(Blocks <= 1) THEN
CALL Fatal('CRS_BlockMatrixPick','No applicable to just one block!')
RETURN
END IF
N = A % NumberOfRows
Nsub = N / Blocks
modNcol = MOD( Ncol,Blocks)
NewMatrix = ( B % NumberOfRows == 0 )
Diagonal = ( Nrow == Ncol )
IF( NewMatrix ) THEN
B % ListMatrix => NULL()
B % FORMAT = MATRIX_CRS
B % NumberOfRows = Nsub
kb = 0
DO isub=1,Nsub
i = Blocks * ( isub - 1 ) + Nrow
DO k= A % Rows(i), A % Rows(i+1)-1
l = A % Cols(k)
IF( MOD(l,Blocks) == modNcol ) THEN
kb = kb + 1
END IF
END DO
END DO
IF( kb == 0 ) THEN
PRINT *,'Nrow:',Nrow,'Ncol:',Ncol
CALL Warn('CRS_BlockMatrixPick','No matrix entries in submatrix')
RETURN
END IF
ALLOCATE(B % Rows(nsub+1),B % Cols(kb), B % Values(kb),STAT=istat )
IF( istat /= 0 ) CALL Fatal('CRS_BlockMatrixPick','memory allocation error 1')
IF( Diagonal ) THEN
ALLOCATE( B % Diag(nsub), B % rhs(nsub), STAT=istat)
IF( istat /= 0 ) CALL Fatal('CRS_BlockMatrixPick','memory allocation error 2')
END IF
END IF
kb = 1
DO isub=1,Nsub
IF( NewMatrix ) B % Rows(isub) = kb
i = Blocks * ( isub - 1 ) + Nrow
DO k = A % Rows(i), A % Rows(i+1)-1
l = A % Cols(k)
IF( MOD( l, Blocks ) == modNcol ) THEN
lsub = ( l - 1) / Blocks + 1
B % Values(kb) = A % Values(k)
IF( NewMatrix ) THEN
B % Cols(kb) = lsub
IF( Diagonal .AND. isub == lsub ) B % Diag(isub) = kb
END IF
kb = kb + 1
END IF
END DO
IF( Diagonal ) B % rhs(isub) = A % rhs(i)
END DO
IF( NewMatrix ) B % Rows(Nsub+1) = kb
!------------------------------------------------------------------------------
END SUBROUTINE CRS_BlockMatrixPick
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
FUNCTION CRS_IncompleteLU(A,ILUn) RESULT(Status)
!------------------------------------------------------------------------------
......@@ -1414,7 +1508,6 @@ END FUNCTION CRS_RowSum
CALL Info( 'CRS_IncompleteLU', Message, Level = 5 )
st = CPUTime()
N = A % NumberOfRows
Diag => A % Diag
Rows => A % Rows
......@@ -1468,6 +1561,7 @@ END FUNCTION CRS_RowSum
C = .FALSE.
S = 0.0d0
IF ( A % Cholesky ) THEN
ALLOCATE( T(n) )
......@@ -1534,6 +1628,7 @@ END FUNCTION CRS_RowSum
ELSE
! The factorization row by row:
! -----------------------------
DO i=1,N
......@@ -1590,6 +1685,7 @@ END FUNCTION CRS_RowSum
END DO
END IF
!------------------------------------------------------------------------------
WRITE(Message,'(a,i1,a,i9)') 'ILU(', ILUn, &
') (Real), NOF nonzeros: ',ILURows(n+1)
......
......@@ -809,10 +809,13 @@ CONTAINS
n_nbr,comm,i_type,d_type,me,mactive,maxnz,proc
logical :: iterative, found
integer,allocatable :: l_nz(:),g_nz(:),l_i(:),g_i(:),l_nbr(:)
real(kind=dp), allocatable :: l_g(:,:),l_gtg(:,:),work(:),g_gtg(:,:)
real(kind=dp), allocatable :: l_g(:,:),l_gtg(:,:),g_x(:),g_b(:),g_gtg(:,:)
save :: g_nz,g_i,g_gtg
type(matrix_t), pointer :: gtg_m => null()
save :: g_nz,g_i,gtg_m,g_gtg
!------------------------------------------------------------------------------
!call resettimer('project')
maxnz = parallelreduction(1._dp*nz,2)
......@@ -839,6 +842,7 @@ CONTAINS
d_type = mpi_double_precision
if (op==1) then
!call resettimer('factorize')
! If not using CG, assemble and factorize the G'G matrix:
!
! NOTE for further development: Instead of using one task to do the
......@@ -952,8 +956,11 @@ CONTAINS
! Get storage for the (global) G'G:
! ---------------------------------
if(allocated(g_gtg)) deallocate(g_gtg)
allocate(g_gtg(g_n,g_n)); g_gtg=0;
! if(allocated(g_gtg)) deallocate(g_gtg)
! allocate(g_gtg(g_n,g_n)); g_gtg=0;
IF (ASSOCIATED(gtg_m)) CALL FreeMatrix(gtg_m)
gtg_m=>AllocateMatrix()
gtg_m % format = MATRIX_LIST
! Assemble our own part...
! ------------------------
......@@ -962,7 +969,7 @@ CONTAINS
n_nbr=nneigh
l_nbr(0)=me; l_nbr(1:n_nbr)=gpnum(1:n_nbr)
call addtogtg(g_gtg,l_gtg,n_nbr,l_nbr,g_nz,g_i)
call addtogtg(gtg_m,l_gtg,n_nbr,l_nbr,g_nz,g_i)
deallocate(l_gtg,l_nbr)
! ...and the rest of it:
......@@ -979,7 +986,7 @@ CONTAINS
allocate(l_gtg(l_n,l_n))
call mpi_recv(l_gtg,l_n**2,d_type,proc,404,comm,stat,ierr)
call addtogtg(g_gtg,l_gtg,n_nbr,l_nbr,g_nz,g_i)
call addtogtg(gtg_m,l_gtg,n_nbr,l_nbr,g_nz,g_i)
deallocate(l_gtg,l_nbr)
end do
......@@ -988,9 +995,11 @@ CONTAINS
! to Lapack full matrix Cholesky factorization, at least for
! the time being...):
! ------------------------------------------------------------
call dpotrf('u',g_n,g_gtg,g_n,ierr)
call list_tocrsmatrix(gtg_m)
! call dpotrf('u',g_n,g_gtg,g_n,ierr)
end if ! task 'mactive' doing the factoring (and later, solving)
!call checktimer('factorize',delete=.true.)
end if ! the first time G'G assembly and factorizing
end if ! direct solver setup
......@@ -1011,33 +1020,36 @@ CONTAINS
! assemble the r.h.s. vector:
! ---------------------------
g_n=size(g_gtg,1)
allocate(work(g_n)); work=0
! g_n=size(g_gtg,1)
g_n=gtg_m % numberofrows
allocate(g_x(g_n),g_b(g_n)); g_x=0; g_b=0
if (nz>0) work(1:nz) = b
if (nz>0) g_b(1:nz) = b
do i=0,parenv % pes-1
if (i==me .or. .not. parenv % active(i+1)) cycle
if (g_nz(i)>0) &
call mpi_recv(work(g_i(i):g_i(i+1)-1), g_nz(i), &
call mpi_recv(g_b(g_i(i):g_i(i+1)-1), g_nz(i), &
d_type,i,405,comm,stat,ierr)
end do
! solve x (lapack routine used here):
! -----------------------------------
call dpotrs('u',g_n,1,g_gtg,g_n,work,g_n,ierr)
! g_x=g_b
! call dpotrs('u',g_n,1,g_gtg,g_n,g_x,g_n,ierr)
call directsolver(gtg_m,g_x,g_b,getsolver())
! distribute the result:
! ----------------------
do i=0,parenv % pes-1
if (i==me .or. .not. parenv % active(i+1)) cycle
if (g_nz(i)>0) &
call mpi_bsend(work(g_i(i):g_i(i+1)-1), &
call mpi_bsend(g_x(g_i(i):g_i(i+1)-1), &
g_nz(i),d_type,i,406,comm,ierr)
end do
! finally my own part:
! --------------------
if (nz>0) x(1:nz)=work(1:nz)
if (nz>0) x(1:nz)=g_x(1:nz)
else
! send r.h.s.:
! ------------
......@@ -1062,13 +1074,17 @@ CONTAINS
T=T-S
END IF
END IF
!call checktimer('project',delete=.true.)
CONTAINS
!------------------------------------------------------------------------------
subroutine addtogtg(g_gtg,l_gtg,n_nbr,l_nbr,g_nz,g_i)
subroutine addtogtg(gtg_m,l_gtg,n_nbr,l_nbr,g_nz,g_i)
! subroutine addtogtg(g_gtg,l_gtg,n_nbr,l_nbr,g_nz,g_i)
!------------------------------------------------------------------------------
real(kind=dp) :: g_gtg(:,:),l_gtg(:,:)
real(kind=dp) :: l_gtg(:,:)
type(matrix_t), pointer :: gtg_m
! real(kind=dp) :: g_gtg(:,:)
integer :: n_nbr, l_nbr(0:), g_nz(0:), g_i(0:)
!------------------------------------------------------------------------------
integer :: i,j,k,l,g_r,g_c,l_r,l_c,l_i(0:n_nbr+1)
......@@ -1086,7 +1102,9 @@ CONTAINS
do l=0,g_nz(l_nbr(j))-1
l_r=l_i(j)+l
g_r=g_i(l_nbr(j))+l
g_gtg(g_r,g_c) = g_gtg(g_r,g_c) + l_gtg(l_r,l_c)
! g_gtg(g_r,g_c)=g_gtg(g_r,g_c)+l_gtg(l_r,l_c)
if (abs(l_gtg(l_r,l_c))>aeps) &
call addtomatrixelement(gtg_m,g_r,g_c,l_gtg(l_r,l_c))
end do
end do
end do
......@@ -1204,6 +1222,7 @@ END SUBROUTINE FetiProject
LOGICAL :: Found,prec
REAL(KIND=dp), ALLOCATABLE :: Ri(:),S(:),T(:),P(:),Ssave(:,:),Psave(:,:),srho(:)
!------------------------------------------------------------------------------
!call resettimer('cpg')
nrows=A % NumberOfRows
Params => GetSolverParams()
......@@ -1247,6 +1266,7 @@ END SUBROUTINE FetiProject
beta=0._dp
DO iter=1,maxit
!call resettimer('iter')
IF (Restart>0) P=T
IF (Precondition) THEN
......@@ -1295,12 +1315,14 @@ END SUBROUTINE FetiProject
END IF
IF (err1<TOL) EXIT
!call checktimer('iter',delete=.true.)
END DO
CALL matvecsubr(x(1:n),T,ipar)
T = -(T-b(1:n))
CALL FetiProject(A,n,T,OP=2,TOL=1d-12)
IF (nz>0) x(n+1:n+nz) = T(1:nz)
!call checktimer('cpg',delete=.true.)
!------------------------------------------------------------------------------
END SUBROUTINE FetiCPG
!------------------------------------------------------------------------------
......@@ -1318,6 +1340,7 @@ END SUBROUTINE FetiProject
INTEGER :: FixInds(:)
TYPE(Solver_t) :: Solver
!------------------------------------------------------------------------------
REAL(KIND=dp), PARAMETER :: floatEps = 1.0d-12
REAL(KIND=dp), ALLOCATABLE :: x(:),tz(:,:)
INTEGER, POINTER :: p(:)
INTEGER :: i,j,k,n,m,dofs,neigs,dim,FixNodes(6)
......@@ -1368,7 +1391,7 @@ END SUBROUTINE FetiProject
ALLOCATE(x(n))
CALL MatrixVectorMultiply(A,z(1,:),x)
Floating = ALL(x<100*AEPS)
Floating = ALL(x<floatEps)
IF (.NOT.Floating) THEN
DEALLOCATE(z); nz=0;
......@@ -1429,7 +1452,7 @@ END SUBROUTINE FetiProject
nz=0
DO i=1,j
CALL MatrixVectorMultiply(A,tz(i,:),x)
IF (ALL(ABS(x)<100*AEPS)) THEN
IF (ALL(ABS(x)<floatEps)) THEN
nz=nz+1
FloatInds(nz)=i
END IF
......@@ -1501,7 +1524,7 @@ END SUBROUTINE FetiProject
! Finally create null(A) from zero freq. eigenvectors:
! ----------------------------------------------------
DO nz=0,neigs-1
IF (ABS(EigValues(nz+1))>100*AEPS) EXIT
IF (ABS(EigValues(nz+1))>floatEps) EXIT
END DO
IF (nz>0) THEN
......@@ -1530,7 +1553,7 @@ END SUBROUTINE FetiProject
!------------------------------------------------------------------------------
INTEGER :: n
REAL(KIND=dp), POINTER CONTIG :: tx(:),tb(:)
!call resettimer('direct')
n = A % NumberOfRows
tx=>x
tb=>b
......@@ -1549,6 +1572,7 @@ END SUBROUTINE FetiProject
x=tx(1:n)
DEALLOCATE(tx,tb)
END IF
!call checktimer('direct')
!------------------------------------------------------------------------------
END SUBROUTINE FetiDirectSolver
!------------------------------------------------------------------------------
......@@ -1766,6 +1790,7 @@ END SUBROUTINE FetiProject
TYPE(Matrix_t), POINTER :: A
TYPE(Solver_t), POINTER :: Solver
!call resettimer('mv')
Solver => GetSolver()
A => GetMatrix()
n = A % NumberOfRows
......@@ -1786,6 +1811,7 @@ END SUBROUTINE FetiProject
v(1:nLC) = v(1:nLC) + w(1:nLC)
IF (nz>0) v(nLC+1:nLC+nz) = MATMUL(z,b)
END IF
!call checktimer('mv',delete=.true.)
!------------------------------------------------------------------------------
END SUBROUTINE FetiMV
!------------------------------------------------------------------------------
......@@ -1808,6 +1834,7 @@ END SUBROUTINE FetiProject
INTEGER :: n, nLC
TYPE(Matrix_t), POINTER :: A
TYPE(Solver_t), POINTER :: Solver
!call resettimer('prec')
IF(.NOT.Precondition) THEN
u=v
......@@ -1829,6 +1856,7 @@ END SUBROUTINE FetiProject
IF (.NOT. CPG .AND. nz>0) THEN
u(nLC+1:nLC+nz)=v(nLC+1:nLC+nz)
END IF
!call checktimer('prec',delete=.true.)
!------------------------------------------------------------------------------
END SUBROUTINE FetiPrec
!------------------------------------------------------------------------------
......
......@@ -81,8 +81,9 @@ SUBROUTINE ForceCompute( Model,Solver,dt,TransientSimulation )
!------------------------------------------------------------------------------
REAL(KIND=dp), ALLOCATABLE :: Pressure(:), Velocity(:,:), Viscosity(:)
REAL(KIND=dp), ALLOCATABLE :: ShearData(:,:)
REAL(KIND=dp), ALLOCATABLE :: MomentAbout(:,:), Forces(:,:), Moments(:,:), Areas(:)
REAL(KIND=dp) :: Force(3), Moment(3), Area, ShearStress
REAL(KIND=dp), ALLOCATABLE :: Forces(:,:), Moments(:,:), Areas(:)
REAL(KIND=dp) :: Force(3), Moment(3), MomentAbout(3),Area, ShearStress
REAL(KIND=dp), POINTER :: mWork(:,:)
LOGICAL :: Stat, CalculateMoment, ViscousForce, Compressible, SumForces
LOGICAL :: ShearOutput
INTEGER :: i,j,k,n,pn,t,dim
......@@ -93,6 +94,7 @@ SUBROUTINE ForceCompute( Model,Solver,dt,TransientSimulation )
TYPE(ValueList_t), POINTER :: Material
TYPE(Nodes_t) :: ElementNodes, ParentNodes
TYPE(Element_t), POINTER :: CurrentElement, Parent
TYPE(Valuelist_t), POINTER :: BC
CHARACTER(LEN=MAX_NAME_LEN) :: ShearFilename, MessageL, ViscosityFlag
CHARACTER(LEN=MAX_NAME_LEN) :: CompressibilityFlag, BoundaryName, VariableName
......@@ -113,7 +115,6 @@ SUBROUTINE ForceCompute( Model,Solver,dt,TransientSimulation )
ALLOCATE( ElementNodes % x(n), ElementNodes % y(n), ElementNodes % z(n))
ALLOCATE( ParentNodes % x(n), ParentNodes % y(n), ParentNodes % z(n) )
ALLOCATE( Pressure( n ), Viscosity( n ), Velocity( 3, n ) )
ALLOCATE( MomentAbout(3,n) )
ALLOCATE( Forces(Model % NumberOfBCs,3) )
ALLOCATE( Moments(Model % NumberOfBCs,3) )
......@@ -177,8 +178,9 @@ SUBROUTINE ForceCompute( Model,Solver,dt,TransientSimulation )
IF ( CurrentElement % TYPE % ElementCode == 101 ) CYCLE
DO k=1, Model % NumberOfBCs
BC => Model % BCs(k) % Values
IF ( Model % BCs(k) % Tag /= CurrentElement % BoundaryInfo % Constraint ) CYCLE
IF ( .NOT. ListGetLogical(Model % BCs(k) % Values,'Calculate Fluidic Force',stat ) ) CYCLE
IF ( .NOT. ListGetLogical(BC,'Calculate Fluidic Force',stat ) ) CYCLE
ElementNodes % x(1:n) = Mesh % Nodes % x(NodeIndexes)
ElementNodes % y(1:n) = Mesh % Nodes % y(NodeIndexes)
......@@ -219,16 +221,16 @@ SUBROUTINE ForceCompute( Model,Solver,dt,TransientSimulation )
Viscosity(1:pn) = ListGetReal( Material, 'Viscosity', pn, Parent % NodeIndexes )
MomentAbout(1,1:n) = ListGetReal( Model % BCs(k) % Values, &
'Moment About 1', n, NodeIndexes, stat )
MomentAbout(2,1:n) = ListGetReal( Model % BCs(k) % Values, &
'Moment About 2', n, NodeIndexes, CalculateMoment )
CalculateMoment = stat .OR. CalculateMoment
MomentAbout(3,1:n) = ListGetReal( Model % BCs(k) % Values, &
'Moment About 3', n, NodeIndexes, stat )
CalculateMoment = stat .OR. CalculateMoment
mWork => ListGetConstRealArray( BC,'Moment About',CalculateMoment)
IF( CalculateMoment ) THEN
MomentAbout(1:dim) = mWork(1:dim,1)
ELSE
MomentAbout(1) = ListGetCReal( BC,'Moment About 1', stat )
MomentAbout(2) = ListGetCReal( BC,'Moment About 2', CalculateMoment )
CalculateMoment = stat .OR. CalculateMoment
MomentAbout(3) = ListGetCReal( BC,'Moment About 3', stat )
CalculateMoment = stat .OR. CalculateMoment
END IF
Velocity = 0.0d0
DO i=1,pn
......@@ -273,11 +275,7 @@ SUBROUTINE ForceCompute( Model,Solver,dt,TransientSimulation )
WRITE( BoundaryName, '("")')
ELSE
IF ( .NOT. ListGetLogical(Model % BCs(k) % Values,'Calculate Fluidic Force',stat ) ) CYCLE
IF(Model % NumberOfBCs < 10) THEN
WRITE( BoundaryName, '("bc ",I1)') k
ELSE
WRITE( BoundaryName, '("bc ",I2)') k
END IF
WRITE( BoundaryName, '("bc ",I0)') k
END IF
nlen = LEN_TRIM(BoundaryName)
......@@ -292,7 +290,7 @@ SUBROUTINE ForceCompute( Model,Solver,dt,TransientSimulation )
IF ( CalculateMoment ) THEN
WRITE( Message, &
'("Moment about (",ES9.3E1,",",ES9.3E1,",",ES9.3E1,") is:",3Es14.6E2)') &
MomentAbout(1:3,1), Moments(k,1:3)
MomentAbout(1:3), Moments(k,1:3)
CALL Info( 'ForceCompute', Message, Level=4 )
CALL ListAddConstReal( Model % Simulation, &
......@@ -377,7 +375,7 @@ CONTAINS
SUBROUTINE ForceIntegrate( Force, Moment, MomentAbout, ViscousForce, &
Compressible, Area, ShearStress )