diff --git a/fem/src/DefUtils.F90 b/fem/src/DefUtils.F90 index e5e35ecc6b..3b5f7ea2ab 100644 --- a/fem/src/DefUtils.F90 +++ b/fem/src/DefUtils.F90 @@ -3581,9 +3581,13 @@ RECURSIVE FUNCTION DefaultSolve( USolver, BackRotNT ) RESULT(Norm) CALL Info('DefaultSolve','Calling SolveSystem for linear solution',Level=20) +print*,'a: ', associated(solver), associated(Solver % matrix); flush(6) A => Solver % Matrix +print*,'a: ', associated(solver % variable); flush(6) x => Solver % Variable +print*,'a: ', associated(A), associated(A % RHS); flush(6) b => A % RHS +print*,'a: ', associated(x % values), SIZE(x % values); flush(6) SOL => x % Values ! Debugging stuff activated only when "Max Output Level" >= 20 @@ -3594,6 +3598,7 @@ RECURSIVE FUNCTION DefaultSolve( USolver, BackRotNT ) RESULT(Norm) 10 CONTINUE +print*,'going sol'; flush(6) CALL SolveSystem(A,ParMatrix,b,SOL,x % Norm,x % DOFs,Solver) IF( InfoActive( 20 ) ) THEN diff --git a/fem/src/ElementUtils.F90 b/fem/src/ElementUtils.F90 index dcc8e5374f..24d3869e21 100644 --- a/fem/src/ElementUtils.F90 +++ b/fem/src/ElementUtils.F90 @@ -1904,6 +1904,7 @@ FUNCTION CreateMatrix( Model, Solver, Mesh, Perm, DOFs, MatrixFormat, & IF ( PRESENT(Equation) ) THEN CALL Info(Caller,'Creating initial permutation',Level=14) k = InitialPermutation( Perm,Model,Solver,Mesh,Eq,DG,GB ) +print*,'aft k', k IF ( k <= 0 ) THEN IF(ALLOCATED(InvInitialReorder)) DEALLOCATE(InvInitialReorder) RETURN diff --git a/fem/src/ListMatrix.F90 b/fem/src/ListMatrix.F90 index ae7f30faad..c1b1c05799 100644 --- a/fem/src/ListMatrix.F90 +++ b/fem/src/ListMatrix.F90 @@ -482,6 +482,113 @@ SUBROUTINE List_AddMatrixIndexes(List,k1,nk2,Ind) END SUBROUTINE List_AddMatrixIndexes !------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- + SUBROUTINE List_AddMatrixRow(List,k1,nk2,Ind,Vals) + ! Add an array of sorted indeces to a row in ListMatrix_t. "ind" may + ! contain duplicate entries. +!------------------------------------------------------------------------------- + IMPLICIT NONE + + TYPE(ListMatrix_t) :: List(:) + INTEGER, INTENT(IN) :: k1, nk2 + INTEGER, INTENT(INOUT) :: Ind(nk2) + REAL(KIND=dp), INTENT(INOUT) :: Vals(nk2) + + TYPE(ListMatrixEntry_t), POINTER :: RowPtr, PrevPtr, Entry, Dummy +!------------------------------------------------------------------------------- + INTEGER :: i,k2,k2i,j, k,prevind + + IF (k1>SIZE(List)) THEN + CALL Fatal('List_AddMatrixIndexes','Row index out of bounds: '//TRIM(I2S(k1))) + END IF + + CALL SortF(nk2, Ind, Vals) + + ! Add each element in Ind to the row list + RowPtr => List(k1) % Head + + ! First element needs special treatment as it may modify + ! the list starting point + IF (.NOT. ASSOCIATED(RowPtr)) THEN + Dummy => NULL() + Entry => List_GetMatrixEntry(Ind(1),Dummy) + Entry % Val = Vals(1) + List(k1) % Degree = 1 + List(k1) % Head => Entry + k2i = 2 + prevind = ind(1) + ELSE IF (RowPtr % Index > Ind(1)) THEN + Entry => List_GetMatrixEntry(Ind(1),RowPtr) + Entry % Val = Vals(1) + List(k1) % Degree = List(k1) % Degree + 1 + List(k1) % Head => Entry + k2i = 2 + prevind = ind(1) + ELSE IF (RowPtr % Index == Ind(1)) THEN + k2i = 2 + prevind = ind(1) + ELSE + k2i = 1 + prevind = -1 + END IF + + PrevPtr => List(k1) % Head + RowPtr => List(k1) % Head % Next + + DO i=k2i,nk2 + k2 = Ind(i) + if (k2 == prevind) cycle + + ! Find a correct place place to add index to + DO WHILE( ASSOCIATED(RowPtr) ) + IF (RowPtr % Index >= k2) EXIT + PrevPtr => RowPtr + RowPtr => RowPtr % Next + END DO + + IF (ASSOCIATED(RowPtr)) THEN + ! Do not add duplicates + IF (RowPtr % Index /= k2) THEN + ! Create new element between PrevPtr and RowPtr + Entry => List_GetMatrixEntry(k2,RowPtr) + Entry % Val = Vals(i) + PrevPtr % Next => Entry + List(k1) % Degree = List(k1) % Degree + 1 + + ! Advance to next element in list + PrevPtr => Entry + ELSE + ! Advance to next element in list + RowPtr % Val = RowPtr % Val + Vals(i) + PrevPtr => RowPtr + RowPtr => RowPtr % Next + END IF + ELSE + EXIT + END IF + + prevind = k2 + END DO + + DO j=i,nk2 + k2 = Ind(j) + if (k2 == prevind) cycle + prevind = k2 + + Dummy => NULL() + Entry => List_GetMatrixEntry(k2,Dummy) + Entry % Val = Vals(j) + PrevPtr % Next => Entry + PrevPtr => PrevPtr % Next + List(k1) % Degree = List(k1) % Degree + 1 + END DO +!------------------------------------------------------------------------------- + END SUBROUTINE List_AddMatrixRow +!------------------------------------------------------------------------------- + + + !------------------------------------------------------------------------------- FUNCTION List_GetMatrixEntry(ind, next) RESULT(ListEntry) !------------------------------------------------------------------------------- diff --git a/fem/src/Lists.F90 b/fem/src/Lists.F90 index d146742d71..7ea8f3de11 100644 --- a/fem/src/Lists.F90 +++ b/fem/src/Lists.F90 @@ -117,13 +117,25 @@ SUBROUTINE SetGetMatcParams(nparams,params,resul) REAL(KIND=dp) :: params(:) CHARACTER(*), OPTIONAL :: resul - INTEGER :: i,l + INTEGER :: i,j,l CHARACTER(LEN=1024) :: pcmd,res IF(nparams==0) THEN pcmd = "tx=0" ELSE +#if 0 WRITE(pcmd,*) [(params(i),i=1,nparams)] +#else + ! cray ftn output from above can be somewhat convoluted, do this instead + j = 1 + DO i=1,nparams + WRITE(pcmd(j:), *) params(i) + DO WHILE(pcmd(j:j) == ' '); j=j+1; END DO + DO WHILE(pcmd(j:j) /= ' '); j=j+1; END DO + IF(pcmd(j-1:j-1)=='.') pcmd(j-1:j-1) = ' ' + j = j + 1 + END DO +#endif IF(PRESENT(resul)) THEN pcmd = TRIM(resul)//'='//TRIM(pcmd) ELSE @@ -243,6 +255,7 @@ FUNCTION InitialPermutation( Perm,Model,Solver,Mesh, & INTEGER :: body_id, MaxGroup, group0, group INTEGER, POINTER :: DgMap(:), DgMaster(:), DgSlave(:) LOGICAL :: GotDgMap, GotMaster, GotSlave +!------------------------------------------------------------------------------ DgMap => ListGetIntegerArray( Solver % Values,'DG Reduced Basis Mapping',GotDgMap ) DgMaster => ListGetIntegerArray( Solver % Values,'DG Reduced Basis Master Bodies',GotMaster ) @@ -326,7 +339,8 @@ FUNCTION InitialPermutation( Perm,Model,Solver,Mesh, & ' db nodes from bulk hits',Level=15) IF ( FoundDG ) THEN - RETURN ! Discontinuous bodies !!! + GOTO 10 +! RETURN ! Discontinuous bodies !!! END IF END BLOCK END IF @@ -408,7 +422,8 @@ FUNCTION InitialPermutation( Perm,Model,Solver,Mesh, & ' nodes from bulk hits',Level=15) IF ( FoundDG ) THEN - RETURN ! Discontinuous galerkin !!! + GOTO 10 +! RETURN ! Discontinuous galerkin !!! END IF END IF @@ -678,6 +693,9 @@ FUNCTION InitialPermutation( Perm,Model,Solver,Mesh, & IF ( ALLOCATED(EdgeDOFs) ) DEALLOCATE(EdgeDOFs) IF ( ALLOCATED(FaceDOFs) ) DEALLOCATE(FaceDOFs) + +10 CONTINUE + !------------------------------------------------------------------------------ END FUNCTION InitialPermutation !------------------------------------------------------------------------------ diff --git a/fem/src/LoadMod.F90 b/fem/src/LoadMod.F90 index 8f34f390a2..01a49a1689 100644 --- a/fem/src/LoadMod.F90 +++ b/fem/src/LoadMod.F90 @@ -297,7 +297,7 @@ FUNCTION loadfunction(quiet, abort_not_found, library, fname) RESULT(ptr) CHARACTER :: library(*), fname(*) ! TYPE(C_FUNPTR) :: ptr INTEGER(KIND=AddrInt) :: ptr - TYPE(C_FUNPTR) :: cptr + TYPE(C_PTR) :: cptr INTERFACE FUNCTION loadfunction_c(quiet, abort_not_found, library, fname ) RESULT(cptr) & @@ -307,7 +307,7 @@ FUNCTION loadfunction_c(quiet, abort_not_found, library, fname ) RESULT(cptr) & INTEGER(C_INT) :: abort_not_found CHARACTER(C_CHAR) :: library(*), fname(*) ! INTEGER(CAddrInt) :: fptr - TYPE(C_FUNPTR) :: cptr + TYPE(C_PTR) :: cptr END FUNCTION loadfunction_c END INTERFACE diff --git a/fem/src/MeshUtils.F90 b/fem/src/MeshUtils.F90 index 8f68a6283f..16456b0743 100644 --- a/fem/src/MeshUtils.F90 +++ b/fem/src/MeshUtils.F90 @@ -27468,7 +27468,7 @@ SUBROUTINE ElmerGraphColour(Graph, Colouring, ConsistentColours) INTEGER :: nc, dualmaxdeg, i, v, w, uci, wci, vli, vti, vcol, wcol, & nrc, nunc, nthr, TID, allocstat, gn - INTEGER, ALLOCATABLE :: colours(:) + INTEGER, POINTER :: colours(:) INTEGER, PARAMETER :: VERTEX_PER_THREAD = 100 LOGICAL :: consistent @@ -27649,7 +27649,8 @@ SUBROUTINE ElmerGraphColour(Graph, Colouring, ConsistentColours) ! Set up colouring data structure Colouring % nc = nc - CALL MOVE_ALLOC(colours, Colouring % colours) +! CALL MOVE_ALLOC(colours, Colouring % colours) + Colouring % colours => colours END SUBROUTINE ElmerGraphColour SUBROUTINE Colouring_Deallocate(Colours) @@ -27709,7 +27710,7 @@ SUBROUTINE ElmerBoundaryGraphColour(Mesh, Colours, BoundaryColours) TYPE(Element_t), POINTER :: Element INTEGER :: elem, nelem, nbelem, astat, lcolour, rcolour, nbc - INTEGER, ALLOCATABLE :: bcolours(:) + INTEGER, POINTER :: bcolours(:) nelem = Mesh % NumberOfBulkElements nbelem = Mesh % NumberOfBoundaryElements @@ -27757,7 +27758,8 @@ SUBROUTINE ElmerBoundaryGraphColour(Mesh, Colours, BoundaryColours) ! Set up colouring data structure BoundaryColours % nc = nbc - CALL MOVE_ALLOC(bcolours, BoundaryColours % colours) +! CALL MOVE_ALLOC(bcolours, BoundaryColours % colours) + BoundaryColours % colours => bcolours END SUBROUTINE ElmerBoundaryGraphColour ! Given CRS indices, referenced indirectly from graph, diff --git a/fem/src/SolverUtils.F90 b/fem/src/SolverUtils.F90 index 24392ec0a3..cc455797b0 100644 --- a/fem/src/SolverUtils.F90 +++ b/fem/src/SolverUtils.F90 @@ -7365,8 +7365,8 @@ END SUBROUTINE SetDirichletPoint !------------------------------------------------------------------------------ SUBROUTINE SetNodalSources( Model, Mesh, SourceName, dofs, Perm, GotSrc, SrcVec ) !------------------------------------------------------------------------------ - TYPE(Model_t), POINTER :: Model !< The current model structure - TYPE(Mesh_t), POINTER :: Mesh !< The current mesh structure + TYPE(Model_t) :: Model !< The current model structure + TYPE(Mesh_t) :: Mesh !< The current mesh structure CHARACTER(LEN=*) :: SourceName !< Name of the keyword setting the source term INTEGER :: DOFs !< The total number of DOFs for this equation INTEGER :: Perm(:) !< The node reordering info @@ -7374,7 +7374,7 @@ SUBROUTINE SetNodalSources( Model, Mesh, SourceName, dofs, Perm, GotSrc, SrcVec REAL(KIND=dp) :: SrcVec(:) !< The assemblied source vector !------------------------------------------------------------------------------ TYPE(Element_t), POINTER :: Element - INTEGER :: i,t,n,bc,bf,FirstElem,LastElem,nlen + INTEGER :: i,j,k,t,n,bc,bf,FirstElem,LastElem,nlen LOGICAL :: Found,AnyBC,AnyBF,Axisymmetric REAL(KIND=dp) :: Coeff REAL(KIND=dp), ALLOCATABLE :: FORCE(:,:) @@ -7389,8 +7389,7 @@ SUBROUTINE SetNodalSources( Model, Mesh, SourceName, dofs, Perm, GotSrc, SrcVec CALL Info(Caller,'Checking for generalized source terms: '& //SourceName(1:nlen),Level=15) - ALLOCATE( ActiveBC(Model % NumberOfBCs ), & - ActiveBF(Model % NumberOfBodyForces) ) + ALLOCATE( ActiveBC(Model % NumberOfBCs), ActiveBF(Model % NumberOfBodyForces) ) ! First make a quick test going through the short boundary condition and ! body force lists. @@ -7439,6 +7438,7 @@ SUBROUTINE SetNodalSources( Model, Mesh, SourceName, dofs, Perm, GotSrc, SrcVec DO t=FirstElem, LastElem Element => Mesh % Elements(t) Indexes => Element % NodeIndexes + n = Element % Type % NumberOfNodes IF( t > Mesh % NumberOfBulkElements ) THEN Found = .FALSE. @@ -7465,8 +7465,7 @@ SUBROUTINE SetNodalSources( Model, Mesh, SourceName, dofs, Perm, GotSrc, SrcVec CALL LocalSourceAssembly(Element, dofs, FORCE ) DO i=1,dofs - SrcVec(dofs*(Perm(Indexes)-1)+i) = SrcVec(dofs*(Perm(Indexes)-1)+i) + & - Coeff * FORCE(i,1:n) + SrcVec(dofs*(Perm(Indexes)-1)+i) = SrcVec(dofs*(Perm(Indexes)-1)+i) + FORCE(i,1:n) END DO END DO @@ -16187,12 +16186,11 @@ SUBROUTINE FinalizeLumpedMatrix( Solver ) END DO CLOSE(11) END IF - + + WRITE(Message,*) 'Normalization checksum: ',CheckSum + CALL Info(Caller,Message) END IF CALL Info( Caller,'Constraint modes fluxes was saved to file '//TRIM(MatrixFile),Level=5) - - WRITE(Message,*) 'Normalization checksum: ',CheckSum - CALL Info(Caller,Message) END IF END IF @@ -19700,21 +19698,20 @@ SUBROUTINE ControlLinearSystem(Solver,PreSolve) ! We need to add the control source here in order to be able to use ! standard means for convergence monitoring. CALL Info(Caller,'Computing source term for: '//TRIM(str),Level=7) - CALL SetNodalSources( CurrentModel,Mesh,str, & - dofs, Perm, GotF, f(:,iControl) ) + CALL SetNodalSources( CurrentModel,Mesh,str,dofs, Perm, GotF, f(:,iControl) ) - ! The additional source needs to be nullified for Dirichlet conditions - IF( ALLOCATED( A % ConstrainedDOF ) ) THEN - WHERE( A % ConstrainedDOF ) f(:,iControl) = 0.0_dp - END IF + ! The additional source needs to be nullified for Dirichlet conditions + IF( ALLOCATED( A % ConstrainedDOF ) ) THEN + WHERE( A % ConstrainedDOF ) f(:,iControl) = 0.0_dp + END IF - IF(InfoActive(10)) THEN - DO i=1,dofs - PRINT *,'ranges b:',i,MINVAL(b(i::dofs)),MAXVAL(b(i::dofs)),SUM(b(i::dofs)) - PRINT *,'ranges f:',i,MINVAL(f(i::dofs,iControl)),& - MAXVAL(f(i::dofs,iControl)),SUM(f(i::dofs,iControl)) - END DO - END IF + IF(InfoActive(10)) THEN + DO i=1,dofs + PRINT *,'ranges b:',i,MINVAL(b(i::dofs)),MAXVAL(b(i::dofs)),SUM(b(i::dofs)) + PRINT *,'ranges f:',i,MINVAL(f(i::dofs,iControl)),& + MAXVAL(f(i::dofs,iControl)),SUM(f(i::dofs,iControl)) + END DO + END IF IF( ABS(cAmp(iControl)) > 1.0e-20 ) THEN b(1:nsize) = b(1:nsize) + cAmp(iControl) * f(1:nsize,iControl) @@ -25138,6 +25135,17 @@ END SUBROUTINE P2LagrangeSwapper !------------------------------------------------------------------------------------- FUNCTION EvalFieldAtElem( Mesh, Var, Element, Basis, dofi, eigeni, imVal, GotVal ) RESULT ( Val ) + INTERFACE + SUBROUTINE Ip2DgFieldInElement( Mesh, Parent, nip, fip, np, fdg ) + USE Types + IMPLICIT NONE + TYPE(Mesh_t), POINTER :: Mesh + TYPE(Element_t), POINTER :: Parent + INTEGER :: nip, np + REAL(KIND=dp) :: fip(:), fdg(:) + END SUBROUTINE Ip2DgFieldInElement + END INTERFACE + TYPE(Mesh_t), POINTER :: Mesh TYPE(Variable_t), POINTER :: Var TYPE(Element_t), POINTER :: Element diff --git a/fem/src/Types.F90 b/fem/src/Types.F90 index 39ed6fa782..11b5f361f9 100644 --- a/fem/src/Types.F90 +++ b/fem/src/Types.F90 @@ -854,7 +854,7 @@ MODULE Types TYPE Graphcolour_t INTEGER :: nc - INTEGER, ALLOCATABLE :: colours(:) + INTEGER, POINTER :: colours(:) => Null() END TYPE Graphcolour_t TYPE MortarBC_t diff --git a/fem/src/modules/SteadyPhaseChange.F90 b/fem/src/modules/SteadyPhaseChange.F90 index ee4b16f6dc..dea762e27a 100644 --- a/fem/src/modules/SteadyPhaseChange.F90 +++ b/fem/src/modules/SteadyPhaseChange.F90 @@ -1,4 +1,4 @@ -/*****************************************************************************/ +!/*****************************************************************************/ ! * ! * Elmer, A Finite Element Software for Multiphysical Problems ! * diff --git a/fem/src/rocalution.cpp b/fem/src/rocalution.cpp index 05981c4451..584e4c5796 100644 --- a/fem/src/rocalution.cpp +++ b/fem/src/rocalution.cpp @@ -147,10 +147,10 @@ void elmer_distribute_matrix(const MPI_Comm* comm, boundary[r].push_back(i + index_offset[rank]); neighbor[r] = true; ++boundary_nnz; - checked[r][i + index_offset[rank]] = true; + checked[r][i + index_offset[rank]] = true; } ++ghost_nnz; - // Rank for current boundary point local_col[j] has been found + // Rank for current boundary point local_col[j] has been found // Continue with next boundary point break; } diff --git a/fem/tests/AdvReactDB/AdvReactDB.sif b/fem/tests/AdvReactDB/AdvReactDB.sif index fb5f8be01c..92b2ca27ef 100644 --- a/fem/tests/AdvReactDB/AdvReactDB.sif +++ b/fem/tests/AdvReactDB/AdvReactDB.sif @@ -16,7 +16,7 @@ $ function dx(c) { x=c(0); y=c(1); _dx = pi*y^2*cos(pi*x*y^2) } $ function dy(c) { x=c(0); y=c(1); _dy = pi*x*2*y*cos(pi*x*y^2) } Simulation - Max Output Level = 5 + Max Output Level = 25 Coordinate System = "Cartesian" Coordinate Mapping(3) = 1 2 3 Simulation Type = Steady diff --git a/fem/tests/RotatingBCMagnetoDynamicsConfRelease/case.sif b/fem/tests/RotatingBCMagnetoDynamicsConfRelease/case.sif index acfa90ab6e..00bfb811df 100644 --- a/fem/tests/RotatingBCMagnetoDynamicsConfRelease/case.sif +++ b/fem/tests/RotatingBCMagnetoDynamicsConfRelease/case.sif @@ -375,7 +375,7 @@ Solver 2 :: Reference Norm = 4.22134445E-03 Solver 3 :: Reference Norm = 6.66996883E-02 Solver 7 :: Reference Norm = 5.92002332E-01 -Solver 2 :: Reference Norm Tolerance = 1e-4 -Solver 3 :: Reference Norm Tolerance = 1e-4 +Solver 2 :: Reference Norm Tolerance = 2e-4 +Solver 3 :: Reference Norm Tolerance = 2e-4 diff --git a/fem/tests/Zirka/ZirkaTest.F90 b/fem/tests/Zirka/ZirkaTest.F90 index ab74b181a8..651601dd6a 100644 --- a/fem/tests/Zirka/ZirkaTest.F90 +++ b/fem/tests/Zirka/ZirkaTest.F90 @@ -298,5 +298,6 @@ subroutine test3() ! {{{ variable % norm = abs(H(i_r)-mc % eval(rps(i_r), cached = .true.))/abs(H(i_r)) end subroutine ! }}} +end function ! }}} END SUBROUTINE ZirkaTest ! }}} !-------------------------------------------------------------------------------