tuning Mood

Src2D2/main_dec.f90

@@ -11,6 +11,7 @@
 !!!  July 10, 2018
 !!!  Correspondance:	remi.abgrall@math.uzh.ch
 !!!  ------------------------------------------
+
 PROGRAM main
   USE param2d
   USE scheme
@@ -299,10 +300,12 @@ PROGRAM main
 
               ! do the test at the end of timestep, not at every iteration!!!!
               ! maybe positivity checks inside every iteration
-
+              ! or do it
               CALL test(DATA%iordret-1,Debug,Var, mesh, DATA,fluxes_mood(nflux))
 
            ENDDO
+           ! or do it only after all time cycles
+
            !              CALL test(DATA%iordret-1,Debug,Var, mesh, DATA,fluxes_mood(nflux))
 
            !print*

Src2D2/update.f90

@@ -55,7 +55,7 @@ MODULE update
 
 
   INTERFACE main_update_corrected
-!     MODULE PROCEDURE main_update_cor_entropy
+     !     MODULE PROCEDURE main_update_cor_entropy
      MODULE PROCEDURE main_update_cor_kinetic
   END INTERFACE main_update_corrected
 
@@ -137,19 +137,19 @@ CONTAINS
   END SUBROUTINE main_update
 
 
-   SUBROUTINE main_update_cor_kinetic(k,jt,dt, e,Var, DATA,alpha,beta,gamma,n_theta,theta)
+  SUBROUTINE main_update_cor_kinetic(k,jt,dt, e,Var, DATA,alpha,beta,gamma,n_theta,theta)
     IMPLICIT NONE
     CHARACTER(LEN = *), PARAMETER :: mod_name = "main_update"
     REAL(DP), INTENT(in):: dt
     INTEGER, INTENT(in):: k,jt
 
-        REAL(DP),DIMENSION(3,3:4), INTENT(in):: alpha
+    REAL(DP),DIMENSION(3,3:4), INTENT(in):: alpha
     REAL(DP),DIMENSION(3,3,3:4), INTENT(in):: beta, gamma
 
     INTEGER, DIMENSION(-4:4), INTENT(in):: n_theta
     REAL(DP),DIMENSION(0:3,3,-4:4), INTENT(in):: theta
 
-    
+
 
 
     TYPE(element), INTENT(in):: e
@@ -221,7 +221,7 @@ CONTAINS
 
   END SUBROUTINE main_update_cor_kinetic
 
-  
+
   SUBROUTINE main_update_cor_kinetic_old(k,jt,dt, e,Var, DATA,alpha,beta,gamma,n_theta,theta)
     IMPLICIT NONE
     CHARACTER(LEN = *), PARAMETER :: mod_name = "main_update"
@@ -411,7 +411,7 @@ CONTAINS
     TYPE(PVar), DIMENSION(:),ALLOCATABLE       :: u_old, v_old, u_new, v_new
     TYPE(Pvar), DIMENSION(:), ALLOCATABLE      :: phys, phys_d
     INTEGER, DIMENSION(Mesh%ndofs)             :: diag, diag2
-    TYPE(gradi), DIMENSION(:),ALLOCATABLE, SAVE:: gradient_d!, gradient
+    TYPE(gradi), DIMENSION(:),ALLOCATABLE, SAVE:: gradient_d, gradient
     REAL(dp), DIMENSION(:,:),ALLOCATABLE       :: grad
     REAL(dp), DIMENSION(:,:), ALLOCATABLE      :: points
     REAL(dp), DIMENSION(:,:), ALLOCATABLE      :: x
@@ -425,10 +425,10 @@ CONTAINS
     INTEGER, DIMENSION(1)                      :: nombre
 
     IF (.NOT.ALLOCATED(gradient_d)) THEN
-       ALLOCATE(gradient_d(Mesh%nt))!, gradient(Mesh%nt))
+       ALLOCATE(gradient_d(Mesh%nt), gradient(Mesh%nt))
        DO jt=1, Mesh%nt
           ALLOCATE(gradient_d(jt)%grad(n_dim,mesh%e(jt)%nvertex) )
-          !         ALLOCATE(gradient(  jt)%grad(n_dim,mesh%e(jt)%nvertex) )
+          ALLOCATE(gradient(  jt)%grad(n_dim,mesh%e(jt)%nvertex) )
        ENDDO
     ENDIF
 
