Commit bd325599 authored by Adam Powell's avatar Adam Powell

Imported Upstream version 5.5.0.svn.5210.dfsg

parent 11a413b2
This diff is collapsed.
......@@ -222,12 +222,12 @@
<material name="Water (frozen)" >
<parameter name="Density" >910.0</parameter>
<parameter name="Heat conductivity" >Variable Temperature; Real MATC "9.828*exp(-5.7E-03*tx)"</parameter>
<parameter name="Heat conductivity" >Variable Temperature; Real MATC "9.828*exp(-5.7E-03*(tx+273.15))"</parameter>
<parameter name="Heat capacity" >Variable Temperature; Real MATC "146.3+(7.253*(tx+273.15))"</parameter>
<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)))^(-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>
......
***** ElmerGrid input file for structured grid generation *****
Version = 210903
Coordinate System = Cartesian 2D
Subcell Divisions in 2D = 3 1
Subcell Limits 1 = 0 0.1 0.2 0.4
Subcell Limits 2 = 0 6.28318531
Material Structure in 2D
1 1 2
End
Materials Interval = 1 2
Boundary Definitions
! type out int double of the boundaries
1 -2 1 1
2 -2 2 1
End
Numbering = Horizontal
Element Degree = 1
Element Innernodes = False
Triangles = False
Reference Density = 0.03
Coordinate Ratios = 1
Element Ratios 1 = 1 1 3.0
Element Ratios 2 = 1
Element Densities 1 = 1 1 1
Element Densities 2 = 0.5
Revolve Blocks = 4
!Revolve Radius = 2
Revolve Improve = 0.5
Revolve Curve Direct = 0.0
Revolve Curve Angle = 360.0
Revolve Curve Radius = 1.0
Merge Nodes = 1.0e-6
\ No newline at end of file
# geometry of a toy glacier on bedrock in in2d format of netgen
# start with keyword splinecurves2dv2
splinecurves2dv2
# next the refinement factor
1.2
points
1 0 0
2 7000.0 100.0
3 7000.0 1200.0
4 2000.0 1040.0
5 0.0 1000.0
6 7000.0 -2000.0
7 -2000.0 -2000.0
8 -2000.0 -40.0
segments
1 2 2 1 2 -bc=1
1 0 2 2 3 -bc=2
1 0 2 3 4 -bc=3
1 0 3 4 5 1 -bc=3
0 2 2 2 6 -bc=4
0 2 2 6 7 -bc=5
0 2 2 7 8 -bc=6
0 2 2 8 1 -bc=7
materials
1 ice -maxh=50.0
2 rock -maxh=200.0
# geometry of a toy glacier in in2d format of netgen
# start with keyword splinecurves2dv2
splinecurves2dv2
# next the refinement factor
20
points
1 0 0
2 7000.0 100.0
3 7000.0 1200.0
4 2000.0 1040.0
5 0.0 1000.0
segments
1 0 2 1 2 -bc=1
1 0 2 2 3 -bc=2
1 0 2 3 4 -bc=3
1 0 3 4 5 1 -bc=3
......@@ -2339,7 +2339,7 @@ fi
# Define the identity of the package.
PACKAGE=fem
VERSION=6.1
VERSION=6.2
cat >>confdefs.h <<_ACEOF
......
This diff is collapsed.
......@@ -59,31 +59,142 @@ CONTAINS
!------------------------------------------------------------------------
SUBROUTINE EdgeBasis( Element, WBasis, RotWBasis, Basis, dBasisdx )
!------------------------------------------------------------------------
TYPE(Element_t) :: Element
TYPE(Element_t),TARGET :: Element
REAL(KIND=dp) :: WBasis(:,:), RotWBasis(:,:), Basis(:), dBasisdx(:,:)
!------------------------------------------------------------------------
TYPE(Element_t),POINTER :: Edge
TYPE(Mesh_t), POINTER :: Mesh
LOGICAL :: Parallel
INTEGER :: i,j,k,nj,nk
TYPE(Nodes_t), SAVE :: Nodes
REAL(KIND=dp) :: u,v,w,dudx(3,3),du(3),Base,dBase(3),tBase(3), &
rBase(3),triBase(3),dtriBase(3,3)
LOGICAL :: Parallel,stat
INTEGER :: i,j,k,n,nj,nk,i1,i2
INTEGER, POINTER :: EdgeMap(:,:)
!------------------------------------------------------------------------
Mesh => GetMesh()
Parallel = ParEnv % PEs>1
SELECT CASE(GetElementFamily(Element))
CASE(4,7,8)
n = Element % Type % NumberOfNodes
u = SUM(Basis(1:n)*Element % Type % NodeU(1:n))
v = SUM(Basis(1:n)*Element % Type % NodeV(1:n))
w = SUM(Basis(1:n)*Element % Type % NodeW(1:n))
dudx(1,:) = MATMUL(Element % Type % NodeU(1:n),dBasisdx(1:n,:))
dudx(2,:) = MATMUL(Element % Type % NodeV(1:n),dBasisdx(1:n,:))
dudx(3,:) = MATMUL(Element % Type % NodeW(1:n),dBasisdx(1:n,:))
triBase(1) = 1-u-v
triBase(2) = u
triBase(3) = v
dtriBase(1,:) = -dudx(1,:)-dudx(2,:)
dtriBase(2,:) = dudx(1,:)
dtriBase(3,:) = dudx(2,:)
END SELECT
EdgeMap => GetEdgeMap(GetElementFamily(Element))
DO i=1,SIZE(Edgemap,1)
j = EdgeMap(i,1); k = EdgeMap(i,2)
WBasis(i,:) = Basis(j)*dBasisdx(k,:) - Basis(k)*dBasisdx(j,:)
RotWBasis(i,1) = 2.0_dp * ( dBasisdx(j,2) * dBasisdx(k,3) - &
dBasisdx(j,3) * dBasisdx(k,2) )
SELECT CASE(GetElementFamily(Element))
CASE(3,5)
WBasis(i,:) = Basis(j)*dBasisdx(k,:) - Basis(k)*dBasisdx(j,:)
RotWBasis(i,1) = 2.0_dp * ( dBasisdx(j,2) * dBasisdx(k,3) - &
dBasisdx(j,3) * dBasisdx(k,2) )
RotWBasis(i,2) = 2.0_dp * ( dBasisdx(j,3) * dBasisdx(k,1) - &
dBasisdx(j,1) * dBasisdx(k,3) )
RotWBasis(i,3) = 2.0_dp * ( dBasisdx(j,1) * dBasisdx(k,2) - &
dBasisdx(j,2) * dBasisdx(k,1) )
CASE(7)
SELECT CASE(i)
CASE(1)
j=1;k=2; Base=1-w; dBase=-dudx(3,:)
CASE(2)
j=2;k=3; Base=1-w; dBase=-dudx(3,:)
CASE(3)
j=3;k=1; Base=1-w; dBase=-dudx(3,:)
CASE(4)
j=1;k=2; Base=1+w; dBase= dudx(3,:)
CASE(5)
j=2;k=3; Base=1+w; dBase= dudx(3,:)
CASE(6)
j=3;k=1; Base=1+w; dBase= dudx(3,:)
CASE(7)
Base=triBase(1); dBase=dtriBase(1,:); du=dudx(3,:)
CASE(8)
Base=triBase(2); dBase=dtriBase(2,:); du=dudx(3,:)
CASE(9)
Base=triBase(3); dBase=dtriBase(3,:); du=dudx(3,:)
END SELECT
IF(i<=6) THEN
tBase = (triBase(j)*dtriBase(k,:)-triBase(k)*dtriBase(j,:))
rBase(1) = 2*Base*(dtriBase(j,2)*dtriBase(k,3)-dtriBase(k,2)*dtriBase(j,3)) + &
dBase(2)*tBase(3) - dBase(3)*tBase(2)
rBase(2) = 2*Base*(dtriBase(j,3)*dtriBase(k,1)-dtriBase(k,3)*dtriBase(j,1)) + &
dBase(3)*tBase(1) - dBase(1)*tBase(3)
RotWBasis(i,2) = 2.0_dp * ( dBasisdx(j,3) * dBasisdx(k,1) - &
dBasisdx(j,1) * dBasisdx(k,3) )
rBase(3) = 2*Base*(dtriBase(j,1)*dtriBase(k,2)-dtriBase(k,1)*dtriBase(j,2)) + &
dBase(1)*tBase(2) - dBase(2)*tBase(1)
RotWBasis(i,3) = 2.0_dp * ( dBasisdx(j,1) * dBasisdx(k,2) - &
dBasisdx(j,2) * dBasisdx(k,1) )
RotWBasis(i,:)=rBase
WBasis(i,:)=tBase*Base
ELSE
WBasis(i,:)=Base*du/2
RotWBasis(i,1)=(dBase(2)*du(3) - dBase(3)*du(2))/2
RotWBasis(i,2)=(dBase(3)*du(1) - dBase(1)*du(3))/2
RotWBasis(i,3)=(dBase(1)*du(2) - dBase(2)*du(1))/2
END IF
CASE(4,8)
SELECT CASE(i)
CASE(1)
du=dudx(1,:); Base=(1-v)*(1-w)
dBase(:)=-dudx(2,:)*(1-w)-(1-v)*dudx(3,:)
CASE(3)
du=dudx(1,:); Base=(1+v)*(1-w)
dBase(:)=dudx(2,:)*(1-w)-(1+v)*dudx(3,:)
CASE(5)
du=dudx(1,:); Base=(1-v)*(1+w)
dBase(:)=-dudx(2,:)*(1+w)+(1-v)*dudx(3,:)
CASE(7)
du=dudx(1,:); Base=(1+v)*(1+w)
dBase(:)=dudx(2,:)*(1+w)+(1+v)*dudx(3,:)
CASE(2)
du=dudx(2,:); Base=(1+u)*(1-w)
dBase(:)=dudx(1,:)*(1-w)-(1+u)*dudx(3,:)
CASE(4)
du=dudx(2,:); Base=(1-u)*(1-w)
dBase(:)=-dudx(1,:)*(1-w)-(1-u)*dudx(3,:)
CASE(6)
du=dudx(2,:); Base=(1+u)*(1+w)
dBase(:)=dudx(1,:)*(1+w)+(1+u)*dudx(3,:)
CASE(8)
du=dudx(2,:); Base=(1-u)*(1+w)
dBase(:)=-dudx(1,:)*(1+w)+(1-u)*dudx(3,:)
CASE(9)
du=dudx(3,:); Base=(1-u)*(1-v)
dBase(:)=-dudx(1,:)*(1-v)-(1-u)*dudx(2,:)
CASE(10)
du=dudx(3,:); Base=(1+u)*(1-v)
dBase(:)=dudx(1,:)*(1-v)-(1+u)*dudx(2,:)
CASE(11)
du=dudx(3,:); Base=(1+u)*(1+v)
dBase(:)=dudx(1,:)*(1+v)+(1+u)*dudx(2,:)
CASE(12)
du=dudx(3,:); Base=(1-u)*(1+v)
dBase(:)=-dudx(1,:)*(1+v)+(1-u)*dudx(2,:)
END SELECT
wBasis(i,:)=Base*du/n
RotWBasis(i,1)=(dBase(2)*du(3) - dBase(3)*du(2))/n
RotWBasis(i,2)=(dBase(3)*du(1) - dBase(1)*du(3))/n
RotWBasis(i,3)=(dBase(1)*du(2) - dBase(2)*du(1))/n
END SELECT
nj = Element % Nodeindexes(j)
IF (Parallel) nj=Mesh % ParallelInfo % GlobalDOFs(nj)
......@@ -1035,7 +1146,7 @@ CONTAINS
!Numerical integration:
!----------------------
IP = GaussPoints( Element )
IP = GaussPoints(Element)
np = n*Solver % Def_Dofs(Element % BodyId,1)
DO t=1,IP % n
stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
......@@ -3342,10 +3453,10 @@ END SUBROUTINE MagnetoDynamicsCalcFields_Init
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
REAL(KIND=dp) :: s,u,v,w,WBasis(16,3), SOL(2,16), PSOL(16), R(16), C(16), Norm
REAL(KIND=dp) :: RotWBasis(16,3), Basis(16), dBasisdx(16,3), B(2,3), E(2,3)
REAL(KIND=dp) :: s,u,v,w,WBasis(32,3), SOL(2,32), PSOL(32), R(32), C(32), Norm
REAL(KIND=dp) :: RotWBasis(32,3), Basis(32), dBasisdx(32,3), B(2,3), E(2,3)
REAL(KIND=dp) :: detJ, C_ip, R_ip, PR_ip, PR(16), ST(3,3), Omega, Power,Energy
COMPLEX(KIND=dp) :: Magnetization(3,16), MG_ip(3)
COMPLEX(KIND=dp) :: Magnetization(3,32), MG_ip(3)
COMPLEX(KIND=dp) :: CST(3,3)
......@@ -3355,7 +3466,7 @@ END SUBROUTINE MagnetoDynamicsCalcFields_Init
TYPE(Variable_t), POINTER :: Var, MFD, MFS, CD, EF, MST, JH
TYPE(Variable_t), POINTER :: EL_MFD, EL_MFS, EL_CD, EL_EF, EL_MST, EL_JH
INTEGER :: Active, Indexes(12), i,j,k,l,m,n,nd,np,p,q,DOFs,vDOFs,dim
INTEGER :: Active,i,j,k,l,m,n,nd,np,p,q,DOFs,vDOFs,dim
TYPE(Solver_t), POINTER :: pSolver
CHARACTER(LEN=MAX_NAME_LEN) :: Pname
......@@ -3416,7 +3527,7 @@ END SUBROUTINE MagnetoDynamicsCalcFields_Init
ALLOCATE( MASS(n,n), FORCE(n,DOFs) )
SOL = 0._dp; PSOL=0._dp
Omega = GetAngularFrequency()
Omega = GetAngularFrequency(Found=Found)
C = 0._dp; R=0._dp; PR=0._dp
Magnetization = 0._dp
......@@ -3469,8 +3580,8 @@ END SUBROUTINE MagnetoDynamicsCalcFields_Init
dim = CoordinateSystemDimension()
! Calculate on nodes on elements fields:
! --------------------------------------
! Calculate "on nodes on elements" fields:
! ----------------------------------------
DO j=1,n
u = Element % Type % NodeU(j)
v = Element % Type % NodeV(j)
......
......@@ -200,7 +200,7 @@ this ise not in USE
CHARACTER(LEN=MAX_NAME_LEN) :: ViscosityFlag
TYPE(ValueList_t), POINTER :: Material
REAL(KIND=dp) :: x, y, z, c1n(n), c2n(n), c3n(n), c4n(n), &
c1, c2, c3, c4, c5, c6, c7, Temp, h
c1, c2, c3, c4, c5, c6, c7, Temp, TempCoeff, h
! Temperature is needed for thermal models
TYPE(Variable_t), POINTER :: TempSol
......@@ -266,8 +266,10 @@ this ise not in USE
! This is the square of shearrate which results to 1/2 in exponent
! Also the derivative is taken with respect to the square
!-------------------------------------------------------------------
ss = SecondInvariant(Velo,dVelodx,Metric,Symb)/2
SELECT CASE( ViscosityFlag )
CASE('power law')
......@@ -345,41 +347,6 @@ this ise not in USE
c1*LOG(c2*s+SQRT(c2**2*ss+1))/(c2*s**3)/2
END IF
CASE ('thermal carreau')
TempSol => VariableGet( CurrentModel % Variables, 'Temperature' )
IF ( ASSOCIATED( TempSol) ) THEN
TempPerm => TempSol % Perm
Temperature => TempSol % Values
ELSE
CALL Warn('EffectiveViscosity','variable Temperature needed in thermal carreau model')
END IF
c1n = ListGetReal( Material, 'Viscosity Difference',n,Element % NodeIndexes )
c2n = ListGetReal( Material, 'Viscosity Exponent', n, Element % NodeIndexes )
c2 = SUM( Basis * c2n(1:n) )
c3n = ListGetReal( Material, 'Viscosity Transition',n,Element % NodeIndexes )
c4 = ListGetConstReal( Material, 'Yasuda Exponent',gotIt)
c5 = ListGetConstReal( Material, 'Viscosity Temp Ref')
c6 = ListGetConstReal( Material, 'Viscosity Temp Exp')
Temp = SUM(Basis(1:n) * Temperature(TempPerm(Element % NodeIndexes(1:n))))
c7 = EXP(c6*(1/Temp-1/c5))
c1 = c7 * SUM( Basis(1:n) * c1n(1:n) )
c3 = c7 * SUM( Basis(1:n) * c3n(1:n) )
IF(gotIt) THEN
s = SQRT(ss)
mu = Viscosity + c1 * (1.0d0 + (c3*s)**c4)**((c2-1)/c4)
IF ( PRESENT(muder) ) muder = &
c1*(1+(c3*s)**c4)**((c2-1)/c4-1)*(c2-1)*(c3*s)**c4/(2*ss)
ELSE
mu = Viscosity + c1 * (1.0d0 + (c3*c3*ss))**((c2-1)/2.0)
IF ( PRESENT(muder) ) muder = &
c1*(1+c3**2*ss)**((c2-1)/2-1)*(c2-1)/2*c3**2
END IF
CASE( 'smagorinsky' )
c2n = ListGetReal( Material, 'Smagorinsky Constant', &
n, Element % NodeIndexes,gotit )
......@@ -519,6 +486,29 @@ this ise not in USE
CALL WARN('EffectiveViscosity','Unknown material model')
END SELECT
! Add a generic temperature coefficient at the integration point
! for backward compatibility this is activated by an existing keyword
!--------------------------------------------------------------------
c1 = ListGetConstReal( Material, 'Viscosity Temp Exp',GotIt)
IF( GotIt ) THEN
TempSol => VariableGet( CurrentModel % Variables, 'Temperature' )
IF ( ASSOCIATED( TempSol) ) THEN
TempPerm => TempSol % Perm
Temperature => TempSol % Values
ELSE
CALL Warn('EffectiveViscosity','variable Temperature needed for thermal viscosity model')
END IF
c2 = ListGetConstReal( Material, 'Viscosity Temp Ref')
c3 = ListGetConstReal( Material, 'Viscosity Temp Offset',GotIt)
Temp = SUM(Basis(1:n) * Temperature(TempPerm(Element % NodeIndexes(1:n))))
TempCoeff = EXP(c1*(1/(Temp+c3)-1/c2))
mu = TempCoeff * mu
END IF
!------------------------------------------------------------------------------
END FUNCTION EffectiveViscosity
!------------------------------------------------------------------------------
......
......@@ -184,7 +184,7 @@ MODULE NavierStokes
INTEGER :: N_Integ, NBasis, deg(100), Order
REAL(KIND=dp), POINTER :: NodalC(:,:,:)
REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ,V_Integ,W_Integ,S_Integ
TYPE(ValueList_t), POINTER :: Material
TYPE(ValueList_t), POINTER :: Material, Params
LOGICAL :: Found, Transient, stat, Bubbles, PBubbles, Stabilize, &
VMS, P2P1, Isotropic, drhodp_found, Compressible, ViscNewtonLin, &
......@@ -193,8 +193,9 @@ MODULE NavierStokes
!------------------------------------------------------------------------------
dim = CoordinateSystemDimension()
Params => GetSolverParams()
hScale = GetCReal( GetSolverParams(), 'H scale', Found )
hScale = GetCReal( Params, 'H scale', Found )
IF ( .NOT. Found ) hScale = 1._dp
!#ifdef LES
......@@ -239,7 +240,7 @@ MODULE NavierStokes
END IF
IF( ViscNonnewtonian ) THEN
ViscConstantCondition = GetCReal( Material, 'Newtonian Viscosity Condition',Found)
ViscConstantCondition = GetCReal( Params, 'Newtonian Viscosity Condition',Found)
IF( Found .AND. ViscConstantCondition > 0.0_dp ) ViscNonnewtonian = .FALSE.
END IF
......
......@@ -445,6 +445,8 @@ END SUBROUTINE ResultOutputSolver
VarName=ListGetString( Solver % Values, Txt, Found )
IF ( .NOT. Found ) EXIT
Var => VariableGet( Solver % Mesh % Variables, Varname )
IF(.NOT.ASSOCIATED(Var)) &
Var => VariableGet( Solver % Mesh % Variables, TRIM(Varname)//' 1' )
IF ( Var % TYPE == Variable_on_nodes_on_elements ) THEN
On_Nodes_on_elements=On_Nodes_on_elements+1
IF ( PRESENT(ONOEfound) ) Dofs=Dofs+3
......@@ -486,6 +488,8 @@ END SUBROUTINE ResultOutputSolver
WRITE(Txt,'(A,I0)') 'Vector Field ',Vari
VarName = ListGetString( Solver % Values, TRIM(Txt),GotIt)
Var => VariableGet( Solver % Mesh % Variables,VarName )
IF(.NOT.ASSOCIATED(Var)) &
Var => VariableGet( Solver % Mesh % Variables, TRIM(Varname)//' 1' )
IF ( Var % TYPE == Variable_on_nodes_on_elements ) THEN
IF ( .NOT.PRESENT(ONOEfound) ) CYCLE
ELSE
......@@ -809,6 +813,12 @@ END SUBROUTINE ResultOutputSolver
! Check for the presence of component vectors given by its components i=1,2,3
Var => VariableGet( Model % Variables, TRIM(VarName)//' 1' )
IF( ASSOCIATED(Var)) THEN
IF ( Var % TYPE == Variable_on_nodes_on_elements ) THEN
IF (.NOT.PRESENT(ONOEfound) ) CYCLE
ELSE IF ( PRESENT(ONOEfound) ) THEN
CYCLE
END IF
k = i
IF ( ASSOCIATED(Var % Perm) ) k = Var % Perm(k)
......@@ -2060,7 +2070,7 @@ CONTAINS
END IF
IF( dim <= 2 ) THEN
WRITE( IOUnit,'(2ES16.7E3,A)' ) Coord(1:dim),' 0.0'
WRITE( IOUnit,'(2ES16.7E3,A)' ) Coord(1:2),' 0.0'
ELSE
WRITE( IOUnit,'(3ES16.7E3)' ) Coord
END IF
......
......@@ -594,6 +594,7 @@ $ for(i=1:nexp) "Solver:String: 'Operator " i2str(i) "'"
$ for(i=1:nexp) "Solver:String: 'Variable " i2str(i) "'"
$ for(i=1:nexp) "Solver:String: 'Coefficient " i2str(i) "'"
$ for(i=1:nexp) "Solver:String: 'Save Variable " i2str(i) "'"
$ for(i=1:nexp) "Solver:String: 'Parameter " i2str(i) "'"
Solver:Logical: 'File Append'
Solver:Logical: 'Filename Numbering'
Solver:Logical: 'Save Axis'
......
This diff is collapsed.
......@@ -52,30 +52,26 @@
!------------------------------------------------------------------------------
TYPE(ValueList_t), POINTER :: SolverParams
INTEGER :: dim
CHARACTER(LEN=MAX_NAME_LEN) :: VarName, VorName, HelpName
CHARACTER(LEN=MAX_NAME_LEN) :: VarName, FieldName
LOGICAL :: GotIt
!------------------------------------------------------------------------------
SolverParams => GetSolverParams()
dim = CoordinateSystemDimension()
VarName = GetString(SolverParams,'Shearrate Variable',GotIt )
IF(gotIt) THEN
VorName = 'Shearrate '//TRIM(VarName)
ELSE
VorName = 'Shearrate'
END IF
VarName = GetString(SolverParams,'Target Variable',GotIt )
IF(.NOT. gotIt) VarName = TRIM('Velocity')
IF ( .NOT. ListCheckPresent( SolverParams,'Variable') ) THEN
VarName = GetString(SolverParams,'Target Variable',GotIt )
IF(gotIt) THEN
VorName = 'Shearrate '//TRIM(VarName)
ELSE
VorName = 'Shearrate'
END IF
CALL ListAddInteger( SolverParams, 'Variable DOFs', 1 )
CALL ListAddString( SolverParams, 'Variable', VorName )
FieldName = 'Shearrate '//TRIM(VarName)
CALL ListAddString( SolverParams, 'Variable', FieldName )
END IF
IF( GetLogical( SolverParams,'Calculate Viscosity',GotIt ) ) THEN
FieldName = 'Viscosity '//TRIM(VarName)
CALL ListAddString( SolverParams,&
NextFreeKeyword('Exported Variable',SolverParams),FieldName)
END IF
CALL ListAddInteger( SolverParams, 'Time derivative order', 0 )
......@@ -136,12 +132,15 @@ SUBROUTINE ShearrateSolver( Model,Solver,dt,Transient )
CHARACTER(LEN=MAX_NAME_LEN) :: VarName
INTEGER :: i,j,dim,DOFs
LOGICAL :: ConstantBulkMatrix, ConstantBulkMatrixInUse
LOGICAL :: GotIt, Visited = .FALSE.
LOGICAL :: CalculateViscosity, GotIt
REAL(KIND=dp) :: Unorm
REAL(KIND=dp), POINTER :: ForceVector(:), SaveRHS(:)
TYPE(Variable_t), POINTER :: ShearrateSol
SAVE Visited
REAL(KIND=dp), POINTER :: ForceVector(:), ViscVector(:), SaveRHS(:), &
ShearrateField(:), ViscField(:)
TYPE(Variable_t), POINTER :: ShearrateSol, ViscSol
LOGICAL :: AllocationsDone = .FALSE.
SAVE ViscVector, AllocationsDone
CALL Info( 'ShearrateSolver', '-------------------------------------',Level=4 )
......@@ -160,20 +159,32 @@ SUBROUTINE ShearrateSolver( Model,Solver,dt,Transient )
ShearrateSol => Solver % Variable
Dofs = ShearrateSol % DOFs
IF(Dofs /= 1) THEN
PRINT *,'DOFS',Dofs
CALL Fatal('ShearrateSolver','Shearrate should have just 1 component')
END IF
VarName = GetString(SolverParams,'Target Variable',GotIt )
IF(.NOT. gotIt) VarName = TRIM('Velocity')
CalculateViscosity = GetLogical( SolverParams,'Calculate Viscosity',GotIt)
IF( CalculateViscosity ) THEN
IF( .NOT. AllocationsDone ) THEN
ALLOCATE(ViscVector(SIZE(Solver % Matrix % RHS)))
AllocationsDone = .TRUE.
END IF
ViscVector = 0.0_dp
ViscSol => VariableGet( Solver % Mesh % Variables, 'Viscosity '//TRIM(VarName) )
ViscField => ViscSol % Values
END IF
ConstantBulkMatrix = GetLogical( SolverParams, 'Constant Bulk Matrix', GotIt )
ConstantBulkMatrixInUse = ConstantBulkMatrix .AND. &
ASSOCIATED(Solver % Matrix % BulkValues)
ForceVector => Solver % Matrix % rhs
ShearrateField => Solver % Variable % Values
IF ( ConstantBulkMatrixInUse ) THEN
Solver % Matrix % Values = Solver % Matrix % BulkValues
Solver % Matrix % RHS = 0.0_dp
ForceVector = 0.0_dp
ELSE
CALL DefaultInitialize()
END IF
......@@ -182,14 +193,20 @@ SUBROUTINE ShearrateSolver( Model,Solver,dt,Transient )
CALL DefaultFinishAssembly()
!------------------------------------------------------------------------------
IF( CalculateViscosity ) THEN
CALL Info( 'ShearrateSolver','Solving for viscosity',Level=5 )
Solver % Matrix % RHS => ViscVector
Solver % Variable % Values => ViscField
UNorm = DefaultSolve()
Solver % Variable % Values => ShearrateField
Solver % Matrix % RHS => ForceVector
END IF
CALL Info( 'ShearrateSolver','Solving for shearrare',Level=5 )
UNorm = DefaultSolve()
!------------------------------------------------------------------------------
WRITE( Message, * ) 'Result Norm: ',UNorm
CALL Info( 'ShearrateSolver', Message, Level=4 )
CONTAINS
......@@ -198,22 +215,32 @@ CONTAINS
!------------------------------------------------------------------------------
INTEGER :: elem,t,i,j,p,q,n,nd, Rank
REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), FORCE(:)
REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), FORCE(:), FORCE2(:)
TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
TYPE(Nodes_t) :: Nodes
TYPE(Element_t), POINTER :: Element
REAL(KIND=dp) :: SqrtMetric,Metric(3,3),Symb(3,3,3),dSymb(3,3,3,3)
REAL(KIND=dp) :: weight,grad(3),detJ,ShearRate,Velo(3),dVelodx(3,3),s,x,y,z
REAL(KIND=dp) :: weight,grad(3),detJ,ShearRate,Velo(3),dVelodx(3,3),s,x,y,z, &
u,v,w,ViscAtIp, RhoAtIP, mu, muder0
REAL(KIND=dp), ALLOCATABLE :: Basis(:), dBasisdx(:,:)
REAL(KIND=dp), ALLOCATABLE :: Vx(:), Vy(:), Vz(:)
REAL(KIND=dp), ALLOCATABLE :: Vx(:), Vy(:), Vz(:), NodalRho(:), NodalVisc(:)
LOGICAL :: Found
TYPE(ValueList_t), POINTER :: Material
LOGICAL :: AllocationsDone = .FALSE.
SAVE Nodes, STIFF, FORCE, FORCE2, Vx, Vy, Vz, Basis, dBasisdx, &
NodalRho, NodalVisc, AllocationsDone
SAVE Nodes
n = MAX( Solver % Mesh % MaxElementDOFs, Solver % Mesh % MaxElementNodes )
ALLOCATE( STIFF(n,n), FORCE(n) )
ALLOCATE( Vx(n), Vy(n), Vz(n), Basis(n), dBasisdx(n,3) )
IF(.NOT. AllocationsDone ) THEN
n = Solver % Mesh % MaxElementNodes
ALLOCATE( STIFF(n,n), FORCE(n), FORCE2(n), Vx(n), Vy(n), Vz(n), &
Basis(n), NodalVisc(n), NodalRho(n), dBasisdx(n,3) )
Vx = 0.0_dp
Vy = 0.0_dp
Vz = 0.0_dp
AllocationsDone = .TRUE.
END IF
DO elem = 1,Solver % NumberOFActiveElements
......@@ -224,20 +251,34 @@ CONTAINS
nd = GetElementNOFDOFs()
n = GetElementNOFNodes()
IF( CalculateViscosity ) THEN
Material => GetMaterial()
NodalVisc(1:n) = GetReal( Material,'Viscosity')
NodalRho(1:n) = GetReal( Material,'Density')
END IF
CALL GetScalarLocalSolution( Vx, ComponentName(VarName,1) )
CALL GetScalarLocalSolution( Vy, ComponentName(VarName,2) )
CALL GetScalarLocalSolution( Vz, ComponentName(VarName,3) )
CALL GetScalarLocalSolution( Vy, ComponentName(VarName,2) )
IF( dim > 2 ) THEN
CALL GetScalarLocalSolution( Vz, ComponentName(VarName,3) )
END IF
! Integrate local stresses:
! -------------------------
IntegStuff = GaussPoints( Element )
STIFF = 0.0_dp
FORCE = 0.0_dp
FORCE2 = 0.0_dp
DO t=1,IntegStuff % n
Found = ElementInfo( Element, Nodes, IntegStuff % u(t), &
IntegStuff % v(t), IntegStuff % w(t), detJ, Basis, dBasisdx )
u = IntegStuff % u(t)
v = IntegStuff % v(t)
w = IntegStuff % w(t)
Found = ElementInfo( Element, Nodes, u, v, w, detJ, Basis, dBasisdx )
!------------------------------------------------------------------------------
! Coordinate system dependent information
!------------------------------------------------------------------------------
......@@ -248,6 +289,17 @@ CONTAINS
s = detJ * SqrtMetric * IntegStuff % s(t)
IF( CalculateViscosity ) THEN
RhoAtIP = SUM( Nodalrho(1:n) * Basis(1:n) )
ViscAtIP = SUM( NodalVisc(1:n) * Basis(1:n) )
IF( ListCheckPresent( Material, 'Viscosity Model' ) ) THEN
mu = EffectiveViscosity( ViscAtIP, RhoAtIp, Vx, Vy, Vz, &
Element, Nodes, n, n, u, v, w, muder0 )
ELSE
mu = ViscAtIP
END IF
END IF
IF ( .NOT. ConstantBulkMatrixInUse ) THEN
DO p=1,nd
DO q=1,nd
......@@ -269,6 +321,10 @@ CONTAINS
ShearRate = SQRT( SecondInvariant(Velo,dVelodx,Metric,Symb)/2 )
FORCE(1:nd) = FORCE(1:nd) + Basis(1:nd) * s * ShearRate
IF( CalculateViscosity ) THEN
FORCE2(1:nd) = FORCE2(1:nd) + Basis(1:nd) * s * mu
END IF
END DO
!------------------------------------------------------------------------------
......@@ -281,9 +337,16 @@ CONTAINS
ELSE
CALL DefaultUpdateForce( FORCE(1:nd) )
END IF
IF( CalculateViscosity ) THEN
Solver % Matrix % RHS => ViscVector
CALL DefaultUpdateForce( FORCE2 )
Solver % Matrix % RHS => ForceVector
END IF
END DO
DEALLOCATE( STIFF, FORCE, Basis, dBasisdx, Vx, Vy, Vz )
!------------------------------------------------------------------------------
END SUBROUTINE BulkAssembly
......
! Ideal 1D flow obtained by periodic BCs for the pressure
! This enables testing of nonnewtonian velocity profiles
! For the first coupled system iteration newtonian model is enforced.
! For the first nonlinear system iteration newtonian model is enforced.
$ p0 = 0.0
$ v0 = 1.0
......@@ -16,7 +16,7 @@ Simulation
Max Output Level = 4
Coordinate System = Cartesian
Simulation Type = Steady State
Steady State Max Iterations = 50
Steady State Max Iterations = 1
Output Intervals = 0
Post File = case.ep
Set Dirichlet BCs by BC Numbering = True
......@@ -46,8 +46,10 @@ Solver 1
Steady State Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-6
Nonlinear System Max Iterations = 1
Nonlinear System Newton After Iterations = 50
Nonlinear System Max Iterations = 50
! activates also the newtonian linearization of viscosity
Nonlinear System Newton After Iterations = 5
Nonlinear System Relaxation Factor = 1
Nonlinear System Newton After Tolerance = 1.0e-3
......@@ -61,6 +63,13 @@ Solver 1
Linear System Convergence Tolerance = 1.0e-10
Linear System Residual Output = 10
Linear System Preconditioning = ILU0
! If this is more than zero the newtonian fluid prevails
Newtonian Viscosity Condition = Variable nonlin iter
Real
1.0 1.0
2.0 -1.0
End
End