1! 2! Description: Solve Ax=b. A comes from an anisotropic 2D thermal problem with Q1 FEM on domain (-1,1)^2. 3! Material conductivity given by tensor: 4! 5! D = | 1 0 | 6! | 0 epsilon | 7! 8! rotated by angle 'theta' (-theta <90> in degrees) with anisotropic parameter 'epsilon' (-epsilon <0.0>). 9! Blob right-hand side centered at C (-blob_center C(1),C(2) <0,0>) 10! Dirichlet BCs on y=-1 face. 11! 12! -out_matlab will generate binary files for A,x,b and a ex54f.m file that reads them and plots them in matlab. 13! 14! User can change anisotropic shape with function ex54_psi(). Negative theta will switch to a circular anisotropy. 15! 16 17! ----------------------------------------------------------------------- 18#include <petsc/finclude/petscksp.h> 19program main 20 use petscksp 21 implicit none 22 23 Vec xvec, bvec, uvec 24 Mat Amat 25 KSP ksp 26 PetscErrorCode ierr 27 PetscViewer viewer 28 PetscInt qj, qi, ne, M, Istart, Iend, geq 29 PetscInt ki, kj, nel, ll, j1, i1, ndf, f4 30 PetscInt f2, f9, f6, one 31 PetscInt :: idx(4) 32 PetscBool flg, out_matlab 33 PetscMPIInt size, rank 34 PetscScalar::ss(4, 4), val 35 PetscReal::shp(3, 9), sg(3, 9) 36 PetscReal::thk, a1, a2 37 PetscReal::theta, eps, h, x, y, xsj 38 PetscReal::coord(2, 4), dd(2, 2), ev(3), blb(2) 39 40 common/ex54_theta/theta 41! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 42! Beginning of program 43! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 44 PetscCallA(PetscInitialize(ierr)) 45 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)) 46 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr)) 47 one = 1 48! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 49! set parameters 50! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 51 f4 = 4 52 f2 = 2 53 f9 = 9 54 f6 = 6 55 ne = 9 56 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-ne', ne, flg, ierr)) 57 h = 2.0/real(ne) 58 M = (ne + 1)*(ne + 1) 59 theta = 90.0 60! theta is input in degrees 61 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-theta', theta, flg, ierr)) 62 theta = theta/57.2957795 63 eps = 1.0 64 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-epsilon', eps, flg, ierr)) 65 ki = 2 66 PetscCallA(PetscOptionsGetRealArray(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-blob_center', blb, ki, flg, ierr)) 67 if (.not. flg) then 68 blb(1) = 0.0 69 blb(2) = 0.0 70 else if (ki /= 2) then 71 print *, 'error: ', ki, ' arguments read for -blob_center. Needs to be two.' 72 end if 73 PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-out_matlab', out_matlab, flg, ierr)) 74 if (.not. flg) out_matlab = PETSC_FALSE 75 76 ev(1) = 1.0 77 ev(2) = eps*ev(1) 78! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 79! Compute the matrix and right-hand-side vector that define 80! the linear system, Ax = b. 81! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 82! Create matrix. When using MatCreate(), the matrix format can 83! be specified at runtime. 84 PetscCallA(MatCreate(PETSC_COMM_WORLD, Amat, ierr)) 85 PetscCallA(MatSetSizes(Amat, PETSC_DECIDE, PETSC_DECIDE, M, M, ierr)) 86 PetscCallA(MatSetType(Amat, MATAIJ, ierr)) 87 PetscCallA(MatSetOption(Amat, MAT_SPD, PETSC_TRUE, ierr)) 88 PetscCallA(MatSetOption(Amat, MAT_SPD_ETERNAL, PETSC_TRUE, ierr)) 89 if (size == 1) then 90 PetscCallA(MatSetType(Amat, MATAIJ, ierr)) 91 else 92 PetscCallA(MatSetType(Amat, MATMPIAIJ, ierr)) 93 end if 94 PetscCallA(MatMPIAIJSetPreallocation(Amat, f9, PETSC_NULL_INTEGER_ARRAY, f6, PETSC_NULL_INTEGER_ARRAY, ierr)) 95 PetscCallA(MatSetFromOptions(Amat, ierr)) 96 PetscCallA(MatSetUp(Amat, ierr)) 97 PetscCallA(MatGetOwnershipRange(Amat, Istart, Iend, ierr)) 98! Create vectors. Note that we form 1 vector from scratch and 99! then duplicate as needed. 100 PetscCallA(MatCreateVecs(Amat, PETSC_NULL_VEC, xvec, ierr)) 101 PetscCallA(VecSetFromOptions(xvec, ierr)) 102 PetscCallA(VecDuplicate(xvec, bvec, ierr)) 103 PetscCallA(VecDuplicate(xvec, uvec, ierr)) 104! Assemble matrix. 105! - Note that MatSetValues() uses 0-based row and column numbers 106! in Fortran as well as in C (as set here in the array "col"). 107 thk = 1.0 ! thickness 108 nel = 4 ! nodes per element (quad) 109 ndf = 1 110 call int2d(f2, sg) 111 do geq = Istart, Iend - 1, 1 112 qj = geq/(ne + 1); qi = mod(geq, (ne + 1)) 113 x = h*qi - 1.0; y = h*qj - 1.0 ! lower left corner (-1,-1) 114 if (qi < ne .and. qj < ne) then 115 coord(1, 1) = x; coord(2, 1) = y 116 coord(1, 2) = x + h; coord(2, 2) = y 117 coord(1, 3) = x + h; coord(2, 3) = y + h 118 coord(1, 4) = x; coord(2, 4) = y + h 119! form stiff 120 ss = 0.0 121 do ll = 1, 4 122 call shp2dquad(sg(1, ll), sg(2, ll), coord, shp, xsj, f2) 123 xsj = xsj*sg(3, ll)*thk 124 call thfx2d(ev, coord, shp, dd, f2, f2, f4, ex54_psi) 125 j1 = 1 126 do kj = 1, nel 127 a1 = (dd(1, 1)*shp(1, kj) + dd(1, 2)*shp(2, kj))*xsj 128 a2 = (dd(2, 1)*shp(1, kj) + dd(2, 2)*shp(2, kj))*xsj 129! Compute residual 130! p(j1) = p(j1) - a1*gradt(1) - a2*gradt(2) 131! Compute tangent 132 i1 = 1 133 do ki = 1, nel 134 ss(i1, j1) = ss(i1, j1) + a1*shp(1, ki) + a2*shp(2, ki) 135 i1 = i1 + ndf 136 end do 137 j1 = j1 + ndf 138 end do 139 end do 140 141 idx(1) = geq; idx(2) = geq + 1; idx(3) = geq + (ne + 1) + 1 142 idx(4) = geq + (ne + 1) 143 if (qj > 0) then 144 PetscCallA(MatSetValues(Amat, f4, idx, f4, idx, reshape(ss, [f4*f4]), ADD_VALUES, ierr)) 145 else ! a BC 146 do ki = 1, 4, 1 147 do kj = 1, 4, 1 148 if (ki < 3 .or. kj < 3) then 149 if (ki == kj) then 150 ss(ki, kj) = .1*ss(ki, kj) 151 else 152 ss(ki, kj) = 0.0 153 end if 154 end if 155 end do 156 end do 157 PetscCallA(MatSetValues(Amat, f4, idx, f4, idx, reshape(ss, [f4*f4]), ADD_VALUES, ierr)) 158 end if ! BC 159 end if ! add element 160 if (qj > 0) then ! set rhs 161 val = h*h*exp(-100*((x + h/2) - blb(1))**2)*exp(-100*((y + h/2) - blb(2))**2) 162 PetscCallA(VecSetValues(bvec, one, [geq], [val], INSERT_VALUES, ierr)) 163 end if 164 end do 165 PetscCallA(MatAssemblyBegin(Amat, MAT_FINAL_ASSEMBLY, ierr)) 166 PetscCallA(MatAssemblyEnd(Amat, MAT_FINAL_ASSEMBLY, ierr)) 167 PetscCallA(VecAssemblyBegin(bvec, ierr)) 168 PetscCallA(VecAssemblyEnd(bvec, ierr)) 169 170! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 171! Create the linear solver and set various options 172! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 173 174! Create linear solver context 175 176 PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr)) 177 178! Set operators. Here the matrix that defines the linear system 179! also serves as the matrix from which the preconditioner is constructed. 180 181 PetscCallA(KSPSetOperators(ksp, Amat, Amat, ierr)) 182 183! Set runtime options, e.g., 184! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol> 185! These options will override those specified above as long as 186! KSPSetFromOptions() is called _after_ any other customization 187! routines. 188 189 PetscCallA(KSPSetFromOptions(ksp, ierr)) 190 191! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 192! Solve the linear system 193! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 194 195 PetscCallA(KSPSolve(ksp, bvec, xvec, ierr)) 196 197! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 198! output 199! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 200 if (out_matlab) then 201 PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'Amat', FILE_MODE_WRITE, viewer, ierr)) 202 PetscCallA(MatView(Amat, viewer, ierr)) 203 PetscCallA(PetscViewerDestroy(viewer, ierr)) 204 205 PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'Bvec', FILE_MODE_WRITE, viewer, ierr)) 206 PetscCallA(VecView(bvec, viewer, ierr)) 207 PetscCallA(PetscViewerDestroy(viewer, ierr)) 208 209 PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'Xvec', FILE_MODE_WRITE, viewer, ierr)) 210 PetscCallA(VecView(xvec, viewer, ierr)) 211 PetscCallA(PetscViewerDestroy(viewer, ierr)) 212 213 PetscCallA(MatMult(Amat, xvec, uvec, ierr)) 214 val = -1.0 215 PetscCallA(VecAXPY(uvec, val, bvec, ierr)) 216 PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, 'Rvec', FILE_MODE_WRITE, viewer, ierr)) 217 PetscCallA(VecView(uvec, viewer, ierr)) 218 PetscCallA(PetscViewerDestroy(viewer, ierr)) 219 220 if (rank == 0) then 221 open (1, file='ex54f.m', FORM='formatted') 222 write (1, *) 'A = PetscBinaryRead(''Amat'');' 223 write (1, *) '[m n] = size(A);' 224 write (1, *) 'mm = sqrt(m);' 225 write (1, *) 'b = PetscBinaryRead(''Bvec'');' 226 write (1, *) 'x = PetscBinaryRead(''Xvec'');' 227 write (1, *) 'r = PetscBinaryRead(''Rvec'');' 228 write (1, *) 'bb = reshape(b,mm,mm);' 229 write (1, *) 'xx = reshape(x,mm,mm);' 230 write (1, *) 'rr = reshape(r,mm,mm)' 231! write (1,*) 'imagesc(bb')' 232! write (1,*) 'title('RHS'),' 233 write (1, *) 'figure,' 234 write (1, *) 'imagesc(xx'')' 235 write (1, 2002) eps, theta*57.2957795 236 write (1, *) 'figure,' 237 write (1, *) 'imagesc(rr'')' 238 write (1, *) 'title(''Residual''),' 239 close (1) 240 end if 241 end if 2422002 format('title(''Solution: esp='',d9.3,'', theta='',g8.3,''),') 243! Free work space. All PETSc objects should be destroyed when they 244! are no longer needed. 245 246 PetscCallA(VecDestroy(xvec, ierr)) 247 PetscCallA(VecDestroy(bvec, ierr)) 248 PetscCallA(VecDestroy(uvec, ierr)) 249 PetscCallA(MatDestroy(Amat, ierr)) 250 PetscCallA(KSPDestroy(ksp, ierr)) 251 PetscCallA(PetscFinalize(ierr)) 252 253contains 254! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 255! thfx2d - compute material tensor 256! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 257! Compute thermal gradient and flux 258 259 subroutine thfx2d(ev, xl, shp, dd, ndm, ndf, nel, dir) 260 261 PetscInt ndm, ndf, nel, i 262 PetscReal ev(2), xl(ndm, nel), shp(3, *), dir 263 PetscReal xx, yy, psi, cs, sn, c2, s2, dd(2, 2) 264 265 xx = 0.0 266 yy = 0.0 267 do i = 1, nel 268 xx = xx + shp(3, i)*xl(1, i) 269 yy = yy + shp(3, i)*xl(2, i) 270 end do 271 psi = dir(xx, yy) 272! Compute thermal flux 273 cs = cos(psi) 274 sn = sin(psi) 275 c2 = cs*cs 276 s2 = sn*sn 277 cs = cs*sn 278 279 dd(1, 1) = c2*ev(1) + s2*ev(2) 280 dd(2, 2) = s2*ev(1) + c2*ev(2) 281 dd(1, 2) = cs*(ev(1) - ev(2)) 282 dd(2, 1) = dd(1, 2) 283 284! flux(1) = -dd(1,1)*gradt(1) - dd(1,2)*gradt(2) 285! flux(2) = -dd(2,1)*gradt(1) - dd(2,2)*gradt(2) 286 287 end 288 289! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 290! shp2dquad - shape functions - compute derivatives w/r natural coords. 291! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 292 subroutine shp2dquad(s, t, xl, shp, xsj, ndm) 293!-----[--.----+----.----+----.-----------------------------------------] 294! Purpose: Shape function routine for 4-node isoparametric quads 295! 296! Inputs: 297! s,t - Natural coordinates of point 298! xl(ndm,*) - Nodal coordinates for element 299! ndm - Spatial dimension of mesh 300 301! Outputs: 302! shp(3,*) - Shape functions and derivatives at point 303! shp(1,i) = dN_i/dx or dN_i/dxi_1 304! shp(2,i) = dN_i/dy or dN_i/dxi_2 305! shp(3,i) = N_i 306! xsj - Jacobian determinant at point 307!-----[--.----+----.----+----.-----------------------------------------] 308 PetscInt ndm 309 PetscReal xo, xs, xt, yo, ys, yt, xsm, xsp, xtm 310 PetscReal xtp, ysm, ysp, ytm, ytp 311 PetscReal s, t, xsj, xsj1, sh, th, sp, tp, sm 312 PetscReal tm, xl(ndm, 4), shp(3, 4) 313 314! Set up interpolations 315 316 sh = 0.5*s 317 th = 0.5*t 318 sp = 0.5 + sh 319 tp = 0.5 + th 320 sm = 0.5 - sh 321 tm = 0.5 - th 322 shp(3, 1) = sm*tm 323 shp(3, 2) = sp*tm 324 shp(3, 3) = sp*tp 325 shp(3, 4) = sm*tp 326 327! Set up natural coordinate functions (times 4) 328 329 xo = xl(1, 1) - xl(1, 2) + xl(1, 3) - xl(1, 4) 330 xs = -xl(1, 1) + xl(1, 2) + xl(1, 3) - xl(1, 4) + xo*t 331 xt = -xl(1, 1) - xl(1, 2) + xl(1, 3) + xl(1, 4) + xo*s 332 yo = xl(2, 1) - xl(2, 2) + xl(2, 3) - xl(2, 4) 333 ys = -xl(2, 1) + xl(2, 2) + xl(2, 3) - xl(2, 4) + yo*t 334 yt = -xl(2, 1) - xl(2, 2) + xl(2, 3) + xl(2, 4) + yo*s 335 336! Compute jacobian (times 16) 337 338 xsj1 = xs*yt - xt*ys 339 340! Divide jacobian by 16 (multiply by .0625) 341 342 xsj = 0.0625*xsj1 343 if (xsj1 == 0.0) then 344 xsj1 = 1.0 345 else 346 xsj1 = 1.0/xsj1 347 end if 348 349! Divide functions by jacobian 350 351 xs = (xs + xs)*xsj1 352 xt = (xt + xt)*xsj1 353 ys = (ys + ys)*xsj1 354 yt = (yt + yt)*xsj1 355 356! Multiply by interpolations 357 358 ytm = yt*tm 359 ysm = ys*sm 360 ytp = yt*tp 361 ysp = ys*sp 362 xtm = xt*tm 363 xsm = xs*sm 364 xtp = xt*tp 365 xsp = xs*sp 366 367! Compute shape functions 368 369 shp(1, 1) = -ytm + ysm 370 shp(1, 2) = ytm + ysp 371 shp(1, 3) = ytp - ysp 372 shp(1, 4) = -ytp - ysm 373 shp(2, 1) = xtm - xsm 374 shp(2, 2) = -xtm - xsp 375 shp(2, 3) = -xtp + xsp 376 shp(2, 4) = xtp + xsm 377 378 end 379 380! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 381! int2d 382! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 383 subroutine int2d(l, sg) 384!-----[--.----+----.----+----.-----------------------------------------] 385! Purpose: Form Gauss points and weights for two dimensions 386 387! Inputs: 388! l - Number of points/direction 389 390! Outputs: 391! sg(3,*) - Array of points and weights 392!-----[--.----+----.----+----.-----------------------------------------] 393 PetscInt l, i, lr(9), lz(9) 394 PetscReal g, third, sg(3, *) 395 data lr/-1, 1, 1, -1, 0, 1, 0, -1, 0/, lz/-1, -1, 1, 1, -1, 0, 1, 0, 0/ 396 data third/0.3333333333333333/ 397 398! 2x2 integration 399 g = sqrt(third) 400 do i = 1, 4 401 sg(1, i) = g*lr(i) 402 sg(2, i) = g*lz(i) 403 sg(3, i) = 1.0 404 end do 405 406 end 407 408! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 409! ex54_psi - anisotropic material direction 410! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 411 PetscReal function ex54_psi(x, y) 412 PetscReal x, y, theta 413 common/ex54_theta/theta 414 ex54_psi = theta 415 if (theta < 0.) then ! circular 416 if (y == 0) then 417 ex54_psi = 2.0*atan(1.0) 418 else 419 ex54_psi = atan(-x/y) 420 end if 421 end if 422 end 423end 424 425! 426!/*TEST 427! 428! testset: 429! nsize: 4 430! args: -ne 39 -theta 30.0 -epsilon 1.e-1 -blob_center 0.,0. -ksp_type cg -pc_type gamg -pc_gamg_type agg -ksp_rtol 1e-4 -pc_gamg_aggressive_square_graph false -ksp_norm_type unpreconditioned 431! requires: !single 432! test: 433! suffix: misk 434! args: -pc_gamg_mat_coarsen_type misk -pc_gamg_aggressive_coarsening 0 -ksp_monitor_short 435! test: 436! suffix: mis 437! args: -pc_gamg_mat_coarsen_type mis -ksp_monitor_short 438! test: 439! suffix: hem 440! args: -pc_gamg_mat_coarsen_type hem -ksp_converged_reason 441! filter: sed -e "s/Linear solve converged due to CONVERGED_RTOL iterations 1[2-3]/Linear solve converged due to CONVERGED_RTOL iterations 11/g" 442! 443!TEST*/ 444