@@ -80,6 +80,9 @@ subroutine dvnlsd(IWM, NFLAG, rwork, vstate)
8080 ! -----------------------------------------------------------------------
8181 !
8282 use vode_rhs_module, only: f_rhs, jac
83+ #ifdef CLEAN_INTEGRATOR_CORRECTION
84+ use vode_type_module, only: clean_state
85+ #endif
8386 use cuvode_dvnorm_module, only: dvnorm ! function
8487
8588 implicit none
@@ -192,12 +195,31 @@ subroutine dvnlsd(IWM, NFLAG, rwork, vstate)
192195 CSCALE = TWO/ (ONE + vstate % RC)
193196 CALL DSCALN (VODE_NEQS, CSCALE, vstate % Y, 1 )
194197 ENDIF
198+
199+ #ifdef CLEAN_INTEGRATOR_CORRECTION
200+ ! Clean the correction to Y. Use vstate % Y as scratch space.
201+
202+ ! Find the corrected Y: Yc = Y_previous + Y_delta
203+ do I = 1 ,VODE_NEQS
204+ vstate % Y(I) = vstate % Y(I) + (rwork % YH(I,1 ) + rwork % ACOR(I))
205+ end do
206+
207+ ! Clean the corrected Y: Yc' = clean(Yc)
208+ call clean_state(vstate % Y, vstate % RPAR)
209+
210+ ! Find the cleaned correction: clean(Y_delta) = Yc' - Y_previous
211+ do I = 1 ,VODE_NEQS
212+ vstate % Y(I) = vstate % Y(I) - (rwork % YH(I,1 ) + rwork % ACOR(I))
213+ end do
214+ #endif
215+
195216 DEL = DVNORM (vstate % Y, rwork % EWT)
196217 call daxpyn(VODE_NEQS, ONE, vstate % Y, 1 , rwork % acor, 1 )
197218
198219 do I = 1 ,VODE_NEQS
199220 vstate % Y(I) = rwork % YH(I,1 ) + rwork % ACOR(I)
200221 end do
222+
201223 ! -----------------------------------------------------------------------
202224 ! Test for convergence. If M .gt. 0, an estimate of the convergence
203225 ! rate constant is stored in CRATE, and this is used in the test.
0 commit comments