#include module shashimodule use petscsnes implicit none contains ! ! ------------------------------------------------------------------------ ! ! FormFunction - Evaluates nonlinear function, F(x). ! ! Input Parameters: ! snes - the SNES context ! x - input vector ! dummy - optional user-defined context (not used here) ! ! Output Parameter: ! f - function vector ! subroutine FormFunction(snes, x, f, dummy, ierr) SNES snes Vec x, f PetscErrorCode ierr integer dummy(*) ! Declarations for use with local arrays PetscScalar, pointer ::lx_v(:), lf_v(:) ! Get pointers to vector data. ! - For default PETSc vectors, VecGetArray() returns a pointer to ! the data array. Otherwise, the routine is implementation dependent. ! - You MUST call VecRestoreArray() when you no longer need access to ! the array. ! - Note that the Fortran interface to VecGetArray() differs from the ! C version. See the Fortran chapter of the users manual for details. PetscCall(VecGetArrayRead(x, lx_v, ierr)) PetscCall(VecGetArray(f, lf_v, ierr)) PetscCall(ShashiFormFunction(lx_v, lf_v)) PetscCall(VecRestoreArrayRead(x, lx_v, ierr)) PetscCall(VecRestoreArray(f, lf_v, ierr)) end ! --------------------------------------------------------------------- ! ! FormJacobian - Evaluates Jacobian matrix. ! ! Input Parameters: ! snes - the SNES context ! x - input vector ! dummy - optional user-defined context (not used here) ! ! Output Parameters: ! A - Jacobian matrix ! B - optionally different matrix used to construct the preconditioner ! subroutine FormJacobian(snes, X, jac, B, dummy, ierr) SNES snes Vec X Mat jac, B PetscErrorCode ierr integer dummy(*) ! Declarations for use with local arrays PetscScalar, pointer ::lx_v(:), lf_v(:, :) ! Get pointer to vector data PetscCall(VecGetArrayRead(x, lx_v, ierr)) PetscCall(MatDenseGetArray(B, lf_v, ierr)) PetscCall(ShashiFormJacobian(lx_v, lf_v)) PetscCall(MatDenseRestoreArray(B, lf_v, ierr)) PetscCall(VecRestoreArrayRead(x, lx_v, ierr)) ! Assemble matrix PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr)) PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr)) end subroutine ShashiLowerBound(an_r) PetscScalar an_r(26) PetscInt i do i = 2, 26 an_r(i) = 1000.0/6.023D+23 end do end subroutine ShashiInitialGuess(an_r) PetscScalar an_c_additive PetscScalar an_h_additive PetscScalar an_o_additive PetscScalar atom_c_init PetscScalar atom_h_init PetscScalar atom_n_init PetscScalar atom_o_init PetscScalar h_init PetscScalar p_init PetscInt nfuel PetscScalar temp, pt PetscScalar an_r(26) PetscInt an_h(1), an_c(1) pt = 0.1 atom_c_init = 6.7408177364816552D-022 atom_h_init = 2.0 atom_o_init = 1.0 atom_n_init = 3.76 nfuel = 1 an_c(1) = 1 an_h(1) = 4 an_c_additive = 2 an_h_additive = 6 an_o_additive = 1 h_init = 128799.7267952987 p_init = 0.1 temp = 1500 an_r(1) = 1.66000D-24 an_r(2) = 1.66030D-22 an_r(3) = 5.00000D-01 an_r(4) = 1.66030D-22 an_r(5) = 1.66030D-22 an_r(6) = 1.88000D+00 an_r(7) = 1.66030D-22 an_r(8) = 1.66030D-22 an_r(9) = 1.66030D-22 an_r(10) = 1.66030D-22 an_r(11) = 1.66030D-22 an_r(12) = 1.66030D-22 an_r(13) = 1.66030D-22 an_r(14) = 1.00000D+00 an_r(15) = 1.66030D-22 an_r(16) = 1.66030D-22 an_r(17) = 1.66000D-24 an_r(18) = 1.66030D-24 an_r(19) = 1.66030D-24 an_r(20) = 1.66030D-24 an_r(21) = 1.66030D-24 an_r(22) = 1.66030D-24 an_r(23) = 1.66030D-24 an_r(24) = 1.66030D-24 an_r(25) = 1.66030D-24 an_r(26) = 1.66030D-24 an_r = 0 an_r(3) = .5 an_r(6) = 1.88000 an_r(14) = 1. #if defined(solution) an_r(1) = 3.802208D-33 an_r(2) = 1.298287D-29 an_r(3) = 2.533067D-04 an_r(4) = 6.865078D-22 an_r(5) = 9.993125D-01 an_r(6) = 1.879964D+00 an_r(7) = 4.449489D-13 an_r(8) = 3.428687D-07 an_r(9) = 7.105138D-05 an_r(10) = 1.094368D-04 an_r(11) = 2.362305D-06 an_r(12) = 1.107145D-09 an_r(13) = 1.276162D-24 an_r(14) = 6.315538D-04 an_r(15) = 2.356540D-09 an_r(16) = 2.048248D-09 an_r(17) = 1.966187D-22 an_r(18) = 7.856497D-29 an_r(19) = 1.987840D-36 an_r(20) = 8.182441D-22 an_r(21) = 2.684880D-16 an_r(22) = 2.680473D-16 an_r(23) = 6.594967D-18 an_r(24) = 2.509714D-21 an_r(25) = 3.096459D-21 an_r(26) = 6.149551D-18 #endif end subroutine ShashiFormFunction(an_r, f_eq) PetscScalar an_c_additive PetscScalar an_h_additive PetscScalar an_o_additive PetscScalar atom_c_init PetscScalar atom_h_init PetscScalar atom_n_init PetscScalar atom_o_init PetscScalar h_init PetscScalar p_init PetscInt nfuel PetscScalar temp, pt PetscScalar an_r(26), k_eq(23), f_eq(26) PetscScalar H_molar(26) PetscInt an_h(1), an_c(1) PetscScalar part_p(26), idiff PetscInt i_cc, i_hh, i_h2o PetscScalar an_t, sum_h PetscScalar a_io2 PetscInt i, ip pt = 0.1 atom_c_init = 6.7408177364816552D-022 atom_h_init = 2.0 atom_o_init = 1.0 atom_n_init = 3.76 nfuel = 1 an_c(1) = 1 an_h(1) = 4 an_c_additive = 2 an_h_additive = 6 an_o_additive = 1 h_init = 128799.7267952987 p_init = 0.1 temp = 1500 k_eq(1) = 1.75149D-05 k_eq(2) = 4.01405D-06 k_eq(3) = 6.04663D-14 k_eq(4) = 2.73612D-01 k_eq(5) = 3.25592D-03 k_eq(6) = 5.33568D+05 k_eq(7) = 2.07479D+05 k_eq(8) = 1.11841D-02 k_eq(9) = 1.72684D-03 k_eq(10) = 1.98588D-07 k_eq(11) = 7.23600D+27 k_eq(12) = 5.73926D+49 k_eq(13) = 1.00000D+00 k_eq(14) = 1.64493D+16 k_eq(15) = 2.73837D-29 k_eq(16) = 3.27419D+50 k_eq(17) = 1.72447D-23 k_eq(18) = 4.24657D-06 k_eq(19) = 1.16065D-14 k_eq(20) = 3.28020D+25 k_eq(21) = 1.06291D+00 k_eq(22) = 9.11507D+02 k_eq(23) = 6.02837D+03 H_molar(1) = 3.26044D+03 H_molar(2) = -8.00407D+04 H_molar(3) = 4.05873D+04 H_molar(4) = -3.31849D+05 H_molar(5) = -1.93654D+05 H_molar(6) = 3.84035D+04 H_molar(7) = 4.97589D+05 H_molar(8) = 2.74483D+05 H_molar(9) = 1.30022D+05 H_molar(10) = 7.58429D+04 H_molar(11) = 2.42948D+05 H_molar(12) = 1.44588D+05 H_molar(13) = -7.16891D+04 H_molar(14) = 3.63075D+04 H_molar(15) = 9.23880D+04 H_molar(16) = 6.50477D+04 H_molar(17) = 3.04310D+05 H_molar(18) = 7.41707D+05 H_molar(19) = 6.32767D+05 H_molar(20) = 8.90624D+05 H_molar(21) = 2.49805D+04 H_molar(22) = 6.43473D+05 H_molar(23) = 1.02861D+06 H_molar(24) = -6.07503D+03 H_molar(25) = 1.27020D+05 H_molar(26) = -1.07011D+05 !============= an_t = 0 sum_h = 0 do i = 1, 26 an_t = an_t + an_r(i) end do f_eq(1) = atom_h_init & - (an_h(1)*an_r(1) + an_h_additive*an_r(2) & + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14) & + an_r(16) + 2*an_r(17) + an_r(19) & + an_r(20) + 3*an_r(22) + an_r(26)) f_eq(2) = atom_o_init & - (an_o_additive*an_r(2) + 2*an_r(3) & + 2*an_r(4) + an_r(5) & + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) & + 2*an_r(15) + 2*an_r(16) + an_r(20) + an_r(22) & + an_r(23) + 2*an_r(24) + 1*an_r(25) + an_r(26)) f_eq(3) = an_r(2) - 1.0d-150 f_eq(4) = atom_c_init & - (an_c(1)*an_r(1) + an_c_additive*an_r(2) & + an_r(4) + an_r(13) + 2*an_r(17) + an_r(18) & + an_r(19) + an_r(20)) do ip = 1, 26 part_p(ip) = (an_r(ip)/an_t)*pt end do i_cc = an_c(1) i_hh = an_h(1) a_io2 = i_cc + i_hh/4.0 i_h2o = i_hh/2 idiff = (i_cc + i_h2o) - (a_io2 + 1) f_eq(5) = k_eq(11)*an_r(1)*an_r(3)**a_io2 & - (an_r(4)**i_cc)*(an_r(5)**i_h2o)*((pt/an_t)**idiff) ! write(6,*)f_eq(5),an_r(1), an_r(3), an_r(4),an_r(5),' II' ! stop f_eq(6) = atom_n_init & - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12) & + an_r(15) & + an_r(23)) f_eq(7) = part_p(11) & - (k_eq(1)*sqrt(part_p(14) + 1d-23)) f_eq(8) = part_p(8) & - (k_eq(2)*sqrt(part_p(3) + 1d-23)) f_eq(9) = part_p(7) & - (k_eq(3)*sqrt(part_p(6) + 1d-23)) f_eq(10) = part_p(10) & - (k_eq(4)*sqrt(part_p(3) + 1d-23)) & *sqrt(part_p(14)) f_eq(11) = part_p(9) & - (k_eq(5)*sqrt(part_p(3) + 1d-23)) & *sqrt(part_p(6) + 1d-23) f_eq(12) = part_p(5) & - (k_eq(6)*sqrt(part_p(3) + 1d-23)) & *(part_p(14)) f_eq(13) = part_p(4) & - (k_eq(7)*sqrt(part_p(3) + 1.0d-23)) & *(part_p(13)) f_eq(14) = part_p(15) & - (k_eq(8)*sqrt(part_p(3) + 1.0d-50)) & *(part_p(9)) f_eq(15) = part_p(16) & - (k_eq(9)*part_p(3)) & *sqrt(part_p(14) + 1d-23) f_eq(16) = part_p(12) & - (k_eq(10)*sqrt(part_p(3) + 1d-23)) & *(part_p(6)) f_eq(17) = part_p(14)*part_p(18)**2 & - (k_eq(15)*part_p(17)) f_eq(18) = (part_p(13)**2) & - (k_eq(16)*part_p(3)*part_p(18)**2) print *, f_eq(18), part_p(3), part_p(18), part_p(13), k_eq(16) f_eq(19) = part_p(19)*part_p(3) - k_eq(17)*part_p(13)*part_p(10) f_eq(20) = part_p(21)*part_p(20) - k_eq(18)*part_p(19)*part_p(8) f_eq(21) = part_p(21)*part_p(23) - k_eq(19)*part_p(7)*part_p(8) f_eq(22) = part_p(5)*part_p(11) - k_eq(20)*part_p(21)*part_p(22) f_eq(23) = part_p(24) - k_eq(21)*part_p(21)*part_p(3) f_eq(24) = part_p(3)*part_p(25) - k_eq(22)*part_p(24)*part_p(8) f_eq(25) = part_p(26) - k_eq(23)*part_p(21)*part_p(10) f_eq(26) = -(an_r(20) + an_r(22) + an_r(23)) & + (an_r(21) + an_r(24) + an_r(25) + an_r(26)) do i = 1, 26 write (44, *) i, f_eq(i) end do end subroutine ShashiFormJacobian(an_r, d_eq) PetscScalar an_c_additive PetscScalar an_h_additive PetscScalar an_o_additive PetscScalar atom_c_init PetscScalar atom_h_init PetscScalar atom_n_init PetscScalar atom_o_init PetscScalar h_init PetscScalar p_init PetscInt nfuel PetscScalar temp, pt PetscScalar an_t, ai_o2 PetscScalar an_tot1_d, an_tot1 PetscScalar an_tot2_d, an_tot2 PetscScalar const5, const3, const2 PetscScalar const_cube PetscScalar const_five PetscScalar const_four PetscScalar const_six PetscInt jj, jb, ii3, id, ib, i PetscScalar pt2, pt1 PetscScalar an_r(26), k_eq(23) PetscScalar d_eq(26, 26), H_molar(26) PetscInt an_h(1), an_c(1) PetscScalar ai_pwr1, idiff PetscInt i_cc, i_hh PetscInt i_h2o, i_pwr2, i_co2_h2o PetscScalar pt_cube, pt_five PetscScalar pt_four PetscScalar pt_val1, pt_val2 PetscInt j pt = 0.1 atom_c_init = 6.7408177364816552D-022 atom_h_init = 2.0 atom_o_init = 1.0 atom_n_init = 3.76 nfuel = 1 an_c(1) = 1 an_h(1) = 4 an_c_additive = 2 an_h_additive = 6 an_o_additive = 1 h_init = 128799.7267952987 p_init = 0.1 temp = 1500 k_eq(1) = 1.75149D-05 k_eq(2) = 4.01405D-06 k_eq(3) = 6.04663D-14 k_eq(4) = 2.73612D-01 k_eq(5) = 3.25592D-03 k_eq(6) = 5.33568D+05 k_eq(7) = 2.07479D+05 k_eq(8) = 1.11841D-02 k_eq(9) = 1.72684D-03 k_eq(10) = 1.98588D-07 k_eq(11) = 7.23600D+27 k_eq(12) = 5.73926D+49 k_eq(13) = 1.00000D+00 k_eq(14) = 1.64493D+16 k_eq(15) = 2.73837D-29 k_eq(16) = 3.27419D+50 k_eq(17) = 1.72447D-23 k_eq(18) = 4.24657D-06 k_eq(19) = 1.16065D-14 k_eq(20) = 3.28020D+25 k_eq(21) = 1.06291D+00 k_eq(22) = 9.11507D+02 k_eq(23) = 6.02837D+03 H_molar(1) = 3.26044D+03 H_molar(2) = -8.00407D+04 H_molar(3) = 4.05873D+04 H_molar(4) = -3.31849D+05 H_molar(5) = -1.93654D+05 H_molar(6) = 3.84035D+04 H_molar(7) = 4.97589D+05 H_molar(8) = 2.74483D+05 H_molar(9) = 1.30022D+05 H_molar(10) = 7.58429D+04 H_molar(11) = 2.42948D+05 H_molar(12) = 1.44588D+05 H_molar(13) = -7.16891D+04 H_molar(14) = 3.63075D+04 H_molar(15) = 9.23880D+04 H_molar(16) = 6.50477D+04 H_molar(17) = 3.04310D+05 H_molar(18) = 7.41707D+05 H_molar(19) = 6.32767D+05 H_molar(20) = 8.90624D+05 H_molar(21) = 2.49805D+04 H_molar(22) = 6.43473D+05 H_molar(23) = 1.02861D+06 H_molar(24) = -6.07503D+03 H_molar(25) = 1.27020D+05 H_molar(26) = -1.07011D+05 ! !======= do jb = 1, 26 do ib = 1, 26 d_eq(ib, jb) = 0.0d0 end do end do an_t = 0.0 do id = 1, 26 an_t = an_t + an_r(id) end do !==== !==== d_eq(1, 1) = -an_h(1) d_eq(1, 2) = -an_h_additive d_eq(1, 5) = -2 d_eq(1, 10) = -1 d_eq(1, 11) = -1 d_eq(1, 14) = -2 d_eq(1, 16) = -1 d_eq(1, 17) = -2 d_eq(1, 19) = -1 d_eq(1, 20) = -1 d_eq(1, 22) = -3 d_eq(1, 26) = -1 d_eq(2, 2) = -1*an_o_additive d_eq(2, 3) = -2 d_eq(2, 4) = -2 d_eq(2, 5) = -1 d_eq(2, 8) = -1 d_eq(2, 9) = -1 d_eq(2, 10) = -1 d_eq(2, 12) = -1 d_eq(2, 13) = -1 d_eq(2, 15) = -2 d_eq(2, 16) = -2 d_eq(2, 20) = -1 d_eq(2, 22) = -1 d_eq(2, 23) = -1 d_eq(2, 24) = -2 d_eq(2, 25) = -1 d_eq(2, 26) = -1 d_eq(6, 6) = -2 d_eq(6, 7) = -1 d_eq(6, 9) = -1 d_eq(6, 12) = -2 d_eq(6, 15) = -1 d_eq(6, 23) = -1 d_eq(4, 1) = -an_c(1) d_eq(4, 2) = -an_c_additive d_eq(4, 4) = -1 d_eq(4, 13) = -1 d_eq(4, 17) = -2 d_eq(4, 18) = -1 d_eq(4, 19) = -1 d_eq(4, 20) = -1 !---------- const2 = an_t*an_t const3 = (an_t)*sqrt(an_t) const5 = an_t*const3 const_cube = an_t*an_t*an_t const_four = const2*const2 const_five = const2*const_cube const_six = const_cube*const_cube pt_cube = pt*pt*pt pt_four = pt_cube*pt pt_five = pt_cube*pt*pt i_cc = an_c(1) i_hh = an_h(1) ai_o2 = i_cc + float(i_hh)/4.0 i_co2_h2o = i_cc + i_hh/2 i_h2o = i_hh/2 ai_pwr1 = 1 + i_cc + float(i_hh)/4.0 i_pwr2 = i_cc + i_hh/2 idiff = (i_cc + i_h2o) - (ai_o2 + 1) pt1 = pt**(ai_pwr1) an_tot1 = an_t**(ai_pwr1) pt_val1 = (pt/an_t)**(ai_pwr1) an_tot1_d = an_tot1*an_t pt2 = pt**(i_pwr2) an_tot2 = an_t**(i_pwr2) pt_val2 = (pt/an_t)**(i_pwr2) an_tot2_d = an_tot2*an_t d_eq(5, 1) = & -(an_r(4)**i_cc)*(an_r(5)**i_h2o) & *((pt/an_t)**idiff)*(-idiff/an_t) do jj = 2, 26 d_eq(5, jj) = d_eq(5, 1) end do d_eq(5, 1) = d_eq(5, 1) + k_eq(11)*(an_r(3)**ai_o2) d_eq(5, 3) = d_eq(5, 3) + k_eq(11)*(ai_o2*an_r(3)**(ai_o2 - 1)) & & *an_r(1) d_eq(5, 4) = d_eq(5, 4) - (i_cc*an_r(4)**(i_cc - 1))* & (an_r(5)**i_h2o)*((pt/an_t)**idiff) d_eq(5, 5) = d_eq(5, 5) & - (i_h2o*(an_r(5)**(i_h2o - 1))) & *(an_r(4)**i_cc)*((pt/an_t)**idiff) d_eq(3, 1) = -(an_r(4)**2)*(an_r(5)**3)*(pt/an_t)*(-1.0/an_t) do jj = 2, 26 d_eq(3, jj) = d_eq(3, 1) end do d_eq(3, 2) = d_eq(3, 2) + k_eq(12)*(an_r(3)**3) d_eq(3, 3) = d_eq(3, 3) + k_eq(12)*(3*an_r(3)**2)*an_r(2) d_eq(3, 4) = d_eq(3, 4) - 2*an_r(4)*(an_r(5)**3)*(pt/an_t) d_eq(3, 5) = d_eq(3, 5) - 3*(an_r(5)**2)*(an_r(4)**2)*(pt/an_t) ! & *(pt_five/const_five) do ii3 = 1, 26 d_eq(3, ii3) = 0.0d0 end do d_eq(3, 2) = 1.0d0 d_eq(7, 1) = pt*an_r(11)*(-1.0)/const2 & - k_eq(1)*sqrt(pt)*sqrt(an_r(14) + 1d-50)*(-0.5/const3) do jj = 2, 26 d_eq(7, jj) = d_eq(7, 1) end do d_eq(7, 11) = d_eq(7, 11) + pt/an_t d_eq(7, 14) = d_eq(7, 14) & - k_eq(1)*sqrt(pt)*(0.5/(sqrt((an_r(14) + 1d-50)*an_t))) d_eq(8, 1) = pt*an_r(8)*(-1.0)/const2 & - k_eq(2)*sqrt(pt)*sqrt(an_r(3) + 1.0d-50)*(-0.5/const3) do jj = 2, 26 d_eq(8, jj) = d_eq(8, 1) end do d_eq(8, 3) = d_eq(8, 3) & - k_eq(2)*sqrt(pt)*(0.5/(sqrt((an_r(3) + 1.0d-50)*an_t))) d_eq(8, 8) = d_eq(8, 8) + pt/an_t d_eq(9, 1) = pt*an_r(7)*(-1.0)/const2 & - k_eq(3)*sqrt(pt)*sqrt(an_r(6))*(-0.5/const3) do jj = 2, 26 d_eq(9, jj) = d_eq(9, 1) end do d_eq(9, 7) = d_eq(9, 7) + pt/an_t d_eq(9, 6) = d_eq(9, 6) & - k_eq(3)*sqrt(pt)*(0.5/(sqrt(an_r(6)*an_t))) d_eq(10, 1) = pt*an_r(10)*(-1.0)/const2 & - k_eq(4)*(pt)*sqrt((an_r(3) + 1.0d-50) & *an_r(14))*(-1.0/const2) do jj = 2, 26 d_eq(10, jj) = d_eq(10, 1) end do d_eq(10, 3) = d_eq(10, 3) & - k_eq(4)*(pt)*sqrt(an_r(14)) & *(0.5/(sqrt(an_r(3) + 1.0d-50)*an_t)) d_eq(10, 10) = d_eq(10, 10) + pt/an_t d_eq(10, 14) = d_eq(10, 14) & - k_eq(4)*(pt)*sqrt(an_r(3) + 1.0d-50) & *(0.5/(sqrt(an_r(14) + 1.0d-50)*an_t)) d_eq(11, 1) = pt*an_r(9)*(-1.0)/const2 & - k_eq(5)*(pt)*sqrt((an_r(3) + 1.0d-50)*an_r(6)) & *(-1.0/const2) do jj = 2, 26 d_eq(11, jj) = d_eq(11, 1) end do d_eq(11, 3) = d_eq(11, 3) & - k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/ & (sqrt(an_r(3) + 1.0d-50)*an_t)) d_eq(11, 6) = d_eq(11, 6) & - k_eq(5)*(pt)*sqrt(an_r(3) + 1.0d-50) & *(0.5/(sqrt(an_r(6))*an_t)) d_eq(11, 9) = d_eq(11, 9) + pt/an_t d_eq(12, 1) = pt*an_r(5)*(-1.0)/const2 & - k_eq(6)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) & *(an_r(14))*(-1.5/const5) do jj = 2, 26 d_eq(12, jj) = d_eq(12, 1) end do d_eq(12, 3) = d_eq(12, 3) & - k_eq(6)*(pt**1.5)*((an_r(14) + 1.0d-50)/const3) & *(0.5/sqrt(an_r(3) + 1.0d-50)) d_eq(12, 5) = d_eq(12, 5) + pt/an_t d_eq(12, 14) = d_eq(12, 14) & - k_eq(6)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3) d_eq(13, 1) = pt*an_r(4)*(-1.0)/const2 & - k_eq(7)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) & *(an_r(13))*(-1.5/const5) do jj = 2, 26 d_eq(13, jj) = d_eq(13, 1) end do d_eq(13, 3) = d_eq(13, 3) & - k_eq(7)*(pt**1.5)*(an_r(13)/const3) & *(0.5/sqrt(an_r(3) + 1.0d-50)) d_eq(13, 4) = d_eq(13, 4) + pt/an_t d_eq(13, 13) = d_eq(13, 13) & - k_eq(7)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3) d_eq(14, 1) = pt*an_r(15)*(-1.0)/const2 & - k_eq(8)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) & *(an_r(9))*(-1.5/const5) do jj = 2, 26 d_eq(14, jj) = d_eq(14, 1) end do d_eq(14, 3) = d_eq(14, 3) & - k_eq(8)*(pt**1.5)*(an_r(9)/const3) & *(0.5/sqrt(an_r(3) + 1.0d-50)) d_eq(14, 9) = d_eq(14, 9) & - k_eq(8)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3) d_eq(14, 15) = d_eq(14, 15) + pt/an_t d_eq(15, 1) = pt*an_r(16)*(-1.0)/const2 & - k_eq(9)*(pt**1.5)*sqrt(an_r(14) + 1.0d-50) & *(an_r(3))*(-1.5/const5) do jj = 2, 26 d_eq(15, jj) = d_eq(15, 1) end do d_eq(15, 3) = d_eq(15, 3) & - k_eq(9)*(pt**1.5)*(sqrt(an_r(14) + 1.0d-50)/const3) d_eq(15, 14) = d_eq(15, 14) & - k_eq(9)*(pt**1.5)*(an_r(3)/const3) & *(0.5/sqrt(an_r(14) + 1.0d-50)) d_eq(15, 16) = d_eq(15, 16) + pt/an_t d_eq(16, 1) = pt*an_r(12)*(-1.0)/const2 & - k_eq(10)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) & *(an_r(6))*(-1.5/const5) do jj = 2, 26 d_eq(16, jj) = d_eq(16, 1) end do d_eq(16, 3) = d_eq(16, 3) & - k_eq(10)*(pt**1.5)*(an_r(6)/const3) & *(0.5/sqrt(an_r(3) + 1.0d-50)) d_eq(16, 6) = d_eq(16, 6) & - k_eq(10)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3) d_eq(16, 12) = d_eq(16, 12) + pt/an_t const_cube = an_t*an_t*an_t const_four = const2*const2 d_eq(17, 1) = an_r(14)*an_r(18)*an_r(18)*(pt**3)*(-3/const_four) & - k_eq(15)*an_r(17)*pt*(-1/const2) do jj = 2, 26 d_eq(17, jj) = d_eq(17, 1) end do d_eq(17, 14) = d_eq(17, 14) + an_r(18)*an_r(18)*(pt**3)/const_cube d_eq(17, 17) = d_eq(17, 17) - k_eq(15)*pt/an_t d_eq(17, 18) = d_eq(17, 18) + 2*an_r(18)*an_r(14) & & *(pt**3)/const_cube d_eq(18, 1) = an_r(13)*an_r(13)*(pt**2)*(-2/const_cube) & - k_eq(16)*an_r(3)*an_r(18)*an_r(18) & *(pt*pt*pt)*(-3/const_four) do jj = 2, 26 d_eq(18, jj) = d_eq(18, 1) end do d_eq(18, 3) = d_eq(18, 3) & - k_eq(16)*an_r(18)*an_r(18)*pt*pt*pt/const_cube d_eq(18, 13) = d_eq(18, 13) & + 2*an_r(13)*pt*pt/const2 d_eq(18, 18) = d_eq(18, 18) - k_eq(16)*an_r(3) & & *2*an_r(18)*pt*pt*pt/const_cube !====for eq 19 d_eq(19, 1) = an_r(3)*an_r(19)*(pt**2)*(-2/const_cube) & - k_eq(17)*an_r(13)*an_r(10)*pt*pt*(-2/const_cube) do jj = 2, 26 d_eq(19, jj) = d_eq(19, 1) end do d_eq(19, 13) = d_eq(19, 13) & - k_eq(17)*an_r(10)*pt*pt/const2 d_eq(19, 10) = d_eq(19, 10) & - k_eq(17)*an_r(13)*pt*pt/const2 d_eq(19, 3) = d_eq(19, 3) + an_r(19)*pt*pt/const2 d_eq(19, 19) = d_eq(19, 19) + an_r(3)*pt*pt/const2 !====for eq 20 d_eq(20, 1) = an_r(21)*an_r(20)*(pt**2)*(-2/const_cube) & - k_eq(18)*an_r(19)*an_r(8)*pt*pt*(-2/const_cube) do jj = 2, 26 d_eq(20, jj) = d_eq(20, 1) end do d_eq(20, 8) = d_eq(20, 8) & - k_eq(18)*an_r(19)*pt*pt/const2 d_eq(20, 19) = d_eq(20, 19) & - k_eq(18)*an_r(8)*pt*pt/const2 d_eq(20, 20) = d_eq(20, 20) + an_r(21)*pt*pt/const2 d_eq(20, 21) = d_eq(20, 21) + an_r(20)*pt*pt/const2 !======== !====for eq 21 d_eq(21, 1) = an_r(21)*an_r(23)*(pt**2)*(-2/const_cube) & - k_eq(19)*an_r(7)*an_r(8)*pt*pt*(-2/const_cube) do jj = 2, 26 d_eq(21, jj) = d_eq(21, 1) end do d_eq(21, 7) = d_eq(21, 7) & - k_eq(19)*an_r(8)*pt*pt/const2 d_eq(21, 8) = d_eq(21, 8) & - k_eq(19)*an_r(7)*pt*pt/const2 d_eq(21, 21) = d_eq(21, 21) + an_r(23)*pt*pt/const2 d_eq(21, 23) = d_eq(21, 23) + an_r(21)*pt*pt/const2 !======== ! for 22 d_eq(22, 1) = an_r(5)*an_r(11)*(pt**2)*(-2/const_cube) & - k_eq(20)*an_r(21)*an_r(22)*pt*pt*(-2/const_cube) do jj = 2, 26 d_eq(22, jj) = d_eq(22, 1) end do d_eq(22, 21) = d_eq(22, 21) & - k_eq(20)*an_r(22)*pt*pt/const2 d_eq(22, 22) = d_eq(22, 22) & - k_eq(20)*an_r(21)*pt*pt/const2 d_eq(22, 11) = d_eq(22, 11) + an_r(5)*pt*pt/(const2) d_eq(22, 5) = d_eq(22, 5) + an_r(11)*pt*pt/(const2) !======== ! for 23 d_eq(23, 1) = an_r(24)*(pt)*(-1/const2) & - k_eq(21)*an_r(21)*an_r(3)*pt*pt*(-2/const_cube) do jj = 2, 26 d_eq(23, jj) = d_eq(23, 1) end do d_eq(23, 3) = d_eq(23, 3) & - k_eq(21)*an_r(21)*pt*pt/const2 d_eq(23, 21) = d_eq(23, 21) & - k_eq(21)*an_r(3)*pt*pt/const2 d_eq(23, 24) = d_eq(23, 24) + pt/(an_t) !======== ! for 24 d_eq(24, 1) = an_r(3)*an_r(25)*(pt**2)*(-2/const_cube) & - k_eq(22)*an_r(24)*an_r(8)*pt*pt*(-2/const_cube) do jj = 2, 26 d_eq(24, jj) = d_eq(24, 1) end do d_eq(24, 8) = d_eq(24, 8) & - k_eq(22)*an_r(24)*pt*pt/const2 d_eq(24, 24) = d_eq(24, 24) & - k_eq(22)*an_r(8)*pt*pt/const2 d_eq(24, 3) = d_eq(24, 3) + an_r(25)*pt*pt/const2 d_eq(24, 25) = d_eq(24, 25) + an_r(3)*pt*pt/const2 !======== !for 25 d_eq(25, 1) = an_r(26)*(pt)*(-1/const2) & - k_eq(23)*an_r(21)*an_r(10)*pt*pt*(-2/const_cube) do jj = 2, 26 d_eq(25, jj) = d_eq(25, 1) end do d_eq(25, 10) = d_eq(25, 10) & - k_eq(23)*an_r(21)*pt*pt/const2 d_eq(25, 21) = d_eq(25, 21) & - k_eq(23)*an_r(10)*pt*pt/const2 d_eq(25, 26) = d_eq(25, 26) + pt/(an_t) !============ ! for 26 d_eq(26, 20) = -1 d_eq(26, 22) = -1 d_eq(26, 23) = -1 d_eq(26, 21) = 1 d_eq(26, 24) = 1 d_eq(26, 25) = 1 d_eq(26, 26) = 1 do j = 1, 26 do i = 1, 26 write (44, *) i, j, d_eq(i, j) end do end do end subroutine ShashiPostCheck(ls, X, Y, W, c_Y, c_W, dummy) SNESLineSearch ls PetscErrorCode ierr Vec X, Y, W PetscObject dummy PetscBool c_Y, c_W PetscScalar, pointer :: xx(:) PetscInt i PetscCall(VecGetArray(W, xx, ierr)) do i = 1, 26 if (xx(i) < 0.0) then xx(i) = 0.0 c_W = PETSC_TRUE end if if (xx(i) > 3.0) then xx(i) = 3.0 end if end do PetscCall(VecRestoreArray(W, xx, ierr)) end end module shashimodule program main use petsc use shashimodule implicit none ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Variable declarations ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! Variables: ! snes - nonlinear solver ! x, r - solution, residual vectors ! J - Jacobian matrix ! SNES snes Vec x, r, lb, ub Mat J PetscErrorCode ierr PetscInt i2 PetscMPIInt size PetscScalar, pointer :: xx(:) PetscScalar zero, big SNESLineSearch ls ! ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Beginning of program ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - PetscCallA(PetscInitialize(ierr)) PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)) PetscCheckA(size == 1, PETSC_COMM_WORLD, 1, 'requires one process') big = 2.88 big = PETSC_INFINITY zero = 0.0 i2 = 26 ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - - ! Create nonlinear solver context ! - - - - - - - - - -- - - - - - - - - - - - - - - - - - - - - - - - - - PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Create matrix and vector data structures; set corresponding routines ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Create vectors for solution and nonlinear function PetscCallA(VecCreateSeq(PETSC_COMM_SELF, i2, x, ierr)) PetscCallA(VecDuplicate(x, r, ierr)) ! Create Jacobian matrix data structure PetscCallA(MatCreateDense(PETSC_COMM_SELF, 26, 26, 26, 26, PETSC_NULL_SCALAR, J, ierr)) ! Set function evaluation routine and vector PetscCallA(SNESSetFunction(snes, r, FormFunction, 0, ierr)) ! Set Jacobian matrix data structure and Jacobian evaluation routine PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, 0, ierr)) PetscCallA(VecDuplicate(x, lb, ierr)) PetscCallA(VecDuplicate(x, ub, ierr)) PetscCallA(VecSet(lb, zero, ierr)) ! PetscCallA(VecGetArray(lb,xx,ierr)) ! PetscCallA(ShashiLowerBound(xx)) ! PetscCallA(VecRestoreArray(lb,xx,ierr)) PetscCallA(VecSet(ub, big, ierr)) ! PetscCallA(SNESVISetVariableBounds(snes,lb,ub,ierr)) PetscCallA(SNESGetLineSearch(snes, ls, ierr)) PetscCallA(SNESLineSearchSetPostCheck(ls, ShashiPostCheck, 0, ierr)) PetscCallA(SNESSetType(snes, SNESVINEWTONRSLS, ierr)) PetscCallA(SNESSetFromOptions(snes, ierr)) ! set initial guess PetscCallA(VecGetArray(x, xx, ierr)) PetscCallA(ShashiInitialGuess(xx)) PetscCallA(VecRestoreArray(x, xx, ierr)) PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Free work space. All PETSc objects should be destroyed when they ! are no longer needed. ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - PetscCallA(VecDestroy(lb, ierr)) PetscCallA(VecDestroy(ub, ierr)) PetscCallA(VecDestroy(x, ierr)) PetscCallA(VecDestroy(r, ierr)) PetscCallA(MatDestroy(J, ierr)) PetscCallA(SNESDestroy(snes, ierr)) PetscCallA(PetscFinalize(ierr)) end