@@ -521,13 +521,16 @@ CONTAINS
              !             DO l=1, e%nsommets
              !                grad(:,l)=e%grad_at_dofs(:,l,k)!gradient(l, e%coorL(:,k))
              !             ENDDO
-             gradient_d(jt)%grad(1,id)=SUM((phys_d(e%nu)%u(lk)&
-                  & -phys_d(e%nu(1))%u(lk))*e%grad_at_dofs(1,1,id) ) !grad(1,:))
-             gradient_d(jt)%grad(2,id)=SUM((phys_d(e%nu)%u(lk)-phys_d(e%nu(1))%u(lk))*&
+             gradient_d(jt)%grad(1,id)=SUM((phys_d(e%nu)%u(lk) - phys_d(e%nu(1))%u(lk) ) &
+                  &*e%grad_at_dofs(1,1,id) ) !grad(1,:))
+             gradient_d(jt)%grad(2,id)=SUM((phys_d(e%nu)%u(lk) - phys_d(e%nu(1))%u(lk) )*&
+                  & e%grad_at_dofs(2,1,id) ) !grad(2,:))
+
+             gradient(jt)%grad(1,id)=SUM((phys(e%nu)%u(lk)  - phys(e%nu(1))%u(lk) ) &
+                  & *e%grad_at_dofs(1,1,id) ) !grad(1,:))
+             gradient(jt)%grad(2,id)=SUM( (phys(e%nu)%u(lk) - phys(e%nu(1))%u(lk) )*&
                   & e%grad_at_dofs(2,1,id) ) !grad(2,:))
 
-             !             gradient(jt)%grad(1,k)=SUM((phys(e%nu)%u(lk)-phys(e%nu(1))%u(lk))*grad(1,:))
-             !             gradient(jt)%grad(2,k)=SUM((phys(e%nu)%u(lk)-phys(e%nu(1))%u(lk))*grad(2,:))
 
 
           ENDDO
@@ -560,12 +563,12 @@ CONTAINS
        ! here test if the min of density or pressure for the debug variables are >0
 
        !IF ((dvmin(is)%u(1).LT.0._dp) .OR. (dvmin(is)%u(n_vars).LT.0._dp )) THEN
-       IF ((phys_d(is)%u(1) .LT.0._dp) .OR. (phys_d(is)%u(n_vars) .LT.0._dp) ) then
+       IF ((phys_d(is)%u(1) .LT.0._dp) .OR. (phys_d(is)%u(n_vars) .LT.0._dp) ) THEN
           diag(is)=MAX(PAD_criteria,diag(is))
           diag2(is)=diag2(is)+diag(is)
        ENDIF
     ENDDO LoopPad
-  
+
 
     DO jt=1, Mesh%nt
        e=Mesh%e(jt)
@@ -607,7 +610,7 @@ CONTAINS
           IF (e%diag.GE.Plateau) CYCLE loopDMP
           ! Extrema detection: compare old and new variables
           !eps=coeff*MAX(( u_max%u(lk)-u_min%u(lk))*eps1, eps2)
-          eps=coeff*MAX(( umax(jt)%u(lk)-umin(jt)%u(lk))*eps1, eps2)
+          eps=-coeff*MAX(( umax(jt)%u(lk)-umin(jt)%u(lk))*eps1, eps2) !+
           IF (dumin(jt)%u(lk).LT.umin(jt)%u(lk)-eps.OR.dumax(jt)%u(lk).GT.umax(jt)%u(lk)+eps) THEN
 
 
@@ -624,7 +627,8 @@ CONTAINS
                    jtt=Mesh%vois(e%nu(l))%nvois(p)
                    IF (jtt.NE.jt) THEN
                       compte=compte+1
-                      points(:,compte)=gradient_d(jtt)%grad(:,Mesh%vois(e%nu(l))%loc(p) )
+                      !points(:,compte)=gradient_d(jtt)%grad(:,Mesh%vois(e%nu(l))%loc(p) )
+                      points(:,compte)=gradient(jtt)%grad(:,Mesh%vois(e%nu(l))%loc(p) )
                       !                      points(:,compte)=gradient(jtt)%grad(:,Mesh%vois(e%nu(l))%loc(p) )
                    ENDIF
                 ENDDO