From 836547a000aaafedc4c7bbe2eb68d06e56682cca Mon Sep 17 00:00:00 2001 From: unknown Date: Fri, 17 Apr 2020 00:25:36 +0200 Subject: [PATCH] version 6 --- changes.txt | 74 ++- macros/all.OCM | 8 +- macros/step1.OCM | 4 +- minimizer/matsmin.F90 | 112 ++-- models/gtp3.F90 | 1 + models/gtp3A.F90 | 118 ++++- models/gtp3C.F90 | 8 + models/gtp3H.F90 | 29 +- models/gtp3Y.F90 | 84 ++- models/gtp3Z.F90 | 114 ++-- pmain1.F90 | 2 +- stepmapplot/smp2.F90 | 8 +- stepmapplot/smp2A.F90 | 1150 +++++++++++++++++++++++++++++++---------- stepmapplot/smp2B.F90 | 28 +- userif/pmon6.F90 | 27 +- 15 files changed, 1362 insertions(+), 405 deletions(-) diff --git a/changes.txt b/changes.txt index 07310c5..c002b83 100644 --- a/changes.txt +++ b/changes.txt @@ -1,38 +1,68 @@ -This file contains information on OpenCalphad (OC) version 5 and earlier +This file contains information on OpenCalphad (OC) version 6 and earlier -For general information please read the readme-general.pdf +For general information please read the intro-OC6.pdf and the +OC6-macros.pdf. -The installation and use of OC requires some general knowledge about -compilation and linking of software. If you are not familiar with -such procedures please ask a local guru. We who are providing this -software do not have time to answer such questions. - -If you want a thermodynamic software that can install itself and which -does not require any understanding of thermodynamics please contact a -commercial vendor. +The installation of OC on Windows is now fully automatic including an +introduction, on-line help and macro examples. On Linux, MacOS and +other OS the installation requires some general knowledge about +compiling and linking of software. There is an installation manual +(and video on You Tube) but if you are not familiar with such +procedures please ask a local guru for help. We who are providing +this software do not have time to answer such questions. In the list below the most recent changes come first. +************************ +* This is OC version 6 * +************************ + +2020.04.14 OC version 6.001: The Coronavirus epidemy is changing the +world as we know it. + +An automatic installation of OC6 for Windows is now available on the +opencalphad.org website. This creates the directories and the OCHOME +environment variable and installs the executable, on-line help files, +macro files, the macro manual and an intoductary text. Many thanks to +Chunhui for the help with his. The macro manual and the informative +text has been updated. + +The calculation of invariant equilibria for isopleths was added in the +previous version but the handling of the many exit lines was not +elaborated. This has been improved, for example in map7 and map16, +but there are still some problems to find the set of fix phases +defining each line exiting rom the invariant. + +I found that "c n" actually called the grid minimizer after the +iterative calculation if the conditions did not allow using the grid +minimizer initially. Now "c n" never call the grid minimizer. + +I have also fixed that the STEP command can handle invariant +equilibria, previously it stopped there because it is necessary to add +and remove a phase at the same node. Now it should discover by itself +which phases to exchange. As always step and map are fragile and may +require several startpoints to be complete. + ************************************************************************ -* This is a prerelease of OC version 6 * +* This is a prerelease of OC version 6 * ************************************************************************ 2020.03.15 OC version 6.000: Note the linkpara and Makefiles are changed to use the name oc6p for the executable. I will update the "stable" version on http://www.opencalphad.org soon. The handling of isopleth diagrams has improved and now I have added colors to the -lines, indicating which phase which haz zero amount, i.e. becomes -stable along the line. The isopleth invariants are now calculated -correctly but there are still problems generating the correct set of -stable phases at the exit lines from the invariants. There is a new -macro, map16, featuring an isothermal invariant in the C-Cr-Fe system. -Macros 6 and 7 also features the colors but there are several -incmplete or missing lines. - -I will add automatic start points for phase diagram mapping but one +lines, indicating which phase has zero amount, i.e. becomes stable +along the line. The isopleth invariants are now calculated correctly +but there are still problems generating the correct set of stable +phases at the exit lines from the invariants. There is a new macro, +map16, featuring an isothermal invariant in the C-Cr-Fe system. +Macros 6 and 7 also features the colors but there are some incomplete +or missing lines. + +I intend add automatic start points for phase diagram mapping but one problem is that sometimes metastable lines are calculated. There are -also numerical problems when following a line and some line may end in -the middle of nowhere. A complete diagram may require several +also numerical problems when following a line and some lines may end +in the middle of nowhere. A complete diagram may require several independent start points but such a routine may generate lines representing metastable equilibria. diff --git a/macros/all.OCM b/macros/all.OCM index 2afe6e6..41d8229 100644 --- a/macros/all.OCM +++ b/macros/all.OCM @@ -183,16 +183,16 @@ mac ./map16 mac ./map17 @$ ********************************************************* -@$ Calculation for 20 elements and 191 phases using COST507 +@$ Testing the UNIQUAC model @& ********************************************************* -mac ./allcost +mac ./uniquac @$ ********************************************************* -@$ Testing the UNIQUAC model +@$ Calculation for 20 elements and 191 phases using COST507 @& ********************************************************* -mac ./uniquac +mac ./allcost @$ ********************************************************* @$ Calculating 21 equilibria in parallel diff --git a/macros/step1.OCM b/macros/step1.OCM index a937643..fcbcb7c 100644 --- a/macros/step1.OCM +++ b/macros/step1.OCM @@ -32,7 +32,9 @@ set c t=1200 p=1e5 n=1 w(c)=.009 w(cr)=.045, w(mo)=.1,w(si)=.001 w(v)=.009 @$ Enter a composition set for the MC carbide (FCC) @$ This is convenient to specify an additional pre/suffix -amend phase fcc comp_set y MC , +amend phase fcc comp_set y +MC + NONE <.1 NONE diff --git a/minimizer/matsmin.F90 b/minimizer/matsmin.F90 index 5e9f2dd..69e5146 100644 --- a/minimizer/matsmin.F90 +++ b/minimizer/matsmin.F90 @@ -51,7 +51,11 @@ MODULE liboceqplus ! BITS in meqrec status word ! MMQUIET means no output for the equilibrium calculation ! MMNOSTARTVAL means grid minimizer not called at start - integer, parameter :: MMQUIET=0, MMNOSTARTVAL=1 + integer, parameter :: MMQUIET=0, MMNOSTARTVAL=1,MMSTEPINV=2 +! NOTE in calceq7 status word is set to zero if more bits used because +! it seemed to have an arbitrary value and it created problems in macro map7 +! I have now correceted the main reason (creating linehead records in SMP) +! but I kept this check ! !\begin{verbatim} TYPE meq_phase @@ -146,7 +150,7 @@ MODULE liboceqplus !\begin{verbatim} TYPE map_fixph ! provides information about phase sets for each line during mapping - integer nfixph,nstabph + integer nfixph,nstabph,status type(gtp_phasetuple), dimension(:), allocatable :: fixph type(gtp_phasetuple), dimension(:), allocatable :: stableph ! most likely some of these variables are redundant stable_phr added 2020.03.05 @@ -156,6 +160,7 @@ MODULE liboceqplus double precision, dimension(:), allocatable :: fixphamap end TYPE map_fixph !\end{verbatim} +! declared as mapfix in call to calceq7 and some other routines ! ! Added for debugging converge problems TYPE meqdebug @@ -254,7 +259,8 @@ subroutine calceq2(mode,ceq) ! For example to ensure a fcc-carbonitrides is always the same compset. ij=1 ! if meqrec%status indicate no initial startvalues set ij<0 to indicate test - if(btest(meqrec%status,MMNOSTARTVAL)) ij=-ij +! DO not test if mode=0 + if(mode.ne.0 .and. btest(meqrec%status,MMNOSTARTVAL)) ij=-ij ! OC went into a loop for a complex alloy calcumation here (once long ago ...) ! write(*,*)'MM calling todo_after: 2',& ! btest(meqrec%status,MMNOSTARTVAL),mode @@ -322,7 +328,7 @@ subroutine calceq3(mode,confirm,ceq) ! For example to ensure a fcc-carbonitrides is always the same compset. ij=1 ! if meqrec%status indicate no initial startvalues set ij<0 to indicate test - if(btest(meqrec%status,MMNOSTARTVAL)) ij=-ij + if(mode.ne.0 .and. btest(meqrec%status,MMNOSTARTVAL)) ij=-ij ! write(*,*)'MM Calling todo_after calceq3' call todo_after_found_equilibrium(ij,addtuple,ceq) if(gx%bmperr.eq.4358) then @@ -393,7 +399,14 @@ subroutine calceq7(mode,meqrec,mapfix,ceq) ntup=nooftup() ! write(*,*)'MM in calceq7',ntup ycond=.FALSE. - meqrec%status=0 +! if(btest(meqrec%status,MMSTEPINV)) then +! this is the problem with map7? only bit 0 and 1 are used!! +! write(*,'(a,z8)')'MM warning **** eqcalc7 meqrec%status: ',& +! meqrec%status,' reset' +! meqrec%status=0 +! meqrec%status may be set by STEPMAPPLOT to indicate +! write(*,*)'MM calceq7 meqrec%status: ',meqrec%status,meqrec%nstph +! endif if(btest(globaldata%status,GSSILENT)) & meqrec%status=ibset(meqrec%status,MMQUIET) if(ocv()) write(*,*)"Entering calceq7",mode @@ -442,7 +455,8 @@ subroutine calceq7(mode,meqrec,mapfix,ceq) ! gx%bmperr.eq.4174 .or. & ! (gx%bmperr.ge.4176 .and. gx%bmperr.le.4185)) goto 1000 ! if mode=0 we should not use grid minimizer - if(mode.ne.0 .or. .not.btest(meqrec%status,MMQUIET)) & +! if(mode.ne.0 .or. .not.btest(meqrec%status,MMQUIET)) & + if(mode.ne.0 .and. .not.btest(meqrec%status,MMQUIET)) & write(*,9) 9 format('Warning: global minimizer cannot be used for the current',& ' set of conditions') @@ -867,8 +881,14 @@ subroutine calceq7(mode,meqrec,mapfix,ceq) meqrec%icsl(meqrec%nv)=meqrec%fixph(2,mjj) meqrec%aphl(meqrec%nv)=meqrec%fixpham(mjj) enddo addfixph -!------------------------------- special for mapping - if(allocated(mapfix)) then +!------------------------------- special for mapping and STEP + mapfixdata: if(allocated(mapfix)) then +! for step only the status word is used to indicate an invarant node +! if(mapfix%nfixph.eq.0) then +! if(btest(mapfix,STEPINVARIANT)) then +! exit mapfixdata +! endif +! endif ! the stable and fix phases copied from mapfix record. do ij=1,meqrec%nv meqrec%iphl(ij)=0 @@ -909,7 +929,19 @@ subroutine calceq7(mode,meqrec,mapfix,ceq) ! mapfix%stableph(1)%ixphase,mapfix%stableph(1)%compset,& ! mapfix%stablepham(1) 64 format(a,i3,2i5,1pe12.4) - endif +! elseif(formap) then +! mapfixrecord not allocated for STEP calculations +! this dis not work for handling invariant nodes for STEP +! write(*,*)'MM calceq7 formap MMSTEPINV:',btest(meqrec%status,MMSTEPINV) +! if(btest(meqrec%status,MMSTEPINV)) then +! The line start at an invariant node for a STEP calculation, +! write(*,*)'MM invariant node with phases: ',meqrec%nstph +! do jj=1,meqrec%nstph +! jq=meqrec%stphl(jj) +! write(*,*)'MM stable: ',jj,jq,meqrec%phr(jq)%curd%amfu +! enddo +! endif + endif mapfixdata !------------------------------- ! zero start of link to phases set temporarily dormant .... meqrec%dormlink=0 @@ -952,10 +984,10 @@ subroutine calceq7(mode,meqrec,mapfix,ceq) call sumprops(props,ceq) if(gx%bmperr.ne.0) then write(*,*)'Convergence error, check your conditions are reasonable' - elseif(props(4).gt.10 .and. & + elseif(props(4).gt.1.0D1 .and. & .not.(saverr.eq.4210 .or. saverr.eq.4364)) then - write(*,*)'Convergence error: *** REDUCE THE SIZE OF YOUR SYSTEM!',& - saverr + write(*,'(a,a,i5,1pe12.4)')'Convergence error, maybe reduce ',& + 'the size of your system!',saverr,props(4) endif gx%bmperr=saverr endif @@ -1071,6 +1103,7 @@ subroutine meq_phaseset(meqrec,formap,mapfix,ceq) ! new: -4 hidden, -3 suspended, -2 dormant, -1,0,1 entered, 2 fixed if(zap.ge.PHDORM) then mph=mph+1 +! this iph is the index in the phlista record meqrec%phr(mph)%iph=iph meqrec%phr(mph)%ics=ics ! compare with these the first time a phase wants to be added or removed @@ -1272,7 +1305,8 @@ subroutine meq_phaseset(meqrec,formap,mapfix,ceq) ! addremloop(iadd,2),addremloop(iadd,3) endif if(addremloop(iadd,2).gt.5) then - write(*,'(a,2i4,"#",i1)')'MM Removing phase: ',iadd,& + if(.not.btest(meqrec%status,MMQUIET)) & + write(*,'(a,2i4,"#",i1)')'MM Removing phase: ',iadd,& meqrec%phr(iadd)%iph,meqrec%phr(iadd)%ics meqrec%phr(iadd)%phasestatus=PHDORM meqrec%phr(iadd)%curd%phstate=PHDORM @@ -1361,7 +1395,29 @@ subroutine meq_phaseset(meqrec,formap,mapfix,ceq) phnames='-' call get_phasetup_name(tuprem,phnames(2:)) endif - if(formap) then + addph: if(formap) then +! if(btest(meqrec%status,MMSTEPINV)) then +! This did not work to handle invariants during STEP +! we are exiting an invariant node for a STEP calculation, allow phase change +! meq_sameset wants to ADD a phase, instead remove the last stable phase +! write(*,*)'MM meq_phaseset invariant node',meqrec%noofits,iadd +! do jj=1,meqrec%nstph +! irem=meqrec%stphl(jj) +! if(iadd.eq.0 .and. & +! meqrec%phr(irem)%curd%amfu.eq.zero) then +! meqrec%phr(irem)%curd%amfu=1.0D-1 +! endif +! write(*,*)'MM stable: ',jj,irem,meqrec%phr(irem)%curd%amfu +! enddo +! if(iadd.gt.0 .and. meqrec%nstph.gt.1) then +! meqrec%nstph=meqrec%nstph-1 +! meqrec%phr(irem)%curd%amfu=zero +! write(*,*)'MM ignore adding ',iadd,' but remove ',irem +! iadd=0 +! goto 200 +! endif +! exit addph +! endif ! This can be too strong, we can have a tie-line betwen two stoichiometric ! phases, i.e. a new phase appears at first attempt to step in two-phase region. ! UNFINISHED handling of many exceptions during mapping @@ -1395,15 +1451,8 @@ subroutine meq_phaseset(meqrec,formap,mapfix,ceq) write(*,281)meqrec%noofits,trim(phnames) 281 format('Phase change iteration: ',i5,2x,a) #endif - endif + endif addph endif -! if(formap) then -! when called during mapping the set of phases must not change! -! if(ocv()) write(*,*)'Phase change not allowed',ceq%tpval(1) -! Phase change not allowed due to step/map constraints -! step/map should handle this by creating a node point -! gx%bmperr=4210; goto 1000 -! endif endif 222 continue remove: if(irem.gt.0) then @@ -1662,11 +1711,13 @@ subroutine meq_phaseset(meqrec,formap,mapfix,ceq) write(*,*)'MM cannot find phasetup name: ',jj,kk,gx%bmperr gx%bmperr=0 endif - if(meqrec%phr(jj)%curd%dgm.gt.zero) then - write(*,1220)jj,kk,trim(phnames),meqrec%phr(jj)%curd%dgm -1220 format('Restoring phase: ',2i5,2x,a,5x,1pe12.4) - else - write(*,1220)jj,kk,trim(phnames) + if(.not.btest(meqrec%status,MMQUIET)) then + if(meqrec%phr(jj)%curd%dgm.gt.zero) then + write(*,1220)jj,kk,trim(phnames),meqrec%phr(jj)%curd%dgm +1220 format('MM Restoring phase: ',2i5,2x,a,5x,1pe12.4) + else + write(*,1220)jj,kk,trim(phnames) + endif endif if(meqrec%phr(jj)%curd%dgm.gt.1.0D-2) jph=jj ! do I have two places for suspendeded ?? YES!! @@ -1677,7 +1728,8 @@ subroutine meq_phaseset(meqrec,formap,mapfix,ceq) goto 1200 endif if(jph.gt.0) then - write(*,*)' *** Warning, a restored phase wants to be stable:',jph + if(.not.btest(meqrec%status,MMQUIET)) & + write(*,*)'MM warning, a restored phase wants to be stable:',jph gx%bmperr=4363 endif ! we may already have had an error ... @@ -7685,6 +7737,8 @@ subroutine meq_state_var_dot_derivative(svr1,svr2,value,ceq) write(*,*)'MM Allocation error 48: ',errall gx%bmperr=4370; goto 1000 endif +! problems with map7 ??? + meqrec1%status=0 meqrec=>meqrec1 ! write(*,88)'MM calling initiate_meqrec',svr2%statevarid,ceq%eqno 88 format(a,2i4) @@ -10378,7 +10432,7 @@ recursive subroutine two_stoich_same_comp(irem,iadd,mapx,meqrec,inmap,ceq) ! How to check if I should use this routine? Only 2 components? ! If we have an activity condition one could have 3 components .... write(*,*)'MM This routine should be used only when tie-lines in plane' - goto 1000 + gx%bmperr=4399; goto 1000 endif ! call get_state_var_value('X(O) ',value,phases,ceq) ! write(*,806)meqrec%fixph(1,1),meqrec%fixph(2,1),mapx,iadd,value diff --git a/models/gtp3.F90 b/models/gtp3.F90 index 0691fbc..49bae6f 100644 --- a/models/gtp3.F90 +++ b/models/gtp3.F90 @@ -929,6 +929,7 @@ MODULE GENERAL_THERMODYNAMIC_PACKAGE ! TPVALUE set if evaluated only explicitly (keeping its value) ! TPEXPORT set if value should be exported to symbol ! TPIMPORT set if value should be imported from symbol (only for constants) +! TPINTEIN set if value should always be calculated integer, parameter :: & TPCONST=0, TPOPTCON=1, TPNOTENT=2, TPVALUE=3, & TPEXPORT=4, TPIMPORT=5 diff --git a/models/gtp3A.F90 b/models/gtp3A.F90 index 7ad1a03..823044b 100644 --- a/models/gtp3A.F90 +++ b/models/gtp3A.F90 @@ -1783,7 +1783,7 @@ end function mass_of subroutine get_phase_name(iph,ics,name) ! Given the phase index and composition set number this subroutine returns ! the name with pre- and suffix for composition sets added and also -! a \# followed by a digit 2-9 for composition sets higher than 1. +! a \# followed by a digit 1-9 if there are more than one composition sets implicit none character name*(*) integer iph,ics @@ -1794,6 +1794,12 @@ subroutine get_phase_name(iph,ics,name) if(gx%bmperr.ne.0) goto 1000 if(ics.eq.1) then name=phlista(lokph)%name + if(phlista(lokph)%noofcs.ge.2) then +! this was added 2020.04.02 because a call to change_many_phase_status +! using a phase name returned from this routine will suspend all compsets + kp=len_trim(name)+1 + name(kp:)='#1' + endif else kp=len_trim(firsteq%phase_varres(lokcs)%prefix) if(kp.gt.0) then @@ -1835,26 +1841,44 @@ end subroutine get_phasetup_name !/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ -!\addtotable subroutine get_phasetup_name_old -!\begin{verbatim} - subroutine get_phasetup_name_old(phtuple,name) -! Given the phase tuple this subroutine returns the name with pre- and suffix -! for composition sets added and also a \# followed by a digit 2-9 for -! composition sets higher than 1. +!\addtotable subroutine get_phasetuple_name +!\begin{verbatim} %- + subroutine get_phasetuple_name(phtuple,name) +! phtuple is a phase tuple +! the name has pre- and suffix for composition sets added and also +! a \# followed by a digit 2-9 for composition sets higher than 1. implicit none character name*(*) type(gtp_phasetuple) :: phtuple !\end{verbatim} %+ -! -! PROBABLY REDUNDANT and wrong ... -! -! call get_phase_name(phtuple%phaseix,phtuple%compset,name) +! integer phx,phy call get_phase_name(phtuple%ixphase,phtuple%compset,name) 1000 continue return - end subroutine get_phasetup_name_old + end subroutine get_phasetuple_name !/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ +! +!-\addtotable subroutine get_phasetup_name_old +!-\begin{verbatim} +! subroutine get_phasetup_name_old(phtuple,name) +! Given the phase tuple this subroutine returns the name with pre- and suffix +! for composition sets added and also a \# followed by a digit 2-9 for +! composition sets higher than 1. +! implicit none +! character name*(*) +! type(gtp_phasetuple) :: phtuple +!-\end{verbatim} %+ +! +! PROBABLY REDUNDANT and wrong ... +! +! call get_phase_name(phtuple%phaseix,phtuple%compset,name) +! call get_phase_name(phtuple%ixphase,phtuple%compset,name) +!1000 continue +! return +! end subroutine get_phasetup_name_old +! +!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ !\addtotable subroutine get_phasetup_record !\begin{verbatim} %- @@ -2694,5 +2718,75 @@ integer function gettupix(iph,ics) return end function gettupix +!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ + +!\addtotable subroutine suspend_somephases +!\begin{verbatim} + subroutine suspend_somephases(mode,invph,dim1,dim2,ceq) +! This was added to handle calculating restricted equilibria during mapping +! to suspend (mode=1) or restore (mode=0) phases not involved +! in an invariant equilibrium. +! invph is array with phases that are involved, it has dimension (dim1,*) +! the current status is saved and restored + implicit none + integer mode,dim1,dim2,invph(dim1,*) + type(gtp_equilibrium_data), pointer :: ceq +!\end{verbatim} + integer, save, allocatable, dimension(:) :: phtupixstatus + integer, save :: ntup + integer ii,jj,kk,lokcs,lokph + character phname*24 + ii=nooftup() + kk=0 + if(mode.eq.1) then +! after saving current status suspend all phases not included in invph +! write(*,*)'3A suspending some phases',ii + ntup=ii + if(allocated(phtupixstatus)) then + write(*,*)'3A calls to suspend_somephases cannot be nested' + gx%bmperr=4399; goto 1000 + else + allocate(phtupixstatus(ntup)) + endif + loop1: do ii=1,ntup + lokcs=phasetuple(ii)%lokvares + phtupixstatus(ii)=ceq%phase_varres(lokcs)%phstate + do jj=1,dim2 +! write(*,*)'3A suspend? ',jj,lokcs,& +! phlista(invph(1,jj))%linktocs(invph(2,jj)),phtupixstatus(ii) +! invph(1,jj) is index in phases (phase and alphabetcal order) +! lokph is the order the phase were entered into phlista (arbitrary) + lokph=phases(invph(1,jj)) + if(lokcs.eq.phlista(lokph)%linktocs(invph(2,jj))) then +! write(*,'(a,6i5)')'3A not suspending',jj,invph(1,jj),& +! invph(2,jj),phlista(lokph)%linktocs(invph(2,jj)) + cycle loop1 + endif + enddo +! this phase should be suspended + kk=kk+1 + ceq%phase_varres(lokcs)%phstate=PHSUS + enddo loop1 +! write(*,'(a,i3,a,i3)')'3A suspededed ',kk,' phases out of ',ntup + elseif(mode.eq.0) then +! restore status of all phases except those in invph +! write(*,*)'3A restoring some phases',ii + if(ii.ne.ntup) then + write(*,*)'3A number of phases and compsets changed!',ntup,ii + stop + endif + do ii=1,ntup + ceq%phase_varres(phasetuple(ii)%lokvares)%phstate=phtupixstatus(ii) + enddo +! write(*,'(a,i3,a)')'3A restored phase status for ',ntup,' phases' + deallocate(phtupixstatus) + else + write(*,*)'3A mode must be 0 or 1' + gx%bmperr=4399 + endif +1000 continue + return + end subroutine suspend_somephases + !/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ diff --git a/models/gtp3C.F90 b/models/gtp3C.F90 index 36c179c..e8a32f8 100644 --- a/models/gtp3C.F90 +++ b/models/gtp3C.F90 @@ -331,6 +331,10 @@ subroutine list_sorted_phases(unit,ceq) ! skip composition set number and pre/suffixes at present .... susph(nsusp:)=trim(phlista(lokph)%name)//', ' nsusp=len_trim(susph)+2 + if(ics.gt.1) then + susph(nsusp-2:)='#'//char(ichar('0')+ics)//',' + nsusp=nsusp+2 + endif endif enddo csloop enddo phloop @@ -499,6 +503,10 @@ subroutine list_all_phases(unit,ceq) susph(nsusp:)=phlista(lokph)%name(1:& len_trim(phlista(lokph)%name))//', ' nsusp=len_trim(susph)+2 + if(ics.gt.1) then + susph(nsusp-2:)='#'//char(ichar('0')+ics)//',' + nsusp=nsusp+2 + endif cycle endif elseif(csrec%phstate.ne.PHDORM) then diff --git a/models/gtp3H.F90 b/models/gtp3H.F90 index d92f26e..a5793cf 100644 --- a/models/gtp3H.F90 +++ b/models/gtp3H.F90 @@ -1453,9 +1453,10 @@ subroutine calc_einsteincp(moded,phres,addrec,lokph,mc,ceq) ! NOTE ALL VALUES CALCULATED AS FOR G/RT ! kvot=theta/T if(phres%gval(1,ith).gt.1.0D2) then - write(*,*)'Most likely wrong value of THET, parameter should be ln(THET)' + write(*,*)'3H Probably wrong value of THET, parameter should be ln(THET)' gx%bmperr=4399; goto 1000 endif +! The exp( ) because the parameter value is LN(THETA) kvot=exp(phres%gval(1,ith))/ceq%tpval(1) ! write(*,70)'3H phres: ',ceq%tpval(1),phres%gval(1,1),phres%gval(2,1),& ! phres%gval(3,1),phres%gval(4,1),kvot @@ -1487,7 +1488,7 @@ subroutine calc_einsteincp(moded,phres,addrec,lokph,mc,ceq) kvotexpkvotm1=kvot/(exp(kvot)-one) ln1mexpkvot=log(one-expmkvot) endif -! +! kvot is +THETA/T; gein is integrated contribution to the Gibbs energy gein=1.5D0*kvot+3.0D0*ln1mexpkvot ! write(*,71)'3H Cp E1:',extreme,ceq%tpval(1),gein,ln1mexpkvot,expmkvot,& ! kvotexpkvotm1 @@ -1993,9 +1994,12 @@ subroutine calc_twostate_model_john(moded,phres,addrec,lokph,mc,ceq) ! write(*,*)'3H found liquid delta-cp: ',dcpl endif enddo findix +! ignore this error message if(ith.eq.0) then - write(*,*)'Cannot find value for amorphous THET' - gein=zero; dgeindt=zero; d2geindt2=zero; goto 300 +! write(*,*)'3H Cannot find value for amorphous THET' + write(*,*)'3H warning no value for amorphous THET' + gein=zero; dgeindt=zero; d2geindt2=zero + goto 300 ! gx%bmperr=4367; goto 1000 endif !---------------------------------- @@ -2236,11 +2240,14 @@ subroutine calc_twostate_model1(moded,phres,addrec,lokph,mc,ceq) endif enddo findix if(ith.eq.0) then - write(*,*)'Cannot find value for amorphous THET' - gx%bmperr=4367; goto 1000 +! write(*,*)'3H Cannot find value for amorphous THET' + write(*,*)'3H warning: no value for amorphous THET' + gein=zero; dgeindt=zero; d2geindt2=zero + goto 300 +! gx%bmperr=4367; goto 1000 endif if(ig2.eq.0) then - write(*,*)'Cannot find value for G2 two-state parameter' + write(*,*)'3H Cannot find value for G2 two-state parameter' gx%bmperr=4367; goto 1000 endif !---------------------------------- @@ -2297,7 +2304,9 @@ subroutine calc_twostate_model1(moded,phres,addrec,lokph,mc,ceq) else d2geindt2=-3.0D0*kvotexpkvotm1**2/(expmkvot*ceq%tpval(1)**2) endif +!-------------------------- jump here if no THET variable ! return the values in phres%gval(*,1) +300 continue phres%gval(1,1)=phres%gval(1,1)+msize*gein phres%gval(2,1)=phres%gval(2,1)+msize*dgeindt ! phres%gval(3,1)=phres%gval(3,1) @@ -2333,7 +2342,7 @@ subroutine calc_twostate_model1(moded,phres,addrec,lokph,mc,ceq) ! if g2val is positive we are in the amorphous region ! if g2val is negative we are in the liquid region ! The if statements here ensure expmg2 is between 1e-60 and 1e+60 -! write(*,'(a,6(1pe12.4))')'3H g2val: ',g2val,dg2dt,-g2val/rt + write(*,'(a,6(1pe12.4))')'3H g2val: ',g2val,dg2dt,-g2val/rt if(-g2val/rt.gt.2.0D2) then ! LIQUID REGION exp(200) >> 1, thus d2g=ln(1+exp(g2val))=g2val ! and the derivatives are those above. DIVIDED BY RT? @@ -2378,6 +2387,10 @@ subroutine calc_twostate_model1(moded,phres,addrec,lokph,mc,ceq) ((g2val/tv)**2+(dg2dt)**2-2.0D0*(g2val/tv)*dg2dt)*expg2/& (rt*(one+expg2)**2))/rt 700 continue + write(*,705)'3H 2SL: ',g2val/rt, dg2, dgfdt, dgfdt, d2g2dt2, tv,& + rt, expg2, dg2dt, msize, d2g2dt2*rt +705 format(a,6(1pe12.4)/8x,6(1pe12.4)) +! THIS IS THE SUBROUTINE USED FOR 2STATE LIQUID ! This should be OK/ 2020.02.27 phres%gval(1,1)=phres%gval(1,1)-msize*dg2 phres%gval(2,1)=phres%gval(2,1)-msize*dgfdt diff --git a/models/gtp3Y.F90 b/models/gtp3Y.F90 index cd45985..bdd127b 100644 --- a/models/gtp3Y.F90 +++ b/models/gtp3Y.F90 @@ -5022,7 +5022,7 @@ logical function global_equil_check1(mode,addtuple,yfr,ceq) double precision totmol,totmass,amount,gmax,dgmax,dgtest integer, allocatable :: kphl(:),iphx(:) integer gmode,iph,ngg,nrel,ny,ifri,firstpoint,sumng,nrph,ii,jj,nz,lokcs - integer ics,pph,nyfas,gpz,iphz,nggz,errall + integer ics,pph,nyfas,gpz,iphz,nggz,errall,haha integer, parameter :: maxgrid=400000 ! ! write(*,*)'3Y In global_equil_check1',mode @@ -5037,6 +5037,12 @@ logical function global_equil_check1(mode,addtuple,yfr,ceq) if(mode.ne.1) addgridpoint=.FALSE. dgmax=zero addtuple=0 +! Problem with invariant when mapping but not here +! if(inveq(haha,ceq)) then +! write(*,*)'3Y equilibrium is invariant when entering',haha +! else +! write(*,*)'3Y equilibrium is not invariant when entering',haha +! endif ! COPY the whole equilibrium record to avoid destroying anything!! ! otherwise I had strange problems with amounts of phases ?? cceq=ceq @@ -5227,6 +5233,11 @@ logical function global_equil_check1(mode,addtuple,yfr,ceq) deallocate(iphx) endif 2000 continue +! if(inveq(haha,ceq)) then +! write(*,*)'3Y equilibrium is invariant when at exit',haha +! else +! write(*,*)'3Y equilibrium is not invariant when at exit',haha +! endif global_equil_check1=global return end function global_equil_check1 @@ -7117,5 +7128,76 @@ integer function vssize(varres) return end function vssize +!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ + +!\addtotable logical function inveq +!\begin{verbatim} + logical function inveq(phases,ceq) +! Only called for mapping tie-lines not in plane. If tie-lines in plane +! then all nodes are invariants. + integer phases + type(gtp_equilibrium_data), pointer :: ceq +!\end{verbatim} + integer nrel,ii,nostph,tpvar,degf,www + type(gtp_condition), pointer :: pcond,lastcond + type(gtp_state_variable), pointer :: stvr +! How to know if the ceq is invariant? Gibbs phase rule, Degrees of freedom +! f = n + z - w - p +! where n is number of components, z=2 if T and P variable, +! z=1 if T or P variable, z=0 if both T and P fixed, +! w is number of other potential conditions (MU, AC) +! p is number of stable phases. +! write(*,*)'3Y in inveq' + nrel=noel() +! sum up nubler of stable phases and check if T and P are fixed + nostph=0 +! ntups=nooftup() +! do ii=1,noofphasetuples() + do ii=1,nooftup() + if(ceq%phase_varres(phasetuple(ii)%lokvares)%phstate.gt.0) & + nostph=nostph+1 + enddo +! loop all conditions + lastcond=>ceq%lastcondition + pcond=>lastcond + tpvar=2 + www=0 +100 continue + if(pcond%active.eq.0) then +! condtion is active + stvr=>pcond%statvar(1) +! statevarid 1 is T and 2 is P + if(stvr%statevarid.eq.1 .or. stvr%statevarid.eq.2) then +! Hm, ceq is not the equilibrium record for the node point ... + tpvar=tpvar-1 + elseif(stvr%statevarid.lt.10) then +! potential/activity condition for a component + www=www+1 + endif + endif + pcond=>pcond%next + if(.not.associated(pcond,lastcond)) goto 100 +! +! Hm again, ignore tpvar? +! degf=nrel+tpvar-www-nostph + degf=nrel-www-nostph +! write(*,'(a,8i4)')'3Y in inveq 2',nrel,tpvar,www,nostph,degf + if(degf.lt.0) then +! We have an invariant equilibrium, return the number of stable phases + phases=nostph + inveq=.true. +! write(*,200)'3Y An invariant equilibrium!',nrel,tpvar,nostph,phases +!200 format(a,5i7) + else +! write(*,210)degef,nrel,tpvar,phases +!210 format('3Y not invariant eq, elements, stable phases: ',4i4) +! if not invariant isoplet node there are 3 exits (2 lines crossing) + phases=nostph + inveq=.false. + endif +1000 continue + return + end function inveq + !/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ diff --git a/models/gtp3Z.F90 b/models/gtp3Z.F90 index b94f290..06c28d5 100644 --- a/models/gtp3Z.F90 +++ b/models/gtp3Z.F90 @@ -432,7 +432,6 @@ subroutine ct1xfn(string,ip,nc,coeff,koder,fromtdb) ! ! for TDB compatibility skip # ! -! DO NOT allow unary functions ABOVE(TB) and BELOW(TB) ! check consistency implicit none integer ip,nc,koder(5,*) @@ -454,7 +453,8 @@ subroutine ct1xfn(string,ip,nc,coeff,koder,fromtdb) ! NEIN is the Einstein function ! MAX1 is 1.0 if argument is larger than 1.0, error if argument negative ! LOG is LOG10 and LN is the natural logarithm!!! - DATA unary/'LOG ','LN ','EXP ','ERF ','XNEIN ','MAX1 '/ + DATA unary/'LOG ','LN ','EXP ','ERF ','INTEIN','MAX1 '/ +! DATA unary/'LOG ','LN ','EXP ','ERF ','XNEIN ','MAX1 '/ ! DATA unary/'LOG ','LN ','EXP ','ERF ','XNEIN '/ ! DATA unary/'LOG ','LN ','EXP ','ERF ','ABOVE ','BELOW '/ ! @@ -618,7 +618,7 @@ subroutine ct1xfn(string,ip,nc,coeff,koder,fromtdb) call store_tpfun_dummy(symbol) else ! otherwise give error message -! write(*,*)'TPFUN Unknown symbol: ',symbol,freetpfun-1 + write(*,*)'TPFUN contain unknown symbol: ',symbol,freetpfun-1 gx%bmperr=4002; goto 1000 endif ! we have found the symbol @@ -995,7 +995,9 @@ end subroutine ct1mexpr subroutine ct1efn(inrot,tpval,val,tpres) !...evaluates a datastructure of an expression. Value returned in val ! inrot is root expression tpfunction record -! tpval is valuse of T and P, symval is values of symbols +! tpval is valuse of T and P, +! val is array of values calculated here +! tpres is array of all calculated functions ! first and second derivatives of T and P also calculated and returned ! in order F, F.T, F.P, F.T.T, F.T.P, F.P.P ! @@ -1092,6 +1094,7 @@ subroutine ct1efn(inrot,tpval,val,tpres) link4=0 ! if(link.ne.0) write(*,201)'funev: ',lrot,ic,link,link3 !201 format(a,10i5) +! write(*,'(a,5i4,1pe12.4)')'3Z intein2: ',ic,ipow,link,link3,link4,ff if(link.lt.0 .and. link3.gt.1000) then ! if link is negative (unary funktion) and link3 is >1000 (link) ! we must evaluate link3 first @@ -1154,6 +1157,7 @@ subroutine ct1efn(inrot,tpval,val,tpres) exprot%wpow(ic)=link4 endif !------------------------------------------------------------- +! write(*,'(a,5i4,1pe12.4)')'3Z intein3: ',ic,ipow,link,link3,link4,ff evlink: if(link.gt.0) then ! link to another symbol, extract its value and use chain rule ! extract the results from the symbol if already calculated @@ -1317,6 +1321,7 @@ subroutine ct1efn(inrot,tpval,val,tpres) elseif(link.lt.0) then !------------------------------------------------------ ! unary function, next term is argument, not very elegant .... +! write(*,'(a,5i4,1pe12.4)')'3Z intein4: ',ic,ipow,link,link3,link4,ff unfun=link cc=exprot%coeffs(ic+1) ! cc should never be zero here, if so bug in the parser @@ -1362,7 +1367,6 @@ subroutine ct1efn(inrot,tpval,val,tpres) if(abs(tpres(link2)%tpused(1)-tpval(1)).lt.& mini*tpres(link2)%tpused(1) .and. & abs(tpres(link2)%tpused(2)-tpval(2)).lt.& -! mini*tpres(link2)%tpused(2)) then mini*tpres(link2)%tpused(2) .and. & ! added this check as it seems new assessment coefficients are nor used!! tpres(link2)%forcenewcalc.eq.tpfuns(link2)%forcenewcalc) then @@ -1430,8 +1434,10 @@ subroutine ct1efn(inrot,tpval,val,tpres) endif ! now combine term1 and term2 using chain rule. link values are ! -1: LOG, -2: LN, -3: EXP, -4: ERF, only LN and EXP implemented below -! -5: NEIN, is the Einstein function +! -5: INTEIN, is the Einstein function, integrated as a Gibbs energy +! the argument is the Einstein T ! -6: MAX1, if argument <0 ERROR, if >1 replace by 1 +! write(*,'(a,5i4,1pe12.4)')'3Z intein5: ',ic,ipow,link,link3,link4,ff evunfun: if(unfun.eq.-1) then ! LOG base 10 ! ff=ff*Log10(gg) added by Sheng Yen Li @@ -1480,9 +1486,21 @@ subroutine ct1efn(inrot,tpval,val,tpres) write(*,*)'Error function not implemented' stop 71 elseif(unfun.eq.-5) then -! EINSTEIN FUNCTION - write(*,*)'Einstein function not implemented' - stop 72 +! INTEGRATED EINSTEIN: INTEIN = 1.5*R + 3*R*T*LN(EXP(THETA/T)+1), THETA=gg + if(dfdt.ne.zero) then + write(*,*)'3Z GEIN must not be multiplied with T!' + gx%bmperr=4399; goto 1000 + endif +! write(*,'(a,5i5,1pe12.4)')'3Z intein6: ',ic,ipow,link,link3,link4,ff +! ff is the constant factor in front of the Einstein function +! It is overwritten by the Einsten function (multiplied by original ff) + call tpfun_geinstein(tpval,gg,ff,dfdt,dfdp,d2fdt2,d2fdtdp,d2fdp2) +! write(*,'(a,i3,6(1pe12.4))')'3Z call Einstein:',link,& +! gg,ff,dfdt,d2fdt2 +! ff is the coefficient for the Einstein Functions, should be a constant ...? +! write(*,*)'Einstein function not implemented' +! stop 72 + if(gx%bmperr.ne.0) goto 1000 elseif(unfun.eq.-6) then ! MAX1 function, used for SRO .... function and derivatives in gg, dgdt etc. ! write(*,*)'MAX1 function',gg @@ -1499,27 +1517,8 @@ subroutine ct1efn(inrot,tpval,val,tpres) d2fdp2=zero; d2fdtdp=zero; d2fdt2=zero dfdp=zero; dfdt=zero; ff=one endif -! elseif(unfun.eq.-5) then -! above(T0): ff=ff*above(gg). gg is t0, breakfun(1..6) are fun and derivatives -! NO P-derivatives means breakfun(3)=breakfun(5)=breakfun(6)=0 -! call above_t0_calc(gg,tpval,breakfun) -! d2fdp2=breakfun(1)*d2fdp2 -! d2fdtdp=breakfun(1)*d2fdtdp+dfdp*breakfun(2) -! d2fdt2=breakfun(1)*d2fdt2+ff*breakfun(4)+2.0D0*dfdt*breakfun(2) -! dfdp=breakfun(1)*dfdp -! dfdt=breakfun(1)*dfdt+ff*breakfun(2) -! ff=ff*breakfun(1) -! elseif(unfun.eq.-6) then -! below T0, ff=ff*above(gg). gg is t0, breakfun(1..6) are fun and derivatives -! NO P-derivatives means breakfun(3)=breakfun(5)=breakfun(6)=0 -! call below_t0_calc(gg,tpval,breakfun) -! d2fdp2=breakfun(1)*d2fdp2 -! d2fdtdp=breakfun(1)*d2fdtdp+dfdp*breakfun(2) -! d2fdt2=breakfun(1)*d2fdt2+ff*breakfun(4)+2.0D0*dfdt*breakfun(2) -! dfdp=breakfun(1)*dfdp -! dfdt=breakfun(1)*dfdt+ff*breakfun(2) -! ff=ff*breakfun(1) else +! undefined function gx%bmperr=4021 goto 1000 endif evunfun @@ -1588,6 +1587,54 @@ end subroutine ct1efn !level %wpow link4 !/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ +!\addtotable subroutine tpfun_geinstein +!\begin{verbatim} + subroutine tpfun_geinstein(tpval,gg,ff,dfdt,dfdp,d2fdt2,d2fdtdp,d2fdp2) +! evaluates the integrated Einstein function (including 1.5*R) INTEIN +! gg is the value of the Einstein THETA +! ff is a constant factor which should be multiplied with all terms +! ff is overwritten with the Einstein function (multiplied with ff in) +! the other parameters are derivatives of the integrated Einstein function + implicit none + double precision tpval(*) + double precision gg,ff,dfdt,dfdp,d2fdt2,d2fdtdp,d2fdp2 +!\end{verbatim} + double precision kvot,kvotexpkvotm1,expmkvot,lnexpkvot,ww,rgas +! return ff = 1.5*R + 3*R*T*LN(EXP(-gg/T) + 1) and derivatives + rgas=globaldata%rgas +! write(*,*)'3Z in Einstein function',gg,tpval(1) + ww=ff + kvot=gg/tpval(1) + if(kvot.gt.2.0d2) then +! handle extreme values of kvot, we divide by kvotexpkvotm1**2 by expmkvot bolw + expmkvot=one + kvotexpkvotm1=zero + lnexpkvot=zero +! write(*,'(a,5(1pe12.4))')'3Z Einetsin 1: ',kvot,expmkvot,& +! kvotexpkvotm1,lnexpkvot + else + expmkvot=exp(-kvot) + kvotexpkvotm1=kvot/(exp(kvot)-one) + lnexpkvot=log(one-expmkvot) +! write(*,'(a,5(1pe12.4))')'3Z Einetsin 2: ',kvot,expmkvot,& +! kvotexpkvotm1,lnexpkvot + endif +! this is the integral G contribution from an Einstein solid + ff=1.5d0*rgas*gg*ww + 3.0D0*rgas*tpval(1)*lnexpkvot*ww + dfdt=3.0d0*rgas*(lnexpkvot-kvotexpkvotm1)*ww +! write(*,10)rgas,kvot,lnexpkvot,kvotexpkvotm1,dfdt +!10 format('3Z bug: ',6(1pe12.4)) + dfdp=zero +! this is the second derivative of G wrt T; i.e. the Einstein solid Cp equation + d2fdt2=-3.0d0*rgas*kvotexpkvotm1**2/(expmkvot*tpval(1))*ww + d2fdtdp=zero + d2fdp2=zero +1000 continue + return + end subroutine tpfun_geinstein + + !/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ + !\addtotable subroutine ct1wfn !\begin{verbatim} subroutine ct1wfn(exprot,tps,string,ip) @@ -1605,7 +1652,8 @@ subroutine ct1wfn(exprot,tps,string,ip) character ch1*1,cht*1,extsym*(lenfnsym),unary(nunary)*6 TYPE(tpfun_expression), pointer :: exprot ! these should be the same as in ct1xfn !!! ?? - DATA unary/'LOG ','LN ','EXP ','ERF ','XNEIN ','MAX1 '/ + DATA unary/'LOG ','LN ','EXP ','ERF ','INTEIN','MAX1 '/ +! DATA unary/'LOG ','LN ','EXP ','ERF ','GEIN ','MAX1 '/ ! if(.not.associated(exprot)) then string(ip:ip+2)='0; ' @@ -1739,6 +1787,7 @@ subroutine ct1wfn(exprot,tps,string,ip) endif symbol ! no symbol or unary function, coefficient with possible powers if(coeff(level).ne.one) then +! if 4th argument >0 then write a sign call wrinum(string,ip,12,6,coeff(level)) mult=-1 else @@ -1750,7 +1799,10 @@ subroutine ct1wfn(exprot,tps,string,ip) if(koder(i,level).ne.0) goto 219 enddo ! without this the Inden magnetic function will miss its initial 1.0 - call wrinum(string,ip,2,0,coeff(level)) +! call wrinum(string,ip,2,0,coeff(level)) +! changed 20.03.17/BoS because EXP(T)+1 missed the + between ) and 1 +! Force wrinum to write positive signs by 4th parameter positive + call wrinum(string,ip,2,1,coeff(level)) goto 220 219 continue ! missing coefficient discovered by Mauro, as the coefficient is unity diff --git a/pmain1.F90 b/pmain1.F90 index 3e3ff73..088c623 100644 --- a/pmain1.F90 +++ b/pmain1.F90 @@ -23,7 +23,7 @@ PROGRAM pmain1 linkdate=date(1:4)//'-'//date(5:6)//'-'//date(7:8) ! for example: linkdate='2019-11-27' ! this is the overall version identifier - version=' 6.000 ' + version=' 6.001 ' ! intvar and dblvar will eventually be used for allocations intvar(1)=30 call init_gtp(intvar,dblvar) diff --git a/stepmapplot/smp2.F90 b/stepmapplot/smp2.F90 index ad10d5d..2cde97c 100644 --- a/stepmapplot/smp2.F90 +++ b/stepmapplot/smp2.F90 @@ -65,7 +65,7 @@ MODULE ocsmp ! with same composition (see for example U-O) integer, parameter :: EXCLUDEDLINE=0, TWOSTOICH=1 ! Bit in the MAP_NODE record - integer, parameter :: MAPINVARIANT=0 + integer, parameter :: MAPINVARIANT=0,STEPINVARIANT=1 !\end{verbatim} ! !\begin{verbatim} @@ -90,12 +90,12 @@ MODULE ocsmp ! if we have 3 or more axis there can be 2 or more fix phases along the line?? ! With 2 axis there can only be one fix phase! type(gtp_phasetuple), dimension(:), allocatable :: linefixph -! also save index to phr!! +! also save index to phr!! (do not trust ...) integer, dimension(:), allocatable :: linefix_phr ! Save the phase tuplet representing the phase fix at start node here ! If it wants to be stable at first step along a line change axis direction -! type(gtp_phasetuple) :: nodfixph -! This is the phase index in the phr array (ncludes both index and compset) +! type(gtp_phasetuple) :: nodfixph <<<< not used +! This is the phase index in the phr array (phr has both phase and compset) integer nodfixph ! We must also save the number and set of stable phases and theit amounts ! as we will have different stable phases for different lines diff --git a/stepmapplot/smp2A.F90 b/stepmapplot/smp2A.F90 index 93c8ab1..b09fe74 100644 --- a/stepmapplot/smp2A.F90 +++ b/stepmapplot/smp2A.F90 @@ -273,8 +273,9 @@ subroutine map_doallines(maptop,nax,axarr,seqxyz,starteq) ! write(*,*)'problems',mapline%problems,ceq%tpval(1) ! endif ! write(*,*)'Calling calceq7 with T=',ceq%tpval(1),mapline%axandir +! write(*,*)'Calling calceq7 with meqrec%status:',meqrec%status call calceq7(mode,meqrec,mapfix,ceq) -! write(*,*)'Back from calceq7A ',gx%bmperr +! write(*,*)'SMP2A Back from calceq7 ',gx%bmperr,meqrec%status if(gx%bmperr.ne.0) then ! error 4187 is to set T or P to less than 0.1 if(gx%bmperr.eq.4187) then @@ -947,6 +948,7 @@ subroutine map_startpoint(maptop,nax,axarr,seqxyz,inactive,ceq) ! (without global minimization). We will save the meq_setup record! ! write(*,*)'meq_startpoint: allocating meqrec' allocate(meqrec) + meqrec%status=0 ! We must use mode=-1 for map_replaceaxis below has to calculate several equil ! and the phr array must not be deallocated. mapfix will be used later to ! indicate fix and stable phases for different lines (maybe ...) @@ -1102,6 +1104,10 @@ subroutine map_startpoint(maptop,nax,axarr,seqxyz,inactive,ceq) !----------------------------------------------------------------------- if(ocv()) write(*,*)'allocating lineheads: ',ieq,maptop%seqy allocate(mapnode%linehead(ieq)) +! meqrec%status + do jp=1,ieq + mapnode%linehead(ieq)%meqrec%status=0 + enddo ! we can have 3 or more exits if starting inside a 3 phase triagle for isotherm if(ieq.eq.2) then ! STEP command: set one exit in each direction of the active axis axactive @@ -1415,7 +1421,7 @@ subroutine map_replaceaxis(meqrec,axactive,ieq,nax,axarr,tmpline,& fixphaseloop: do jph=1,3 ! all phases are already set as stable kj=meqrec%stphl(jph) - write(*,*)'tmpline 1: ',jph,kj +! write(*,*)'tmpline 1: ',jph,kj if(jph.gt.1) then allocate(tmpline(jph)%linefixph(1)) allocate(tmpline(jph)%linefix_phr(1)) @@ -1429,15 +1435,15 @@ subroutine map_replaceaxis(meqrec,axactive,ieq,nax,axarr,tmpline,& ! meqrec%fixpham(1)=zero tmpline(jph)%nfixphases=1 tmpline(jph)%linefixph(1)=zerotup - write(*,*)'tmpline 2A: ',jph,kj - write(*,*)'tmpline 2C: ',allocated(meqrec%phr) - write(*,*)'tmpline 2B: ',meqrec%phr(kj)%iph - write(*,*)'tmpline 2C: ',allocated(tmpline(jph)%linefixph) +! write(*,*)'tmpline 2A: ',jph,kj +! write(*,*)'tmpline 2C: ',allocated(meqrec%phr) +! write(*,*)'tmpline 2B: ',meqrec%phr(kj)%iph +! write(*,*)'tmpline 2C: ',allocated(tmpline(jph)%linefixph) tmpline(jph)%linefix_phr(1)=kj tmpline(jph)%linefixph(1)%ixphase=meqrec%phr(kj)%iph - write(*,*)'SMP2A tmpline 3: ',jph,kj +! write(*,*)'SMP2A tmpline 3: ',jph,kj tmpline(jph)%linefixph(1)%compset=meqrec%phr(kj)%ics - write(*,*)'SMP2A tmpline 4: ',jph,kj +! write(*,*)'SMP2A tmpline 4: ',jph,kj tmpline(jph)%nstabph=1 kph=jph+1 if(kph.gt.3) kph=1 @@ -2218,7 +2224,7 @@ subroutine map_changeaxis(mapline,nyax,oldax,nax,axarr,axval,bytax,ceq) mapline%meqrec%tpindep(1)=.FALSE. if(ocv()) write(*,*)'Marking that T is variable' elseif(pcond%statev.eq.2) then - mapline%meqrec%tpindep(1)=.FALSE. + mapline%meqrec%tpindep(2)=.FALSE. endif !------------------------------------------- ! this is the old axis with active condition, look for its condition @@ -2239,7 +2245,7 @@ subroutine map_changeaxis(mapline,nyax,oldax,nax,axarr,axval,bytax,ceq) if(ocv()) write(*,*)'Marking that T is variable',ceq%tpval(1) ceq%tpval(1)=pcond%prescribed elseif(pcond%statev.eq.2) then - mapline%meqrec%tpindep(1)=.TRUE. + mapline%meqrec%tpindep(2)=.TRUE. endif !---------------------------------------------------------- ! now we calculate the same equilibrium but with different axis condition! @@ -3580,7 +3586,7 @@ subroutine map_bytfixphase(mapline,axis,meqrec,xxx,ceq) meqrec%fixpham(1)=zero meqrec%iphl(1)=mapline%linefixph(1)%ixphase meqrec%icsl(1)=mapline%linefixph(1)%compset -!------------- now the stable phase +!------------- now the stable phase ?? value of stable_phr=?? mapline%stableph(1)=phtup1 mapline%stable_phr(1)=phrix ! write(*,*)'SMP2A phrix switching fix/stable phase: ',phrix @@ -3694,7 +3700,7 @@ subroutine map_calcnode(irem,iadd,maptop,mapline,meqrec,axarr,ceq) if(ocv()) write(*,*)'Marking that T is variable' maxstph=1 elseif(pcond%statev.eq.2) then - meqrec%tpindep(1)=.TRUE. + meqrec%tpindep(2)=.TRUE. maxstph=1 endif !-------------------------------------------- @@ -3872,7 +3878,7 @@ subroutine map_calcnode(irem,iadd,maptop,mapline,meqrec,axarr,ceq) axval,ceq%tpval(1) 55 format(a,6(1pe14.6)) elseif(pcond%statev.eq.2) then - meqrec%tpindep(1)=.FALSE. + meqrec%tpindep(2)=.FALSE. ! ceq%tpval(2)=value endif ! write(*,*)'error in map_calcnode, remove phfix: ',phfix @@ -3908,7 +3914,7 @@ subroutine map_calcnode(irem,iadd,maptop,mapline,meqrec,axarr,ceq) ! When we are there we have successfully calculated an equilibrium with a ! new phase set create a node with this equilibrium and a new line records 500 continue - write(*,*)'SMP2 Successful calculation of a node point',phfix +! write(*,*)'SMP2 Successful calculation of a node point',phfix ! phfix is set negative if phase should be removed ! NOTE the phase set fix in the node may not be the same which ! wanted to disappear/appear when calling the map_calcnode!! @@ -3948,7 +3954,7 @@ subroutine map_calcnode(irem,iadd,maptop,mapline,meqrec,axarr,ceq) ceq%tpval(1)=value if(ocv()) write(*,*)'Marking that T is a condition again',value elseif(pcond%statev.eq.2) then - meqrec%tpindep(1)=.FALSE. + meqrec%tpindep(2)=.FALSE. ceq%tpval(2)=value endif ! Save this as the last equilibrium of the line @@ -4028,7 +4034,7 @@ subroutine map_calcnode(irem,iadd,maptop,mapline,meqrec,axarr,ceq) ! the number of lines ending at an invariant isopleth is 2*haha ! current number of stable phases is meqrec%nstph. ! sign(1,phfix) is 1 if phfix>0; -1 if phfix<0 - write(*,21)meqrec%nstph,haha,phfix,meqrec%nstph-haha+sign(1,phfix) +! write(*,21)meqrec%nstph,haha,phfix,meqrec%nstph-haha+sign(1,phfix) 21 format('SMP2A stable phases mm: ',3i5,i10) endif endif @@ -4091,9 +4097,11 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& !\end{verbatim} type(gtp_equilibrium_data), pointer :: newceq,tmpceq type(map_node), pointer :: mapnode,newnode - type(map_line), pointer :: linenode + type(map_line), pointer :: linenode,tmpline type(gtp_condition), pointer :: pcond type(gtp_state_variable), pointer :: svrrec,svr2 + type(map_fixph), allocatable :: mapfix + type(meq_setup), pointer :: meqrec2 type(gtp_state_variable), target :: svrtarget integer remph,addph,nel,iph,ics,jj,seqx,nrel,jphr,stabph,kph,kcs,kk,lfix integer zph,stepax,kpos,seqy,jp,nopotax,lokcs,lokph,haha2,linefphr @@ -4103,14 +4111,18 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& character eqname*24,phases*60 double precision stepaxval,middle,testv,xxx ! mark that line ended with two stoichometric phases and data for isopleth inv - integer twostoichset,errall,jfix,jstab,jlast,kstab,zp,phrix - integer onlyone,notone,jused,zz,tz,qy + integer twostoichset,errall,jfix,jstab,jlast,kstab,zp,phrix,infix3 + integer onlyone,notone,jused,zz,tz,qy,savefix1,savefix2,nodein(2),nodeut(2) +! lifexix, nodefix and prevfix are used to fix pair of phases that have zero +! amount at the exit points of lines from an invariant equilibrium. + integer linefix,nodefix,infix,infix2,doubline,twice,firstoutfix,outfix integer, allocatable, dimension(:,:) :: invph + double precision, allocatable, dimension(:) :: exitcomp,eqcopy + logical stepinvariantnode ! lfix=0 ! the phase kept fix with zero amount at the node is phfix It can be ! negative at STEP if it is a phase that will dissapear. -! write(*,*)'In map_newnode with phase change',phfix,mapnode%seqx,ceq%tpval(1) if(ocv()) write(*,87)'We are in map_newnode with a fix phase: ',& phfix,ceq%tpval(1) 87 format(a,i4,2x,1pe12.4) @@ -4125,7 +4137,9 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& mapnode=>maptop nrel=meqrec1%nrel 100 continue +!--------------------------------------------------------------------- ! loop all mapnodes to check if any has the same chemical potentials +!--------------------------------------------------------------------- ! write(*,*)'Comparing with node: ',mapnode%seqx,nrel ! write(*,105)'T diff: ',ceq%tpval(1),mapnode%tpval(1),& ! abs(ceq%tpval(1)-mapnode%tpval(1)),abs(vz*mapnode%tpval(1)) @@ -4143,9 +4157,9 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& 105 format(a,5(1pe16.8)) if(abs(ceq%cmuval(nel)-mapnode%chempots(nel)).gt.& abs(2.0D1*vz*mapnode%chempots(nel))) then - write(*,'(a,3(1pe12.4))')'Not same chempots, compare next node',& - abs(ceq%cmuval(nel)-mapnode%chempots(nel)),& - abs(2.0D1*vz*mapnode%chempots(nel)) +! write(*,'(a,3(1pe12.4))')'SMP not same chempots, at this node',& +! abs(ceq%cmuval(nel)-mapnode%chempots(nel)),& +! abs(2.0D1*vz*mapnode%chempots(nel)) goto 110 endif enddo @@ -4262,10 +4276,8 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& ! save index of the phase set fix at the node ! write(*,*)'SMP Saving index of new fix phase: ',abs(phfix) if(phfix.lt.0) then -! newnode%nodefix%phaseix=-meqrec1%phr(abs(phfix))%iph newnode%nodefix%ixphase=-meqrec1%phr(abs(phfix))%iph else -! newnode%nodefix%phaseix=meqrec1%phr(abs(phfix))%iph newnode%nodefix%ixphase=meqrec1%phr(abs(phfix))%iph endif newnode%nodefix%compset=meqrec1%phr(abs(phfix))%ics @@ -4326,11 +4338,10 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& ! an invariant equlibrium with haha stable phases and 2*haha-1 exiting lines ! newnode%lines=2*jj-1 if(inveq(haha2,ceq)) then - newnode%lines=2*haha-1 - write(*,*)'SMP2A *** invariant node with exits: ',newnode%lines -! write(*,*)'SMP2A Currently this is ignored' -! uncommenting next line means isopleth inveriants are ignored -! newnode%lines=3 +! newnode%lines=2*haha-1 +! Only few of the exit lines will be in the plane f the diagram. Assume 8 +! i.e. there will be 7 exits + newnode%lines=7 else newnode%lines=3 endif @@ -4395,10 +4406,12 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& write(*,*)'SMP2A Allocation error 1: ',errall gx%bmperr=4370; goto 1000 endif +! do jp=1,newnode%lines !--------------------- code moved from map_findline ! COPY of the equilibrium record from newnode to newnode%linehead(jp)%lineceq if(ocv()) write(*,*)'We found a line from node: ',mapnode%seqx + newnode%linehead(jp)%meqrec%status=0 eqname='_MAPLINE_' kpos=10 seqy=maptop%seqy+1 @@ -4424,7 +4437,7 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& ! MAP tie-line in plane 2; isopleth non-invariant 3; isopleth invariant >3 kpos=min(newnode%lines,4) ! select case(newnode%lines) - select case(kpos) + exits: select case(kpos) !========================================================================== case default write(*,*)'SMP2A node error: exit lines= ',newnode%lines @@ -4432,6 +4445,7 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& !========================================================================== case(1)! step node with just one exit ! If phfix negative the fix phase wants to dissapear + if(inveq(jj,ceq)) write(*,'(a)')'This is an invariant node' changephaseset: if(phfix.lt.0) then ! remove a phase --------------------------- remph=-phfix @@ -4466,15 +4480,90 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& ! we have to add phase phfix to the stable phase set, that is no problem ! as it is already in all lists, just remove that it should be fix ! write(kou,88)phfix,' appears, ',meqrec1%nstph -! meqrec1%nfixph seems not to be used .... +! meqrec1%nfixph seems not to be used .... ?? if(meqrec1%nfixph.gt.0) then meqrec1%fixph(1,meqrec1%nfixph)=0 meqrec1%fixph(2,meqrec1%nfixph)=0 meqrec1%phr(phfix)%phasestatus=PHENTSTAB meqrec1%nfixph=meqrec1%nfixph-1 endif + stepinvariantnode=.FALSE. + if(inveq(jj,ceq)) then +! +! if node is invariant we must remove one phase, which? NOT phfix + stepinvariantnode=.TRUE. + newnode%status=ibset(newnode%status,STEPINVARIANT) +! write(*,*)'SMP2A invariant node at step',phfix,newnode%status,& +! meqrec1%nstph +! set the invariant bit in the node and calculate en equilibrium at +! a very small step above the invariant to find the new set of phases + newnode%linehead(1)%meqrec=meqrec1 + tmpline=>newnode%linehead(1) +! do kk=1,meqrec1%nstph +! jj=meqrec1%stphl(kk) +! write(*,294)'SMP initial set of phases: ',kk,jj,& +! meqrec1%phr(jj)%iph,meqrec1%phr(jj)%curd%amfu +! enddo + meqrec2=>tmpline%meqrec +! do kk=1,meqrec2%nstph +! jj=meqrec2%stphl(kk) +! write(*,294)'SMP same initial set of phases: ',kk,jj,& +! meqrec2%phr(jj)%iph,meqrec2%phr(jj)%ics,& +! meqrec2%phr(jj)%curd%amfu +! enddo +294 format(a,3i5,i2,1pe14.6) + call locate_condition(axarr(1)%seqz,pcond,tmpline%lineceq) + if(gx%bmperr.ne.0) goto 100 +! call list_conditions(kou,tmpline%lineceq) +! call list_sorted_phases(kou,tmpline%lineceq) +! if(gx%bmperr.ne.0) goto 100 + pcond%prescribed=pcond%prescribed+& + 1.0D-3*stepax*axarr(1)%axinc +! call list_conditions(kou,tmpline%lineceq) +! write(*,*)'SMP small step invariant to find phase which disappear' + call calceq3(0,.FALSE.,tmpline%lineceq) +! first argument -1 to keep the datastructure in meqrec2 +! call calceq7(-1,meqrec2,mapfix,tmpline%lineceq) +! write(*,*)'Back from calceqx',gx%bmperr,meqrec1%nstph +! call list_sorted_phases(kou,tmpline%lineceq) +! NOTE the content of meqrec2 has not been updated as calceq3 creates a new +! independent meqrec structure. We must copy the values of phase amounts +! using a pointer directly to the phase_varres record +! list amount of phases after this small step. However, the layout of +! the meqrec records are the same, we can use phase indices and other things + do kk=1,meqrec2%nstph-1 + jj=meqrec2%stphl(kk) + call get_phase_compset(meqrec2%phr(jj)%iph,meqrec2%phr(jj)%ics,& + lokph,lokcs) + xxx=tmpline%lineceq%phase_varres(lokcs)%amfu +! write(*,294)'SMP new set of phases: ',kk,jj,& +! meqrec2%phr(jj)%iph,meqrec2%phr(jj)%ics,xxx + if(xxx.gt.zero) then + meqrec2%phr(jj)%curd%amfu=xxx + else + do zz=kk,meqrec2%nstph-1 + meqrec2%stphl(zz)=meqrec2%stphl(zz+1) + enddo + endif + enddo + meqrec2%nstph=meqrec2%nstph-1 +! do kk=1,meqrec2%nstph +! jj=meqrec2%stphl(kk) +! write(*,294)'SMP final initial set of phases: ',kk,jj,& +! meqrec2%phr(jj)%iph,meqrec2%phr(jj)%ics,& +! meqrec2%phr(jj)%curd%amfu +! enddo +! finally copy to meqrec1 ... + meqrec1%nstph=meqrec2%nstph + do kk=1,meqrec2%nstph + meqrec1%stphl(kk)=meqrec2%stphl(kk) +! I assume the amounts are not needed, they should already be in lineceq ...?? + enddo +! rearrange the array of stable phases, one should be removed +! stop 'all OK?' + endif else - write(*,*)'This is another never never error' + write(*,*)'This is another never never error',phfix gx%bmperr=4239; goto 1000 endif changephaseset ! set values in linhead record @@ -4492,16 +4581,29 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& newnode%linehead(1)%axfact=1.0D-2 newnode%linehead(1)%nfixphases=0 ! try to get a nice output of stable phases below - allocate(newnode%linehead(1)%stableph(meqrec1%nstph)) - allocate(newnode%linehead(1)%stable_phr(meqrec1%nstph)) - newnode%linehead(1)%nstabph=0 - do iph=1,meqrec1%nstph - newnode%linehead(1)%nstabph=newnode%linehead(1)%nstabph+1 - jj=meqrec1%stphl(iph) - newnode%linehead(1)%stableph(iph)%ixphase=meqrec1%phr(jj)%iph - newnode%linehead(1)%stableph(iph)%compset=meqrec1%phr(jj)%ics - newnode%linehead(1)%stable_phr(iph)=jj - enddo +! if(stepinvariantnode) then +! allocate(newnode%linehead(1)%stableph(meqrec2%nstph)) +! allocate(newnode%linehead(1)%stable_phr(meqrec2%nstph)) +! newnode%linehead(1)%nstabph=0 +! do iph=1,meqrec2%nstph +! newnode%linehead(1)%nstabph=newnode%linehead(1)%nstabph+1 +! jj=meqrec2%stphl(iph) +! newnode%linehead(1)%stableph(iph)%ixphase=meqrec2%phr(jj)%iph +! newnode%linehead(1)%stableph(iph)%compset=meqrec2%phr(jj)%ics +! newnode%linehead(1)%stable_phr(iph)=jj +! enddo +! else + allocate(newnode%linehead(1)%stableph(meqrec1%nstph)) + allocate(newnode%linehead(1)%stable_phr(meqrec1%nstph)) + newnode%linehead(1)%nstabph=0 + do iph=1,meqrec1%nstph + newnode%linehead(1)%nstabph=newnode%linehead(1)%nstabph+1 + jj=meqrec1%stphl(iph) + newnode%linehead(1)%stableph(iph)%ixphase=meqrec1%phr(jj)%iph + newnode%linehead(1)%stableph(iph)%compset=meqrec1%phr(jj)%ics + newnode%linehead(1)%stable_phr(iph)=jj + enddo +! endif ! end attempt newnode%linehead(1)%firstinc=1.0D-2*axarr(1)%axinc*mapline%axandir ! newnode%linehead(1)%firstinc=1.0D-3*axarr(1)%axinc*mapline%axandir @@ -4645,6 +4747,7 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& ! the previously fix phase is set as entered with stablepham as initial amount newnode%linehead(1)%stableph(1)%ixphase=iph newnode%linehead(1)%stableph(1)%compset=ics +! value of %stable_phr=?? newnode%linehead(1)%stable_phr(1)=phrix newnode%linehead(1)%stablepham(1)=one ! store the phase number that must not become stable in nodfixph @@ -5015,14 +5118,14 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& ! write(*,*)'Created node with 2 exits: ',newnode%seqx,ceq%tpval(1) 390 continue !--------------------------------------------------------------- -! invariant isopleth, more then 3 exits +! invariant isopleth, more than 3 exits case(4) ! isopleth invariants for isopleths, inveq ! number of stable phases equal to components+1 ! number of adjacent regions with "components" stable phases is "components+1" -! number of exit lines are 2*(components+1) -! each line has a fix phase and one of the phases stable at the invariant -! as not stable. The remaining phases are entered. -! Each phase is fix for two lines and not stable for two others +! number of exit lines are 2*(components+1) ?? limit to 8 (minus 1 for entry) +! each line has a fix phase and one of the phases is stable at the invariant +! (set as not "nodefix"). The remaining phases are entered. +! Each phase is fix for two lines and "nodefix" for two others ! This is the way to generate the exit lines: ! - loop for all phases to set a phase fix (for two lines) ! - loop for the next two phases to set one phase not stable @@ -5035,7 +5138,6 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& ! 0 if both T and P fixed, p is number of stable phases. ! write(*,*)'SMP2A Generating exits from isopleth invariant',newnode%lines ! Two crossing lines, one in and 3 exits -! THERE IS NO CASE WHEN FINDING AN INVARIANT ! this is probably redundant, fixph already reset ! phfix is the new stable phase! Must be positive ! mapline is the just finished line @@ -5059,7 +5161,6 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& exit flfix2 endif enddo flfix2 -! write(*,*)'SMP2A found lfix?',lfix if(lfix.eq.0) stop 'ERROR' ! this is total number of phases at each the invariant ! 1 fix and stabph-2 should be stable at each exit @@ -5070,62 +5171,89 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& endif ! Collect all stable phases to be used as different exits. ! invph(1,jj) is iph,, invph(2,jj) is ics; invph(3,jj) is index in meqrec1%phr -! invph(4,jj) is to count number of times this phase is forbidden, max 2 - allocate(invph(4,stabph+2)) +! invph(4,jj) is to count number of times jj has been linefix +! invph(5,jj) is to count number of times jj has been nodefix +! invph(6,jj) is index to phase_varres + allocate(invph(6,stabph+2)) invph=0 do jj=1,stabph ! stableph is a phase_tuple invph(1,jj)=mapline%stableph(jj)%ixphase invph(2,jj)=mapline%stableph(jj)%compset -! NOW ADDED stable_phr to find the index in phr !! invph(3,jj)=mapline%stable_phr(jj) +! stable_phr is used to find the index in phr and index to phase_varres +! I DO NOT TRUST THE VALUE, "stable_phr" +! I SHOULD REORGANIZE PHE TO BE IN PHASE TUPLE ORDER. +! THERE ARE ALWAYS PROBLEM IS IF NEW COMPOSIION SETS ARE CREATED DURING MAPPING + do zz=1,meqrec1%nphase + if(meqrec1%phr(zz)%iph.eq.invph(1,jj) .and.& + meqrec1%phr(zz)%ics.eq.invph(2,jj)) then + invph(3,jj)=zz +! if(zz.ne.mapline%stable_phr(jj)) & +! write(*,*)'SMP correction: ',jj,zz,mapline%stable_phr(jj) + endif + enddo + call get_phase_compset(invph(1,jj),invph(2,jj),lokph,lokcs) + if(gx%bmperr.ne.0) goto 1000 + invph(6,jj)=lokcs enddo ! at the end of loop jj=stabph+1; store phfix, the new fix phase, invph(1,jj)=meqrec1%phr(phfix)%iph invph(2,jj)=meqrec1%phr(phfix)%ics invph(3,jj)=phfix + call get_phase_compset(invph(1,jj),invph(2,jj),lokph,lokcs) + if(gx%bmperr.ne.0) goto 1000 + invph(6,jj)=lokcs ! this is the phase fix at incomming line, only one exit line with this fix invph(1,jj+1)=iph invph(2,jj+1)=ics invph(3,jj+1)=lfix + call get_phase_compset(invph(1,jj+1),invph(2,jj+1),lokph,lokcs) + if(gx%bmperr.ne.0) goto 1000 + invph(6,jj+1)=lokcs jlast=stabph+2 -! do jj=1,jlast -! phases=' ' -! call get_phase_name(invph(1,jj),invph(2,jj),phases) -! write(*,420)'SMP2A invrnt phase: ',jj,invph(1,jj),invph(2,jj),& -! invph(3,jj),trim(phases) -! enddo ! STABLE PHASES HAS TO BE IN PHR ORDER!! SORT invph - call sort_invph(jlast,invph) +! do kk=1,stabph+2 +! write(*,'(a,i3,2x,i3,i2,4i4)')'SMP invph:',kk,(invph(zz,kk),zz=1,6) +! enddo +! second argument is first dimenstion of invph!! + call sort_invph(jlast,6,invph) if(gx%bmperr.ne.0) goto 1000 +! write(*,*)'SMP sorted invph: ' +! do kk=1,stabph+2 +! write(*,'(a,i3,2x,i3,i2,4i4)')'SMP invph:',kk,(invph(zz,kk),zz=1,6) +! enddo do jj=1,jlast phases=' ' call get_phase_name(invph(1,jj),invph(2,jj),phases) ! keep track of the phase found at the invariant and the linefix phase - if(invph(3,jj).eq.lfix) onlyone=jj - if(invph(3,jj).eq.phfix) notone=jj -! write(*,420)'SMP2A sorted phase: ',jj,invph(1,jj),invph(2,jj),& -! invph(3,jj),trim(phases) +! nodein(1) is linefix + if(invph(3,jj).eq.lfix) then + nodein(1)=jj + linefix=jj + invph(4,nodein(1))=1 +! nodein(2) is nodefix + elseif(invph(3,jj).eq.phfix) then + nodein(2)=jj + nodefix=jj + invph(5,nodein(2))=1 + endif enddo -420 format(a,4i5,2x,a) -! the entering line had notone "forbidden", mark it is used - invph(4,notone)=1 +! the entering line found at node, nodefix, mark it is used +! the entering line had this phase fix with zero amount ! write(*,'(a,5i4)')'SMP2: linefix and nodefix: ',& ! onlyone,lfix,notone,phfix -! NOTE lfix should only be used once as fixed phase ! all the others should be fixed one two exits !-------------- ! We have to generate newnode%lines exits!! jphr=mapline%nstabph jfix=1 - qy=0 ! set bit in mapnode! if(newnode%status.ne.0) write(*,*)'SMP2 nodestatus: ',newnode%status newnode%status=ibset(newnode%status,MAPINVARIANT) - invexit: do jj=1,newnode%lines -! initiate data in map_line record -! write(*,410)'SMP2A exit: ',jj,jphr -410 format(a,10i4) +! write(*,*)'SMP2 number of exit lines: ',newnode%lines,jphr + allexit: do jj=1,newnode%lines +! initiate common data in map_line record in all exit lines newnode%linehead(jj)%number_of_equilibria=0 newnode%linehead(jj)%first=0 newnode%linehead(jj)%last=0 @@ -5154,119 +5282,203 @@ subroutine map_newnode(mapline,meqrec1,maptop,axval,lastax,axarr,& nullify(newnode%linehead(jj)%end) ! number of stable phases along all lines. Additionally a fix and a forbidden newnode%linehead(jj)%nstabph=stabph +! possible problem with meqrec%status + if(newnode%linehead(jj)%meqrec%status.ne.0) then + write(*,*)'SMP zero meqrec%status for newnode%linehead',& + newnode%linehead(jj)%meqrec%status + newnode%linehead(jj)%meqrec%status=0 + endif + enddo allexit ! -! ------------------------------------------------------------------ -! something like this has to be added to define the stable and fix phase -! different for each line. But we have to specify one fix and one forbidden - phases=' ' - first: if(jj.eq.1) then -! For jj=1 we create the only exit with LFIX (in onlyone) as fix and -! in this case another phase than PHFIX (in notone) should be forbiddem - newnode%linehead(jj)%linefixph%ixphase=invph(1,onlyone) - newnode%linehead(jj)%linefixph%compset=invph(2,onlyone) - newnode%linehead(jj)%linefix_phr=invph(3,onlyone) - call get_phase_name(invph(1,onlyone),invph(2,onlyone),phases) - zp=len_trim(phases)+2 -! write(*,*)'smp2: fix: ',trim(phases),onlyone,invph(3,onlyone) -! set a phase forbidden to become stable when line starts -! must not be the fix phase and not the phase found at the node (phfix) -! THE CODE HERE IS NOT CORRECT - klop1: do kk=jlast,1,-1 - if(kk.ne.notone .and. kk.ne.onlyone) then - newnode%linehead(jj)%nodfixph=invph(3,kk); jused=kk; - invph(4,kk)=1 - exit klop1 - endif - enddo klop1 - call get_phase_name(invph(1,jused),invph(2,jused),phases(zp:)) - zp=len_trim(phases)+2 -! write(*,*)'smp2: not: ',trim(phases),jused,invph(3,jused),kk -! NOW add the stable phases excluding onlyone and jused - zz=1 - do kk=1,stabph - if(zz.eq.onlyone) zz=zz+1 - if(zz.eq.jused) zz=zz+1 - if(zz.eq.onlyone) zz=zz+1 - newnode%linehead(jj)%stableph(kk)%ixphase=invph(1,zz) - newnode%linehead(jj)%stableph(kk)%compset=invph(2,zz) - newnode%linehead(jj)%stablepham(kk)=1.0D-2 - newnode%linehead(jj)%stable_phr(kk)=invph(3,zz) - call get_phase_name(invph(1,zz),invph(2,zz),phases(zp:)) - zp=len_trim(phases)+2 - zz=zz+1 - enddo - write(*,430)jj,phases -430 format('SMP2A invexit: ',i3,2x,a) - else -! UNIFINISHED else is for all other exit lines, jfix is updated globally -! The same phase should be twice as LINEFIX with a different phase as NODEFIXPH -! the use of jfix for selecting LINEFIX is OK but NODEFIXPH is not correct -! The pair has to be determined better +! ------------------------ ISOPLETHAL INVARIANTS EXITS ----------- +! we have to find sets of LINEFIX and NODEFIX phases for each line +! We know the LINEFIX and NODEFIX phases for the line INTO THE INVARIAT +! For the first exit line we just swich these as they are at the same point +! with zero amount of both phases +! +! (C is ?) Cfix Afix Dfix Bfix (B is ?) +! \ ABCE.. \BCE./ BCDE.. / +! \ \ / / +! \____________\/_________/ +! ABDE.. !_________ ABCDE..______! CDE... +! / /\ \ +! / ABDE.. / \ ACDE.. \ +! / /ADE.\ \ +! Dfix Bfix Cfix Afix +! +! Assume we have A as linefix and find the invariant when D becomes stable +! (the top middle point). Then there is an exit with D as linefix and A as +! nodefix. There is certainly one more exit line with A and D as linefix +! but we do not know which nodefix phases they have. And there is a final +! exit with two line where these two nodefix are linefix/nodefix. +! The subroutine FIND_INV loop each phase (not A or D) set as fix=0 and fix T +! to find an equilibrium with either A or D stable with zero amount. The +! phase fix represet B or C. This has to be done twice to geerate 4 lines. +! The final two lines have B or C as linefix/nodefix. +! To simplify these calculations all uninvolved phases are suspended +! ! -! probably one has to calculate combination of fix phases to determine -! combinations of LINEFIX and NODEFIXPHASE that is correct. -! If the invariant has n+1 phases stable (with n elements) there are -! n+1 adjacent regions with n stable phases and in between these there are -! n+1 regions with n-1 stable phases. -! But there mixing are (n+1)*n regions. For n=3 that is 3*2=6 two-phase -! regions but only 4 of these are connected to the invariant. -! In order to determine which one has to calculate equilibria with different -! set of fixed phases at the invariant. - tz=jfix - if(tz.eq.onlyone) then - jfix=jfix+1; tz=tz+1 + jj=1 + phases=' ' +! now code to create correct combination of linefix and nodefix +! first a line with nodefix as linefix and vice versa +! and old linefix as nodefix phase (sorry very confusing for me too) + newnode%linehead(jj)%linefixph%ixphase=invph(1,nodefix) + newnode%linehead(jj)%linefixph%compset=invph(2,nodefix) + newnode%linehead(jj)%linefix_phr=invph(3,nodefix) +! this is needed just to understand what is happening + call get_phase_name(invph(1,nodefix),invph(2,nodefix),phases) + zp=len_trim(phases)+3 + phases(zp-1:zp-1)='(' +! write(*,*)'smp2: fix: ',trim(phases),onlyone,invph(3,onlyone) +! set linefix=LFIX as phase forbidden to become stable when line starts +! enclose the nodefix phase with ( .... ) + newnode%linehead(jj)%nodfixph=invph(3,linefix) + call get_phase_name(invph(1,linefix),invph(2,linefix),phases(zp:)) + zp=len_trim(phases)+3 + phases(zp-2:zp-2)=')' +! at this line we have switched nodefix/linefix; nodefix is fix along the line +! the incomming line had the same fix phases so set both 4 and 5 to 1 + invph(4,nodefix)=1; invph(5,nodefix)=1 + invph(4,linefix)=1; invph(5,linefix)=1 +! NOW add the stable phases excluding linefix and nodefix + kk=0 + names: do zz=1,stabph+2 + if(zz.eq.linefix .or. zz.eq.nodefix) cycle names + kk=kk+1 + newnode%linehead(jj)%stableph(kk)%ixphase=invph(1,zz) + newnode%linehead(jj)%stableph(kk)%compset=invph(2,zz) + newnode%linehead(jj)%stablepham(kk)=1.0D-2 + newnode%linehead(jj)%stable_phr(kk)=invph(3,zz) + call get_phase_name(invph(1,zz),invph(2,zz),phases(zp:)) + zp=len_trim(phases)+2 + enddo names +! note again: nodefix here is fix along the line, linefix is stable at invarant + write(*,430)jj,nodefix,linefix,trim(phases) +430 format('SMP2A invexit: ',3i3,2x,a) +! The code above is for the FIRST exit line +! We will use nodein(1) and nodein(2) to find 4 more lines +! save the constitution of the node equilibrium, copy it to eqcopy + if(allocated(exitcomp)) deallocate(exitcomp) + call save_constitutions(newnode%nodeceq,exitcomp) +! this is used to save the stable phases and constituion for each line found + allocate(eqcopy(size(exitcomp))) +! This will be used for the second call to find_inv ?? +! call restore_constitutions(newnode%linehead(jj)%lineceq,exitcomp) + if(gx%bmperr.ne.0) goto 1000 +! -------------- +! Now we use FIND_INV to find a phase with zero amount together with +! either nodein(1) or nodein(2) and we come back again to find a second phase +! with zero amount with the other phase fix at the entering line + jj=2 + firstoutfix=0 +391 continue +! use the lineceq record to be used in the linehead for the calculation + eqcopy=exitcomp + tmpline=>newnode%linehead(jj) +! call list_sorted_phases(kou,tmpline%lineceq) +! write(*,*)'SMP before calling find_inv with: ',nodein +! the parir of phases for the line is returned in nodeut +! nodeut(1) is a new phase, nodeut(2) is either nodein(1) or nodein(2) + call find_inv_nodephase(nodeut,nodein,stabph,invph,6,axarr,& + eqcopy,tmpline) + nolinefix: if(gx%bmperr.ne.0) then +! Calculation failed, take any unused phase in invph to use as nodefix + gx%bmperr=0 + failed: do kk=1,stabph+2 + if(invph(4,kk).eq.0) then +! select this unused phase as linefix. Note invph(4,linefix) incremented below + nodeut(1)=kk endif - newnode%linehead(jj)%linefixph%ixphase=invph(1,tz) - newnode%linehead(jj)%linefixph%compset=invph(2,tz) - newnode%linehead(jj)%linefix_phr=invph(3,tz) - call get_phase_name(invph(1,tz),invph(2,tz),phases) - zp=len_trim(phases)+2 -! write(*,*)'smp2: ',trim(phases),tz,qy -! The LINEFIX phase is easy because each phase should appear twice as FIX -! but which to be NODEFIXPHASE? -! Rule: Each pair LINEFIX/NODEFIX should appear inverted. -! here set (arbitraily) the phase forbidden to become stable when line starts - klop2: do kk=jlast,1,-1 -! with same fix phase we must not have the same nodfixphase -! the same phase must not be forbidden more than twice .... -! write(*,440)kk,tz,qy,invph(4,kk) -440 format('smp2A forbidden: ',6i4) - if(kk.ne.tz .and. kk.ne.qy .and. & - invph(4,kk).lt.2) then - newnode%linehead(jj)%nodfixph=invph(3,kk); qy=kk; - invph(4,kk)=invph(4,kk)+1 - exit klop2 - endif - enddo klop2 - call get_phase_name(invph(1,qy),invph(2,qy),phases(zp:)) - zp=len_trim(phases)+2 -! write(*,*)'smp2: ',trim(phases),tz,qy -! now set the stable phases avoiding tz - zz=1 - do kk=1,stabph - ! maybe this if sequence is too simple ... - if(zz.eq.tz) zz=zz+1 - if(zz.eq.qy) zz=zz+1 - if(zz.eq.tz) zz=zz+1 - newnode%linehead(jj)%stableph(kk)%ixphase=invph(1,zz) - newnode%linehead(jj)%stableph(kk)%compset=invph(2,zz) - newnode%linehead(jj)%stablepham(kk)=1.0D-2 - newnode%linehead(jj)%stable_phr(kk)=invph(3,zz) - call get_phase_name(invph(1,zz),invph(2,zz),phases(zp:)) - zp=len_trim(phases)+2 - zz=zz+1 - enddo - write(*,430)jj,phases - if(mod(jj,2).ne.0) then -! increment jfix by one when jj is odd, i.e. 3, 5, etc and clear qy - jfix=jfix+1; qy=0 +! as nodefix use an unused invph(5,nodein(1) or nodein(2)) + if(invph(5,invph(5,nodein(1))).eq.1) then + nodeut(2)=nodein(1) + else + nodeut(2)=nodein(2) endif -! write(*,*)'smp2a jfix: ',jj,jfix,tz,qy - endif first - enddo invexit + exit nolinefix + enddo failed + write(*,*)'SMP cannot find exit line with nodefix:',nodefix + stop + endif nolinefix + linefix=nodeut(1) + nodefix=nodeut(2) +! write(*,'(a,6i4)')'SMP back from find_inv: ',jj,nodeut +! call list_sorted_phases(kou,tmpline%lineceq) + if(jj.eq.2) then + firstoutfix=nodeut(1) + endif +! write(*,'(a,i4,i2,5x,i5,i2)')'SMP phases at exit line: ',& +! invph(1,linefix),invph(2,linefix),invph(1,nodefix),invph(2,nodefix) +! doubline is used to jump back to label 392 for second line from same point + doubline=0 +! Here we will create 2 exit lines with linefix/nodefix in both positions +! We should remember the linefix phase to generate the next set of 2 lines +! repeat this with switched linefix/nodefix +392 continue + newnode%linehead(jj)%linefixph%ixphase=invph(1,linefix) + newnode%linehead(jj)%linefixph%compset=invph(2,linefix) + newnode%linehead(jj)%linefix_phr=invph(3,linefix) + invph(4,linefix)=invph(4,linefix)+1 + call get_phase_name(invph(1,linefix),invph(2,linefix),phases) + zp=len_trim(phases)+3 + phases(zp-1:zp-1)='(' + newnode%linehead(jj)%nodfixph=invph(3,nodefix) + invph(5,nodefix)=invph(5,nodefix)+1 + call get_phase_name(invph(1,nodefix),invph(2,nodefix),phases(zp:)) + zp=len_trim(phases)+3 + phases(zp-2:zp-2)=')' + kk=0 + names2: do zz=1,stabph+2 + if(zz.eq.linefix .or. zz.eq.nodefix) cycle names2 + kk=kk+1 + if(kk.le.stabph) then + newnode%linehead(jj)%stableph(kk)%ixphase=invph(1,zz) + newnode%linehead(jj)%stableph(kk)%compset=invph(2,zz) + newnode%linehead(jj)%stablepham(kk)=1.0D-2 + newnode%linehead(jj)%stable_phr(kk)=invph(3,zz) + else + write(*,'(a,10i5)')'SMP2 too many stable phases: ',jj,kk,zz,& + invph(1,zz),invph(2,zz),linefix,nodefix + endif + call get_phase_name(invph(1,zz),invph(2,zz),phases(zp:)) + zp=len_trim(phases)+2 + enddo names2 + write(*,430)jj,linefix,nodefix,trim(phases) + jj=jj+1 + if(doubline.eq.0) then +! we have switch linefix and nodefis to create one more exit line + doubline=linefix + linefix=nodefix; nodefix=doubline + phases=' ' +! copy also the phase amounts and constitutions to this lineceq + call restore_constitutions(newnode%linehead(jj)%lineceq,eqcopy) + goto 392 + endif + if(jj.eq.4) then +! If jj is 4 we generated 3 lines, we have to call find_inv again +! to find another phase fix with zero amount at the invariant +! write(*,*)'SMP generate 2 more exit lines' + goto 391 + endif + if(jj.eq.6) then +! If jj=6 we use the 2 phases we found above, firstoutfix and last nodeut(1) +! to generate the last 2 exit lines +! write(*,*)'SMP last exits:',firstoutfix,nodeut(1),linefix + nodefix=firstoutfix; linefix=nodeut(1); doubline=0 +! we cannot set any constitions for these lines + goto 392 + endif +! clean up ... + deallocate(eqcopy) + deallocate(exitcomp) +!--------------------------------------- +! when doubline is 4 all is done! +! stop 'this will work' +! error exit from above 490 continue ! stop ' *** Unfinished invariant isopleth node exits *** ' - end select + end select exits !========================================================================= goto 1000 !------------------------------------------- @@ -5288,10 +5500,10 @@ end subroutine map_newnode ! redefined argument mecreq to mecreq1 !\addtotable subroutine sort_invph !\begin{verbatim} - subroutine sort_invph(nitems,array) + subroutine sort_invph(nitems,ndim,array) ! primitive sorting of array implicit none - integer nitems,array(4,*) + integer nitems,ndim,array(ndim,*) !\end{verbatim} ! sort array in acending order of value in array(3,*) integer ia,ib,ic,more @@ -5302,7 +5514,7 @@ subroutine sort_invph(nitems,array) do ia=2,nitems if(array(3,ia-1).gt.array(3,ia)) then more=more+1 - do ic=1,4 + do ic=1,ndim ib=array(ic,ia-1); array(ic,ia-1)=array(ic,ia); array(ic,ia)=ib enddo endif @@ -5313,6 +5525,465 @@ end subroutine sort_invph !/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ +!\addtotable subroutine find_inv_nodephase_old +!\begin{verbatim} + subroutine find_inv_nodephase_old(phnew,phnode,stabph,invph,dim1,axarr,& + eqcopy,linerec) +! Find the phase to set as linefix when nodefix has zero amount +! phnew is 2 phase index at point of entering the invariant +! phnode is two phase indices for a point of exit from the invariant +! stabph is number of stable phases along the lines, at node point 2 more stable +! invph is matrix with all phases at the node, dim1 is its first dim +! axarr has axis information (needed to maniplulate conditions) +! eqcopy is array with phase anounts and constitutions at new point +! linerec is a line record with all necessary data to calculate en equil + implicit none + integer phnew,phnode,stabph,dim1,invph(dim1,*) + type(map_axis), dimension(*) :: axarr + double precision eqcopy(*) + type(map_line), pointer :: linerec +!\end{verbatim} +! linefix and nodefix is 2nd index in invph(5,*) + type(gtp_condition), pointer :: lastcond,pcond,pcondx,pcondt + type(gtp_equilibrium_data), pointer :: thisceq +! integer, parameter :: inmap=1 +! double precision, allocatable, dimension(:) :: eqcopy +! type(meq_setup), pointer :: meqrec +! not really needed + integer linefix,nodefix,reset + integer ii,jj,kk,nodeph,mapx,iadd,irem,iz,jax,mode,okcond,lokph,lokcs,lokcs2 + character*24 phname1,phname2,phname3 + double precision, parameter :: mamfu=1.0D-6 +! +! we must set the axis condtion on T and +! remove the axis condition on the composition on the other axis +!-------------------- copied from somewhere ................... +! write(*,10)phnew,phnode,stabph,linerec%lineceq%tpval(1) +10 format(/'SMP find_inv generating an exit from isopleth invariant',3i4,F10.2) +! call list_conditions(kou,linerec%lineceq) + reset=globaldata%status +! supress messages from calceq3 done inside find_inv + globaldata%status=ibset(globaldata%status,GSSILENT) + thisceq=>linerec%lineceq +! loop both axis + okcond=0 + do jax=1,2 +! Set the condition on T and remember the condition on composition + lastcond=>thisceq%lastcondition + if(.not.associated(lastcond)) then + write(*,*)'in find_inv, no conditions: ',jax + gx%bmperr=4221; goto 1000 + endif + pcond=>lastcond +60 continue + pcond=>pcond%next + if(pcond%seqz.eq.axarr(jax)%seqz) goto 70 + if(.not.associated(pcond,lastcond)) goto 60 + write(*,*)'in find_inv the axis condition not found: ',jax + gx%bmperr=4221; goto 1000 +! +70 continue + if(pcond%statev.ge.10) then +! save pointer to extensive condition and remove it + pcond%active=1 + okcond=okcond+1 + pcondx=>pcond + elseif(pcond%statev.eq.1) then +! set condition on T as active + okcond=okcond+1 + pcondt=>pcond + pcond%active=0 + endif + enddo + if(okcond.ne.2) then + write(*,*)'Conditions not T and X, quitting' + gx%bmperr=4399 + goto 1000 + endif +!----------------- end copy from somewhere................. +! write(*,*)'SMP find_inv conditions with no fix phase:' +! call list_conditions(kou,thisceq) + linefix=phnew + nodefix=phnode +! extract the name of the linefix and nodefix phases + call get_phase_name(invph(1,nodefix),invph(2,nodefix),phname1) + call get_phase_name(invph(1,linefix),invph(2,linefix),phname2) +! suspend the linefix phase + call change_many_phase_status(phname2,PHSUS,zero,thisceq) +! set nodefix fix with zero amount + call change_many_phase_status(phname1,PHFIXED,zero,thisceq) + if(gx%bmperr.ne.0) then + write(*,*)'SMP error suspending or fixing ',& + trim(phname1)//' or '//trim(phname2) + gx%bmperr=4399; goto 1000 + endif +! Calculate an equilibrium at the invariant T with nodefix as fix + write(*,*)'SMP find_inv conditions before calling calceq3' + call list_conditions(kou,thisceq) +! mode=0 no grid minimizer and no data structures + mode=0 + call calceq3(mode,.FALSE.,thisceq) + if(gx%bmperr.ne.0) then + write(*,*)'SMP fail calculate invariant',gx%bmperr + goto 1000 + endif +! debug listing + call list_sorted_phases(kou,thisceq) + write(*,*)'SMP find_inv: phases at invariant' +! set nodefix entered with zero amount and linefix as enetered with small + call change_many_phase_status(phname1,PHENTERED,zero,thisceq) + call change_many_phase_status(phname2,PHENTERED,1.0D-2,thisceq) + if(gx%bmperr.ne.0) then + write(*,*)'SMP error restoring phases ',trim(phname1//phname2) + gx%bmperr=4399; goto 1000 + endif +! remove the condition on compostion + pcondx%active=1 +!------------------------------------------------------- +! start loop to find new linefix phase +! deallocate=.false. +! save current amount of phases and constitutions +! call save_constitutions(thisceq,eqcopy) +! if(gx%bmperr.ne.0) goto 1000 +! we must check that the amount of the nodefix phase is zero !! + call get_phase_compset(invph(1,nodefix),invph(2,nodefix),lokph,lokcs) + if(gx%bmperr.ne.0) then + write(*,*)'SMP find_inv failed get lokcs' + goto 1000 + endif +! and the amount of old linefix must not be zero at the same time! + call get_phase_compset(invph(1,linefix),invph(2,linefix),lokph,lokcs2) + if(gx%bmperr.ne.0) then + write(*,*)'SMP find_inv failed get lokcs' + goto 1000 + endif +! suspend/restore phases not involved in the invariant +! maybe not needed ? +! call suspend_somephases(1,invph,6,stabph+2,thisceq) +! this is to restore +! call suspend_somephases(0,invph,6,stabph+2,thisceq) +! +! THIS IS OLD VERSION !! +! +! all calculations below are made at fixed T with different fix phases +! afterwards we check that the amount of nodefix is still zero (or very small) +! total number of phases at nodepoint is stabph+2 + loop: do ii=1,stabph+2 +! inverted order ... works for map7 but fails for map16 SUCK +! loop: do ii=stabph+2,1,-1 + if(invph(4,ii).eq.0) then +! only phases with invph(4,ii)=0 can be tested as linefix + call get_phase_name(invph(1,ii),invph(2,ii),phname3) +! write(*,*)'SMP find_inv testing: ',ii,': ',trim(phname1//phname3) +! test ii as new linefix, set it fix with zero amount + call change_many_phase_status(phname3,PHFIXED,zero,thisceq) + call list_conditions(kou,thisceq) +! debug listing + call list_sorted_phases(kou,thisceq) + write(*,*)'SMP find_inv: phases before calculations' + mode=0 + call calceq3(mode,.FALSE.,thisceq) + if(gx%bmperr.ne.0) then +! calculation error, remove phase ii as fix and try another + write(*,*)'SMP find_inv error: ',gx%bmperr + gx%bmperr=0 + call change_many_phase_status(phname3,PHENTERED,zero,& + thisceq) + call restore_constitutions(thisceq,eqcopy) + if(gx%bmperr.ne.0) goto 1000 + cycle loop + endif +! Calculation converged, check if amount of nodefix, it should be zero + write(*,*)'SMP find_inv amfu: ',nodefix,lokcs,& + thisceq%phase_varres(lokcs)%amfu + if(thisceq%phase_varres(lokcs)%amfu.gt.mamfu) then +! amount of nodefix phase too large, remove phase ii as fix and try another + call change_many_phase_status(phname3,PHENTERED,zero,& + thisceq) + call restore_constitutions(thisceq,eqcopy) + if(gx%bmperr.ne.0) goto 1000 + cycle loop + endif + if(thisceq%phase_varres(lokcs2)%amfu.lt.mamfu) then +! the amount of old linefix phase must be larger then 0! + write(*,*)'SMP both linefix and nodefix must not be zero!' + call change_many_phase_status(phname3,PHENTERED,zero,& + thisceq) + call restore_constitutions(thisceq,eqcopy) + if(gx%bmperr.ne.0) goto 1000 + cycle loop + endif + if(invph(1,ii).eq.invph(1,nodefix)) then +! this is a problem for map7 macro ... + write(*,*)'SMP ignore 2nd composition set of fix: ',ii,nodefix + call change_many_phase_status(phname3,PHENTERED,zero,& + thisceq) + call restore_constitutions(thisceq,eqcopy) + write(*,*)'SMP error? ',gx%bmperr,ii,stabph+2 + if(gx%bmperr.ne.0) then + gx%bmperr=0 + endif + cycle loop + endif +! Heureka! we have found a linefix phase, clean up the conditions +! if we found a phase which is a second composition set of a phase +! that is set fix, ignore!! (like FCC and FCC-carbide + write(*,*)'SMP find_inv success, linefix: ',phname3,ii,nodefix + phnew=ii; goto 1000 + endif + enddo loop +! we did not find any linefix, return error and use default + write(*,*)'SMP find_inv failed to find phase',ii + gx%bmperr=4399 +1000 continue + ii=gx%bmperr; gx%bmperr=0 +! restore axis conditions, set x condition + pcondx%active=0 +! remove T condition + pcondt%active=1 +! +! write(*,*)'SMP finc_inv exit conditions (same as on entering?):' +! call list_conditions(kou,thisceq) +! write(*,*) + gx%bmperr=ii +! remove quite mode! + globaldata%status=reset + return + end subroutine find_inv_nodephase_old + +!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ + +!\addtotable subroutine find_inv_nodephase +!\begin{verbatim} + subroutine find_inv_nodephase(phut,phin,stabph,invph,dim1,axarr,& + eqcopy,linerec) +! Find the phase to set as linefix when nodefix has zero amount +! phut will on exit be two phases with zero amount at the invariant +! phin is on enter the two phase with zero amount when the invariant is found +! stabph is number of stable phases along the lines, at node point 2 more stable +! invph is matrix with all phases at the node, dim1 is its first dim +! axarr has axis information (needed to maniplulate conditions) +! eqcopy is array with phase anounts and constitutions at new point +! linerec is a line record with all necessary data to calculate en equil + implicit none + integer phut(2),phin(2),stabph,dim1,invph(dim1,*) + type(map_axis), dimension(*) :: axarr + double precision eqcopy(*) + type(map_line), pointer :: linerec +!\end{verbatim} +! linefix and nodefix is 2nd index in invph(5,*) + type(gtp_condition), pointer :: lastcond,pcond,pcondx,pcondt + type(gtp_equilibrium_data), pointer :: thisceq +! integer, parameter :: inmap=1 +! double precision, allocatable, dimension(:) :: eqcopy +! type(meq_setup), pointer :: meqrec +! not really needed + integer linefix,nodefix,reset,ph1,ph2,lokcs1,lokcs2,enteredphase + integer ii,jj,kk,nodeph,mapx,iadd,irem,iz,jax,mode,okcond,lokph + character*24 phname1,phname2,phname3 + double precision, parameter :: mamfu=1.0D-6 +! +! we must set the axis condition on T and +! remove the axis condition on the composition on the other axis +!-------------------- copied from somewhere ................... +! write(*,10)phin,stabph,linerec%lineceq%tpval(1) +10 format(/'SMP find_inv exits from isopleth invariant: ',2i4,2x,i2,F10.2) + reset=globaldata%status +! supress messages from calceq3 done inside find_inv +! globaldata%status=ibset(globaldata%status,GSSILENT) + thisceq=>linerec%lineceq +! constitution is OK here +! call list_conditions(kou,linerec%lineceq) +! call list_sorted_phases(kou,thisceq) +! write(*,*)'SMP constitution entering find_inv' +! suspend all phases not involved in the invariant, 1 means suspend + call suspend_somephases(1,invph,6,stabph+2,thisceq) + if(gx%bmperr.ne.0) then + write(*,*)'SMP error calling suspend_somephases'; goto 1000 + endif +! loop both axis to extract condition pointer + okcond=0 + do jax=1,2 +! Set the condition on T and remember the condition on composition + lastcond=>thisceq%lastcondition + if(.not.associated(lastcond)) then + write(*,*)'in find_inv, no conditions: ',jax + gx%bmperr=4221; goto 1000 + endif + pcond=>lastcond +60 continue + pcond=>pcond%next + if(pcond%seqz.eq.axarr(jax)%seqz) goto 70 + if(.not.associated(pcond,lastcond)) goto 60 + write(*,*)'in find_inv the axis condition not found: ',jax + gx%bmperr=4221; goto 1000 +! +70 continue + if(pcond%statev.ge.10) then +! save pointer to extensive condition and remove it + pcond%active=1 + okcond=okcond+1 + pcondx=>pcond + elseif(pcond%statev.eq.1) then +! set condition on T as active + okcond=okcond+1 + pcondt=>pcond +! try setting two fix phases .... remove condition on T +! pcond%active=1 + endif + enddo + if(okcond.ne.2) then + write(*,*)'Conditions not T and X, quitting' + gx%bmperr=4399 + goto 1000 + endif +!----------------- end copy from somewhere................. +! most of these variables are just for debugging + ph1=phin(1) + ph2=phin(2) +! extract the name of the linefix and nodefix phases + call get_phase_name(invph(1,ph1),invph(2,ph1),phname1) + call get_phase_name(invph(1,ph2),invph(2,ph2),phname2) +! set small amounts of ph1 ?? + call change_many_phase_status(phname1,PHENTSTAB,mamfu,thisceq) + call change_many_phase_status(phname1,PHENTSTAB,zero,thisceq) + if(gx%bmperr.ne.0) then + write(*,*)'SMP error setting small amounts of line/nodefix' + goto 1000 + endif +! we check the amounts and driving forces of both linefix and nodefix + call get_phase_compset(invph(1,ph1),invph(2,ph1),lokph,lokcs1) + call get_phase_compset(invph(1,ph2),invph(2,ph2),lokph,lokcs2) + if(gx%bmperr.ne.0) then + write(*,*)'SMP find_inv failed get phase_varres index' + goto 1000 + endif +! remove the condition on compostion and set fix T + pcondx%active=1 + pcondt%active=0 +! set phname1 as fix with zero amount +! call change_many_phase_status(phname1,PHFIX,zero,thisceq) +! if(gx%bmperr.ne.0) then +! write(*,*)'SMP error setting '//trim(phname1)//' as fix' +! goto 1000 +! endif +!------------------------------------------------------- +! start loop to find a phase with zero amount with either phin +! deallocate=.false. +! save current amount of phases and constitutions +! call save_constitutions(thisceq,eqcopy) +! if(gx%bmperr.ne.0) goto 1000 +! +! THIS IS NEW VERSION +! +! the loop below if first run with phases1 entered with small amount +! if that failes we jump hack to label 50 with phase2 entered with small amount + enteredphase=1 +50 continue +! all calculations below are made at fixed T with different fix phases +! afterwards we check that the amount of nodefix is still zero (or very small) +! total number of phases at nodepoint is stabph+2 + loop: do ii=1,stabph+2 + if(invph(4,ii).eq.0 .and. invph(5,ii).eq.0) then +! only test phases which has not already been used, extract its name + call restore_constitutions(thisceq,eqcopy) + call get_phase_name(invph(1,ii),invph(2,ii),phname3) +! write(*,76)trim(phname3),ii +76 format(/'SMP find_inv testing: ',a,i4) +! set ii fix with zero amount + call change_many_phase_status(phname3,PHFIXED,zero,thisceq) +! call list_conditions(kou,thisceq) + if(enteredphase.eq.1) then +! set a small amount on phase1 (if both phases small error TOO MANY STABLE PHA) + call change_many_phase_status(phname1,PHENTSTAB,1.0D-3,thisceq) + else +! set a small amount on phase2 +! write(*,*)'Trying harder with ',phname2 + call change_many_phase_status(phname2,PHENTSTAB,1.0D-3,thisceq) + endif +! debug listing +! call list_sorted_phases(kou,thisceq) +! write(*,*)'SMP find_inv: phases BEFORE calculations',gx%bmperr + mode=0 + call calceq3(mode,.FALSE.,thisceq) + if(gx%bmperr.ne.0) then +! calculation error, remove phase ii as fix and try another +! write(*,*)'SMP find_inv error: ',gx%bmperr + gx%bmperr=0 + call change_many_phase_status(phname3,PHENTERED,zero,thisceq) + if(gx%bmperr.ne.0) goto 1000 + cycle loop + endif +! debug listing +! call list_sorted_phases(kou,thisceq) +! write(*,*)'SMP find_inv: phases AFTER calculations',gx%bmperr +! there should be stabph-1 stable phases + jj=0 + do kk=1,stabph+2 + if(thisceq%phase_varres(invph(6,kk))%amfu.gt.zero) jj=jj+1 + enddo +! write(*,*)'SMP number of stable phases',jj,stabph + if(jj.ne.stabph) goto 120 +! Calculation converged, check if the amount of phin(1) or phin(2) is zero +! write(*,111)'SMP find_inv amfu: ',& +! trim(phname1),invph(4,phin(1)),thisceq%phase_varres(lokcs1)%amfu,& +! trim(phname2),invph(4,phin(2)),thisceq%phase_varres(lokcs2)%amfu +111 format(a,a,i2,1pe12.4,5x,a,i2,1pe12.4) + if(invph(4,phin(1)).eq.1 .and.& + thisceq%phase_varres(lokcs1)%amfu.lt.mamfu) then +! phase invph(*,ii) and invph(*,phin(1)) have zero amount, use as exit +! +! ADD CHECK IF amfu for BOTH phin(1) and ph(2) are zero we must check dgm ... +! + phut(1)=ii; phut(2)=phin(1) +! write(*,112)trim(phname3)//'+'//trim(phname1),phut(1),phut(2) +112 format('SMP find_inv **** success: ',a,2i4) + goto 200 + elseif(invph(4,phin(2)).eq.1 .and. & + thisceq%phase_varres(lokcs2)%amfu.lt.mamfu) then +! phase invph(*,ii) and invph(*,phin(1)) have zero amount, use as exit + phut(1)=ii; phut(2)=phin(2) +! write(*,112)trim(phname3)//'+'//trim(phname2),phut(1),phut(2) + goto 200 + endif +! else +! write(*,*)'SMP skipping ',ii + endif +120 continue +! set phase as entered and restore constitutions for another phase as fix + call change_many_phase_status(phname3,PHENTERED,zero,thisceq) +! call restore_constitutions(thisceq,eqcopy) + if(gx%bmperr.ne.0) goto 1000 + enddo loop +! we have not found any set of phases for an exit line + if(enteredphase.eq.1) then + enteredphase=2; goto 50 + endif + ! write(*,*)'SMP find_inv failed to find two phases with zero amount' + gx%bmperr=4399; goto 1000 +200 continue +! we have a pair of phases in phut, reset phname3 as entered + call change_many_phase_status(phname3,PHENTERED,zero,thisceq) +!----------------------------- exit +1000 continue + ii=gx%bmperr; gx%bmperr=0 +! restore axis conditions, set x condition + pcondx%active=0 +! remove T condition + pcondt%active=1 +! this is to restore phases earlier suspended, 0 menas restore + call suspend_somephases(0,invph,6,stabph+2,thisceq) + if(gx%bmperr.ne.0) then + write(*,*)'SMP error calling suspend_somephases'; goto 1000 + endif + gx%bmperr=ii +! remove quite mode! + globaldata%status=reset + return + end subroutine find_inv_nodephase + +!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ + !\addtotable subroutine reserve_saveceq !\begin{verbatim} subroutine reserve_saveceq(location,saveceq) @@ -5473,6 +6144,7 @@ subroutine map_findline(maptop,axarr,mapfix,mapline) allocate(mapfix) ! with only 2 axis we have just 1 fix phase for mapping, fixph is a tuple allocate(mapfix%fixph(1)) + mapfix%status=0 mapfix%nfixph=1 mapfix%fixph=mapline%linefixph(1) ! we can have several stable phases when no tie-lines in plane @@ -5483,7 +6155,7 @@ subroutine map_findline(maptop,axarr,mapfix,mapline) ! write(*,*)'Findline: Tie-lines not in plane: ',nyline,ip ! create a heading text for the line phaseset=' ' - call get_phasetup_name_old(mapfix%fixph(1),phaseset) + call get_phasetuple_name(mapfix%fixph(1),phaseset) if(gx%bmperr.ne.0) goto 1000 ip=len_trim(phaseset) phaseset(ip+1:ip+10)=', stable: ' @@ -5502,7 +6174,7 @@ subroutine map_findline(maptop,axarr,mapfix,mapline) ! this is stored only for "real" nodes mapfix%stableph(jp)=mapnode%linehead(nyline)%stableph(jp) mapfix%stable_phr(jp)=mapnode%linehead(nyline)%stable_phr(jp) - call get_phasetup_name_old(mapfix%stableph(jp),phaseset(ip:)) + call get_phasetuple_name(mapfix%stableph(jp),phaseset(ip:)) if(gx%bmperr.ne.0) goto 1000 ! this values hould perhaps be in linehead?? ! mapfix%stablepham(jp)=mapnode%linehead(nyline)%stablepham(jp) @@ -5544,6 +6216,7 @@ subroutine map_findline(maptop,axarr,mapfix,mapline) allocate(mapfix%stablepham(1)) allocate(mapfix%stable_phr(1)) mapfix%nfixph=1 + mapfix%status=0 !......................................................... ! trying to impove mapping with two extensive axis variables nofecond=0 @@ -5679,7 +6352,7 @@ subroutine map_findline(maptop,axarr,mapfix,mapline) endif ! create a heading text for the line phaseset=' ' - call get_phasetup_name_old(mapfix%fixph(1),phaseset) + call get_phasetuple_name(mapfix%fixph(1),phaseset) if(gx%bmperr.ne.0) goto 1000 ip=len_trim(phaseset)+4 phaseset(ip-2:ip-2)='+' @@ -5695,7 +6368,7 @@ subroutine map_findline(maptop,axarr,mapfix,mapline) mapfix%stableph(1)=mapnode%linehead(nyline)%stableph(1) mapfix%stable_phr(1)=mapnode%linehead(nyline)%stable_phr(1) endif - call get_phasetup_name_old(mapfix%stableph(1),phaseset(ip:)) + call get_phasetuple_name(mapfix%stableph(1),phaseset(ip:)) if(gx%bmperr.ne.0) goto 1000 ! set positive amount both in mapfix and in phase_varres ...?? mapfix%stablepham(1)=one-fixpham @@ -5737,7 +6410,7 @@ subroutine map_findline(maptop,axarr,mapfix,mapline) phaseset=' ' ip=1 do jp=1,mapnode%linehead(1)%nstabph - call get_phasetup_name_old(mapnode%linehead(1)%stableph(jp),& + call get_phasetuple_name(mapnode%linehead(1)%stableph(jp),& phaseset(ip:)) if(gx%bmperr.ne.0) goto 1000 ip=len_trim(phaseset)+2 @@ -5754,8 +6427,19 @@ subroutine map_findline(maptop,axarr,mapfix,mapline) write(*,*)'Line with unkonwn stable phases: ',& mapnode%linehead(1)%nstabph endif -! for step calculation mapfix is not needed -! nullify(mapfix) +! write(*,*)'SMP is mapfix allocated? ',allocated(mapfix) +! if(.not.allocated(mapfix)) then +! for STEP calculations mapfix was normally not allocated but I need the status +! but instead of adding this set a bit in the meqrec record after first +! call to calceq7 +! allocate(mapfix) +! mapfix%nfixph=0 +! mapfix%status=0 +! if(btest(mapnode%status,STEPINVARIANT)) then +! write(*,*)'SMP invarant step node',mapnode%status +! mapfix%status=ibset(mapfix%status,STEPINVARIANT) +! endif +! endif endif 1000 continue return @@ -5896,6 +6580,7 @@ subroutine map_findline_old(maptop,axarr,mapfix,mapline) deallocate(mapfix) endif allocate(mapfix) + mapfix%status=0 ! with only 2 axis we have just 1 fix phase for mapping allocate(mapfix%fixph(1)) mapfix%nfixph=1 @@ -5908,7 +6593,7 @@ subroutine map_findline_old(maptop,axarr,mapfix,mapline) ! write(*,*)'Findline: Tie-lines not in plane: ',nyline,ip ! create a heading text for the line phaseset=' ' - call get_phasetup_name_old(mapfix%fixph(1),phaseset) + call get_phasetuple_name(mapfix%fixph(1),phaseset) if(gx%bmperr.ne.0) goto 1000 ip=len_trim(phaseset)+4 phaseset(ip-1:ip-2)='+' @@ -5926,7 +6611,7 @@ subroutine map_findline_old(maptop,axarr,mapfix,mapline) ! this is stored only for "real" nodes mapfix%stableph(jp)=mapnode%linehead(nyline)%stableph(jp) mapfix%stable_phr(jp)=mapnode%linehead(nyline)%stable_phr(jp) - call get_phasetup_name_old(mapfix%stableph(jp),phaseset(ip:)) + call get_phasetuple_name(mapfix%stableph(jp),phaseset(ip:)) if(gx%bmperr.ne.0) goto 1000 ! this values hould perhaps be in linehead?? mapfix%stablepham(jp)=one @@ -5967,6 +6652,7 @@ subroutine map_findline_old(maptop,axarr,mapfix,mapline) allocate(mapfix%stablepham(1)) allocate(mapfix%stable_phr(1)) mapfix%nfixph=1 + mapfix%status=0 !......................................................... ! trying to impove mapping with two extensive axis variables nofecond=0 @@ -6102,7 +6788,7 @@ subroutine map_findline_old(maptop,axarr,mapfix,mapline) endif ! create a heading text for the line phaseset=' ' - call get_phasetup_name_old(mapfix%fixph(1),phaseset) + call get_phasetuple_name(mapfix%fixph(1),phaseset) if(gx%bmperr.ne.0) goto 1000 ip=len_trim(phaseset)+4 phaseset(ip-2:ip-2)='+' @@ -6118,7 +6804,7 @@ subroutine map_findline_old(maptop,axarr,mapfix,mapline) mapfix%stableph(1)=mapnode%linehead(nyline)%stableph(1) mapfix%stable_phr(1)=mapnode%linehead(nyline)%stable_phr(1) endif - call get_phasetup_name_old(mapfix%stableph(1),phaseset(ip:)) + call get_phasetuple_name(mapfix%stableph(1),phaseset(ip:)) if(gx%bmperr.ne.0) goto 1000 ! set positive amount both in mapfix and in phase_varres ...?? mapfix%stablepham(1)=one-fixpham @@ -6160,7 +6846,7 @@ subroutine map_findline_old(maptop,axarr,mapfix,mapline) phaseset=' ' ip=1 do jp=1,mapnode%linehead(1)%nstabph - call get_phasetup_name_old(mapnode%linehead(1)%stableph(jp),& + call get_phasetuple_name(mapnode%linehead(1)%stableph(jp),& phaseset(ip:)) if(gx%bmperr.ne.0) goto 1000 ip=len_trim(phaseset)+2 @@ -6336,78 +7022,6 @@ end function tieline_inplane !/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ -!\addtotable logical function inveq -!\begin{verbatim} - logical function inveq(phases,ceq) -! Only called for tie-lines not in plane. If tie-lines in plane then all -! nodes are invariants. -! UNFINISHED - integer phases - type(gtp_equilibrium_data), pointer :: ceq -!\end{verbatim} - integer nrel,ii,nostph,tpvar,degf,www - type(gtp_condition), pointer :: pcond,lastcond - type(gtp_state_variable), pointer :: stvr -! How to know if the ceq is invariant? Gibbs phase rule, Degrees of freedom -! f = n + z - w - p -! where n is number of components, z=2 if T and P variable, -! z=1 if T or P variable, z=0 if both T and P fixed, -! w is number of other potential conditions (MU, AC) -! p is number of stable phases. -! write(*,*)'in inveq' - nrel=noel() -! sum up nubler of stable phases and check if T and P are fixed - nostph=0 -! ntups=nooftup() -! do ii=1,noofphasetuples() - do ii=1,nooftup() - if(ceq%phase_varres(phasetuple(ii)%lokvares)%phstate.gt.0) & - nostph=nostph+1 - enddo -! loop all conditions - lastcond=>ceq%lastcondition - pcond=>lastcond - tpvar=2 - www=0 -100 continue - if(pcond%active.eq.0) then -! condtion is active - stvr=>pcond%statvar(1) -! statevarid 1 is T and 2 is P - if(stvr%statevarid.eq.1 .or. stvr%statevarid.eq.2) then -! Hm, ceq is not the equilibrium record for the node point ... - tpvar=tpvar-1 - elseif(stvr%statevarid.lt.10) then -! potential/activity condition for a component - www=www+1 - endif - endif - pcond=>pcond%next - if(.not.associated(pcond,lastcond)) goto 100 -! -! Hm again, ignore tpvar? -! degf=nrel+tpvar-www-nostph - degf=nrel-www-nostph -! write(*,*)'in inveq 2',nrel,tpvar,www,nostph,degf - if(degf.lt.0) then -! We have an invariant equilibrium, return the number of stable phases - phases=nostph - inveq=.true. -! write(*,200)'We have an invariant equilibrium!',nrel,tpvar,nostph,phases -!200 format(a,5i7) - else -! write(*,210)degef,nrel,tpvar,phases -!210 format('Not inveq, elements, stable phases: ',4i4) -! if not invariant isoplet node there are 3 exits (2 lines crossing) - phases=nostph - inveq=.false. - endif -1000 continue - return - end function inveq - -!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\!/!\ - !\addtotable subroutine list_map_equilibrium !\begin{verbatim} subroutine list_map_equilibrium(maptop,mapline,axarr,xxx,typ) @@ -6859,7 +7473,7 @@ subroutine step_separate(maptop,noofaxis,axarr,seqxyz,starteq) ! this call calculates the value of the axis condition with default composition call state_variable_val(svr,val,ceq) if(gx%bmperr.ne.0) goto 500 - call get_phasetup_name_old(entphcs(itup),name) + call get_phasetuple_name(entphcs(itup),name) ! axis variable is composition, skip hases with no variance ! call get_phase_variance(entphcs(itup)%phaseix,nv) call get_phase_variance(entphcs(itup)%ixphase,nv) diff --git a/stepmapplot/smp2B.F90 b/stepmapplot/smp2B.F90 index b095e4a..33e1dfb 100644 --- a/stepmapplot/smp2B.F90 +++ b/stepmapplot/smp2B.F90 @@ -1246,7 +1246,7 @@ subroutine ocplot2B(np,nrv,nlinesep,linesep,pltax,xax,anpax,anpdim,anp,lid,& 'set size ',F8.4', ',F8.4/& 'set xlabel "',a,'"'/'set ylabel "',a,'"'/& ! Help with stackoverflow to fix nice logo independent of plot size! - 'set label "~O{.0 C}" at graph -0.1, -0.1 font "Garamond Bold,20"'/& + 'set label "~O{.0 C}" at graph -0.1, -0.1 font "Garamond Bold,28"'/& 'set key ',a/& 'set linetype ',i2,' lc rgb "#000000" lw 2 pt 10'/& 'set linetype ',i2,' lc rgb "#4169E1" lw 2 pt 6'/& @@ -1269,7 +1269,7 @@ subroutine ocplot2B(np,nrv,nlinesep,linesep,pltax,xax,anpax,anpdim,anp,lid,& 'set size ',F8.4', ',F8.4/& 'set xlabel "',a,'"'/'set ylabel "',a,'"'/& ! Help with stackoverflow to fix nice logo independent of plot size! - 'set label "~O{.0 C}" at graph -0.1, -0.1 font "Garamond Bold,20"'/& + 'set label "~O{.0 C}" at graph -0.1, -0.1 font "Garamond Bold,28"'/& 'set key ',a/& 'set style line ',i2,' lt ',i2,' lc rgb "#000000" lw 2 pt 10'/& 'set style line ',i2,' lt ',i2,' lc rgb "#4169E1" lw 2 pt 6'/& @@ -1288,7 +1288,7 @@ subroutine ocplot2B(np,nrv,nlinesep,linesep,pltax,xax,anpax,anpdim,anp,lid,& write(21,8000) 8000 format(/'# Some useful GNUPLOT commands for editing the figure'/& '# This is a dashed line (on pdf/wxt):'/& - '# set style line 15 lt 0 lc rgb "C8C800" lw 2 pt 2'//& + '# set style line 15 lt 0 lc rgb "#C8C800" lw 2 pt 2'//& '# set pointsize 0.6'/& '# set label "text" at 0.5, 0.5 rotate by 60 font "arial,12"'/& '# set xrange [0.5 : 0.7] '/& @@ -1468,8 +1468,10 @@ subroutine ocplot2B(np,nrv,nlinesep,linesep,pltax,xax,anpax,anpdim,anp,lid,& ii=len_trim(graphopt%lowerleftcorner) if(ii.gt.0) then ! in square diagram below figure +! write(21,209)trim(graphopt%lowerleftcorner) +!209 format('set label "',a,'" at graph -0.10, -0.08 ') write(21,209)trim(graphopt%lowerleftcorner) -209 format('set label "',a,'" at graph -0.10, -0.08 ') +209 format('set label "',a,'" at graph -0.05, -0.05 ') endif ! if lowerleftcorner is empty ignore it !--------------------------------------------------------------- @@ -2271,21 +2273,15 @@ subroutine ocplot3B(same,nofinv,lineends,nx1,xval,ny1,yval,nz1,zval,plotkod,& ! THIS IS THE Y-AXIS WITH 60 degrees angle write(21,131)trim(pltax(2)), 0.15*xmax, 0.37*xmax 131 format('set label "',a,'" at ',F8.4,',',F8.4,' rotate by 60 '/& -! 'set label "O" at screen 0.130, 0.027 font "Garamond Bold,20"'/& -! 'set label "C" at screen 0.139, 0.027 font "Garamond Bold,20"') ! Help with stackoverflow to fix nice logo independent of plot size! - 'set label "~O{.0 C}" at graph -0.1, -0.1 font "Garamond Bold,20"') -! 'set label "O" at graph -0.090, -0.100 font "Garamond Bold,20"'/& -! 'set label "C" at graph -0.077, -0.100 font "Garamond Bold,20"') + 'set label "~O{.0 C}" at graph -0.1, -0.1 font "Garamond Bold,28"') ! we should also enforce same length of X and Y axis !!! else ! SQUARE DIAGRAM write(21,132)trim(pltax(2)) 132 format('set ylabel "',a,'"'/& ! Help with stackoverflow to fix nice logo independent of plot size! - 'set label "~O{.0 C}" at graph -0.1, -0.1 font "Garamond Bold,20"') -! 'set label "O" at graph -0.090, -0.100 font "Garamond Bold,20"'/& -! 'set label "C" at graph -0.080, -0.100 font "Garamond Bold,20"') + 'set label "~O{.0 C}" at graph -0.1, -0.1 font "Garamond Bold,28"') endif lz=graphopt%linetype write(21,133)lz,lz,lz,lz,lz,lz,lz,lz,lz,lz @@ -2305,7 +2301,7 @@ subroutine ocplot3B(same,nofinv,lineends,nx1,xval,ny1,yval,nz1,zval,plotkod,& write(21,8000) 8000 format(/'# Some useful GNUPLOT commands for editing the figure'/& '# This is a dashed line (on pdf/wxt):'/& - '# set style line 15 lt 0 lc rgb "C8C800" lw 2 pt 2'//& + '# set style line 15 lt 0 lc rgb "#C8C800" lw 2 pt 2'//& '# set pointsize 0.6'/& '# set label "text" at 0.5, 0.5 rotate by 60 font "arial,12"'/& '# set xrange [0.5 : 0.7] '/& @@ -2655,8 +2651,10 @@ subroutine ocplot3B(same,nofinv,lineends,nx1,xval,ny1,yval,nz1,zval,plotkod,& ii=len_trim(graphopt%lowerleftcorner) if(graphopt%gibbstriangle) then if(ii.gt.3) then - write(21,208)trim(graphopt%lowerleftcorner),-0.14 -208 format('set label "',a,'" at graph ',F10.4,', -0.05 ') +! write(21,208)trim(graphopt%lowerleftcorner),-0.14 +!208 format('set label "',a,'" at graph ',F10.4,', -0.05 ') + write(21,208)trim(graphopt%lowerleftcorner),-0.05 +208 format('set label "',a,'" at graph ',F10.4,', -0.03 ') elseif(ii.gt.0) then write(21,208)trim(graphopt%lowerleftcorner),-0.08 endif diff --git a/userif/pmon6.F90 b/userif/pmon6.F90 index 7567726..129d78b 100644 --- a/userif/pmon6.F90 +++ b/userif/pmon6.F90 @@ -310,7 +310,7 @@ subroutine oc_command_monitor(version,linkdate,narg,argline) 'UNIQUAC_MODEL ','DIFFUSION ','DEFAULT_CONSTIT ',& ' ','FCC_PERMUTATIONS','BCC_PERMUTATIONS',& ' ','GADDITION ','AQUEUS_MODEL ',& - 'QUASICHEM_MODEL ','FCC_CVM_TETRAHDR','FLORY_HUGG_MODEL',& + 'QUASICHEM_MODEL ','FCC_CVM_TETRAHDR',' ',& ' ',' ','QUIT '] !------------------- ! subsubsubcommands to PHASE ADDITION @@ -924,7 +924,7 @@ subroutine oc_command_monitor(version,linkdate,narg,argline) ! ' ',' ','DEFAULT_CONSTIT ',& ! ' ','FCC_PERMUTATIONS','BCC_PERMUTATIONS',& ! ' ','GADDITION ','AQUEUS_MODEL ',& -! 'QUASICHEM_MODEL ','FCC_CVM_TETRAHDR','FLORY_HUGG_MODEL',& +! 'QUASICHEM_MODEL ','FCC_CVM_TETRAHDR',' ',& ! ' ',' ','QUIT '] !.................................................... CASE DEFAULT @@ -999,6 +999,13 @@ subroutine oc_command_monitor(version,linkdate,narg,argline) !.................................................... case(4) ! amend phase addition twostate_liquid model call add_addrecord(lokph,' ',twostatemodel1) + call gparcdx('Is the addition calculated for one mole atoms? ',& + cline,last,1,ch1,'Y','?Add per formula unit') +! The CP model calculates a molar Gibbs energy, must be multiplied with +! the number of atoms in the phase. j2 set above to the addition type + if(ch1.eq.'Y' .or. ch1.eq.'y') then + call setpermolebit(lokph,twostatemodel1) + endif write(kou,667) 667 format('This addition require THET parameters for the',& ' Einstein T of the amorphous state'/'and G2 parameters',& @@ -1017,7 +1024,7 @@ subroutine oc_command_monitor(version,linkdate,narg,argline) case(7) ! amend phase LowT_CP_model call add_addrecord(lokph,' ',einsteincp) write(*,*)'This addition requires the THET parameter' - call gparcdx('Is the addition calculated for one mole? ',& + call gparcdx('Is the addition calculated for one mole atoms? ',& cline,last,1,ch1,'Y','?Add per formula unit') ! The CP model calculates a molar Gibbs energy, must be multiplied with ! the number of atoms in the phase. j2 set above to the addition type @@ -1038,13 +1045,17 @@ subroutine oc_command_monitor(version,linkdate,narg,argline) !.................................................... case(12) ! amend phase ... smooth-Cp-step call add_addrecord(lokph,' ',secondeinstein) - write(*,672) -672 format('This addition recures the THT2 and DCP2 parameters') -! The smooth CP model calculates a molar Gibbs energy, must be multiplied with + call gparcdx('Is the addition calculated for one mole? ',& + cline,last,1,ch1,'Y','?Add per formula unit') +! The CP model calculates a molar Gibbs energy, must be multiplied with ! the number of atoms in the phase. j2 set above to the addition type if(ch1.eq.'Y' .or. ch1.eq.'y') then call setpermolebit(lokph,secondeinstein) endif + write(*,672) +672 format('This addition recures the THT2 and DCP2 parameters') +! The smooth CP model calculates a molar Gibbs energy, must be multiplied with +! the number of atoms in the phase. j2 set above to the addition type end select amendphaseadd !************************************ end of amend phase ... addition !.................................................... @@ -1159,10 +1170,8 @@ subroutine oc_command_monitor(version,linkdate,narg,argline) write(kou,*)'Not implemented yet' ! call set_phase_status_bit(lokph,PHCVMCE) !.................................................... - case(15) ! amend phase ... FLORY_HUGGINS_MODEL + case(15) ! amend phase ... unused write(kou,*)'Not implemented, use UNIQUAC' -! call set_phase_status_bit(lokph,PHFHV) -! call clear_phase_status_bit(lokph,PHID) !.................................................... case(16) ! moved !....................................................