From 42470eb419e10cd0ebf6f60906a20a90f78bf147 Mon Sep 17 00:00:00 2001 From: wme7 Date: Mon, 29 Jul 2013 17:15:06 +0800 Subject: [PATCH] after fortran class --- myf90/ESBGK/Stu_DG2D/Makefile | 53 ++ myf90/ESBGK/Stu_DG2D/README.md | 4 + myf90/ESBGK/Stu_DG2D/dg_c_solvers.cpp | 460 +++++++++ myf90/ESBGK/Stu_DG2D/dg_cuda_solvers.cu | 550 +++++++++++ myf90/ESBGK/Stu_DG2D/dg_input.txt | 7 + myf90/ESBGK/Stu_DG2D/dg_solvers.f90 | 304 ++++++ myf90/ESBGK/Stu_DG2D/lgl_interp.f90 | 119 +++ myf90/ESBGK/Stu_DG2D/lin_alg.f90 | 166 ++++ myf90/ESBGK/Stu_DG2D/mesh_tools_2d_dg.f90 | 881 ++++++++++++++++++ myf90/ESBGK/Stu_DG2D/problem_setup_2d_dg.f90 | 96 ++ myf90/ESBGK/Stu_DG2D/stuDG.f90 | 378 ++++++++ myf90/ESBGK/Stu_DG2D/stuDG_F90 | Bin 0 -> 94586 bytes myf90/ESBGK/Stu_DG2D/stuDG_driver.m | 46 + .../TestingDGWENO/TestingDGWENO.depend | 2 + .../TestingDGWENO/TestingDGWENO.layout | 16 +- .../TestingDGWENO/configuration.in | 2 +- .../TestingDGWENO/TestingDGWENO/main.f90 | 29 +- .../TestingDGWENO/advection1d/advection1d.cbp | 40 + .../advection1d/advection1d.layout | 4 + .../ESBGK/TestingDGWENO/advection1d/main.f95 | 8 + myf90/advection1d/extramodule.f90 | 20 + myf90/advection1d/main.f90 | 51 + myf90/advection1d/output.plt | 11 + myf90/advection1d/runme | Bin 0 -> 13422 bytes 24 files changed, 3226 insertions(+), 21 deletions(-) create mode 100644 myf90/ESBGK/Stu_DG2D/Makefile create mode 100644 myf90/ESBGK/Stu_DG2D/README.md create mode 100644 myf90/ESBGK/Stu_DG2D/dg_c_solvers.cpp create mode 100644 myf90/ESBGK/Stu_DG2D/dg_cuda_solvers.cu create mode 100644 myf90/ESBGK/Stu_DG2D/dg_input.txt create mode 100644 myf90/ESBGK/Stu_DG2D/dg_solvers.f90 create mode 100644 myf90/ESBGK/Stu_DG2D/lgl_interp.f90 create mode 100644 myf90/ESBGK/Stu_DG2D/lin_alg.f90 create mode 100644 myf90/ESBGK/Stu_DG2D/mesh_tools_2d_dg.f90 create mode 100644 myf90/ESBGK/Stu_DG2D/problem_setup_2d_dg.f90 create mode 100644 myf90/ESBGK/Stu_DG2D/stuDG.f90 create mode 100755 myf90/ESBGK/Stu_DG2D/stuDG_F90 create mode 100644 myf90/ESBGK/Stu_DG2D/stuDG_driver.m create mode 100644 myf90/ESBGK/TestingDGWENO/advection1d/advection1d.cbp create mode 100644 myf90/ESBGK/TestingDGWENO/advection1d/advection1d.layout create mode 100644 myf90/ESBGK/TestingDGWENO/advection1d/main.f95 create mode 100644 myf90/advection1d/extramodule.f90 create mode 100644 myf90/advection1d/main.f90 create mode 100644 myf90/advection1d/output.plt create mode 100755 myf90/advection1d/runme diff --git a/myf90/ESBGK/Stu_DG2D/Makefile b/myf90/ESBGK/Stu_DG2D/Makefile new file mode 100644 index 0000000..9c84d57 --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/Makefile @@ -0,0 +1,53 @@ +FC=gfortran +FC_FLAGS=-O3 +FC_LIBS=-lstdc++ +# -lcudart -lcublas +#FC_LD_LIB=/usr/local/cuda-5.0/lib64 +#FC_INCL=/usr/local/cuda-5.0/include + +CC=g++ +CC_FLAGS=-O3 +CC_LIBS= +CC_LD_LIB= +CC_INCL= + + +#CU_CC=nvcc +#CU_FLAGS=-arch compute_20 +#CU_LIBS=-lcublas + + +TARGET=stuDG_F90 +SOURCE=stuDG.f90 + + +$(TARGET): $(SOURCE) lgl_interp.o mesh_tools_2d_dg.o problem_setup_2d_dg.o \ +lin_alg.o dg_solvers.o dg_c_solvers.o + $(FC) -o $(TARGET) $(SOURCE) $(FC_FLAGS) lgl_interp.o \ + mesh_tools_2d_dg.o problem_setup_2d_dg.o lin_alg.o \ + dg_solvers.o dg_c_solvers.o -L$(FC_LD_LIB) -I$(FC_INCL) $(FC_LIBS) + +lgl_interp.o: lgl_interp.f90 + $(FC) -c lgl_interp.f90 $(FC_FLAGS) + +mesh_tools_2d_dg.o: mesh_tools_2d_dg.f90 + $(FC) -c mesh_tools_2d_dg.f90 $(FC_FLAGS) + +problem_setup_2d_dg.o: problem_setup_2d_dg.f90 + $(FC) -c problem_setup_2d_dg.f90 $(FC_FLAGS) + +lin_alg.o: lin_alg.f90 + $(FC) -c lin_alg.f90 $(FC_FLAGS) + +dg_solvers.o: dg_solvers.f90 + $(FC) -c dg_solvers.f90 $(FC_FLAGS) + +dg_c_solvers.o: dg_c_solvers.cpp + $(CC) -c dg_c_solvers.cpp $(CC_FLAGS) + +#dg_cuda_solvers.o: dg_cuda_solvers.cu +# $(CU_CC) -c dg_cuda_solvers.cu $(CU_LIBS) $(CU_FLAGS) + + +clean: + rm *.o *.mod $(TARGET) diff --git a/myf90/ESBGK/Stu_DG2D/README.md b/myf90/ESBGK/Stu_DG2D/README.md new file mode 100644 index 0000000..1033ef9 --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/README.md @@ -0,0 +1,4 @@ +advection-2D-DG-CUDA +==================== + +a simple 2D linear advection code written in F90 with C, F90 and CUDA solvers. Provided with a MATLAB interface script for problem setup. \ No newline at end of file diff --git a/myf90/ESBGK/Stu_DG2D/dg_c_solvers.cpp b/myf90/ESBGK/Stu_DG2D/dg_c_solvers.cpp new file mode 100644 index 0000000..5a3b5fd --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/dg_c_solvers.cpp @@ -0,0 +1,460 @@ +#include +#include + + +//contains C/C++ subroutines for performing the DG time-stepping for the +//advection equation. Since the core routine will only be called from +//a FORTRAN main program, all array data will be arranged in FORTRAN +//column-major style. + +using namespace std; + +//-------------------- subroutine for computing RHS --------------------- + +void compute_rhs(float * rhs, float * q, float * u, float * v, + float * ksi_x, float * ksi_y, float * eta_x, + float * eta_y, float * jac, float * psi, + float * psi_ksi, float * psi_eta, int nelem, + int npts, int nqs){ + + //local variables + float wq, e_x, e_y, n_x,n_y,h_k,u_k,v_k,q_k,h_e,h_n,dhdx_k,dhdy_k; + + + + //initialize rhs + for(int l = 0;l +#include +#include "device_functions.h" +#include "cublas_v2.h" +#include +#define MAX_NGL 5 +#define TPB_1 96 +#define TPB_2 32 + + +using namespace std; + +//------- compute rhs kernel ---------------------------------------------------------- + + + +// ----- new version with additional ranks of threads at the nqs level and atomic operations + +__global__ void compute_rhs_cuda(float * rhs, float * q, float * u, float * v, + float * ksi_x, float * ksi_y, float * eta_x, + float * eta_y, float * jac, float * psi, + float * psi_ksi, float * psi_eta, int nelem, + int npts, int nqs){ + + int l_tid = threadIdx.x; + int l_block = blockDim.x; + int tid = threadIdx.x+blockIdx.x*blockDim.x; + int k = threadIdx.y; + + if(tid clamr){ + clam = claml; + }else{ + clam = clamr; + } + + // //flux variables + // fxl = qlq_k*ul_k; + // fyl = qlq_k*vl_k; + // fxr = qrq_k*ur_k; + // fyr = qrq_k*vr_k; + + // //normal flux component + // flux_ql = nxl*fxl+nyl*fyl; + // flux_qr = nxr*fxr+nyr*fyr; + // flux_q = flux_ql-flux_qr; + + + flux_q = nxl*qlq_k*ul_k+nyl*qlq_k*vl_k-(nxr*qrq_k*ur_k+nyr*qrq_k*vr_k); + //dissipation term + diss_q = clam*(qrq_k-qlq_k); + + //construct Rusanov flux + flux_q = 0.5*(flux_q-diss_q); + + // //loop through side interpolation points + // cout << "side = " << is+1 << endl; + + for(int i=0;i>>(rhs_d,qp_d,u0_d,v0_d,ksi2d_x_d, + ksi2d_y_d,eta2d_x_d,eta2d_y_d, + jac2d_d,psi2d_d,psi2d_ksi_d, + psi2d_eta_d,nelem,npts,nqs); + + // cout << "After compute_rhs, rhs_d = "<< endl; + // print_q(rhs_d,q_print,nelem,npts); + //compute/apply flux + + + compute_flux_cuda<<>>(rhs_d,qp_d,u0_d,v0_d,psideh_d,nx_d,ny_d, + jac_side_d,psi_d,nside,ngl,nq,imapl_d, + imapr_d,nelem); + + // cout << "After compute_flux, rhs_d = " << endl; + // print_q(rhs_d,q_print,nelem,npts); + + + //apply inverse mass matrix (possibly re-structure inverse mass matrix to allow using cublas/cusparse) + apply_invM_cuda<<>>(rhs_d,invm_d,nelem,npts); + + + + // cout << "After applying inverse mass matrix, rhs_d = " << endl; + // print_q(rhs_d,q_print,nelem,npts); + + // cout << "Before inter-stage update, q0_d = " << endl; + // print_q(q0_d,q_print,nelem,npts); + + // cout << "Before inter-stage update, q1_d = " << endl; + // print_q(q1_d,q_print,nelem,npts); + + // cout << "Before inter-stage update, qp_d = " << endl; + // print_q(qp_d,q_print,nelem,npts); + + // cout << "dtt = " << dtt; + // cout << "a0 = " << a0; + // cout << "a1 = " << a1; + + //perform inter-stage updates + //qp_d = a0*q0_d + a1*q1_d + dtt*rhs_d + stat=cublasScopy(handle,nelem*npts,rhs_d,1,qp_d,1);//qp_d=rhs_d + stat=cublasSscal(handle,nelem*npts,&dtt,qp_d,1);//qp_d=dtt*qp_d + stat=cublasSaxpy(handle,nelem*npts,&a1,q1_d,1,qp_d,1); //qp_d = qp_d+a1*q1_d + stat=cublasSaxpy(handle,nelem*npts,&a0,q0_d,1,qp_d,1);//qp_d = qp_d+a0*q0_d + + // interstage_update<<>>(q0_d,q1_d,qp_d,rhs_d,a0,a1,dtt,nelem,npts); + + + // cout << "After inter-stage update, qp_d = " << endl; + // print_q(qp_d,q_print,nelem,npts); + + //q1_d = qp_d + stat=cublasScopy(handle,nelem*npts,qp_d,1,q1_d,1); + + }//for(int ik... + + //q0_d = qp_d + stat=cublasScopy(handle,nelem*npts,qp_d,1,q0_d,1); + + // cout << "After all stages, q0_d = " << endl; + // print_q(q0_d,q_print,nelem,npts); + + + } + //end time-stepping section, transfer needed output data back from GPU + cudaMemcpy(q0,q0_d,nelem*npts*sizeof(float),cudaMemcpyDeviceToHost); + + //de-allocate memory from GPU + + cublasDestroy(handle); + + cudaFree(rhs_d); + cudaFree(qp_d); + cudaFree(q1_d); + cudaFree(q0_d); + cudaFree(u0_d); + cudaFree(v0_d); + cudaFree(ksi2d_x_d); + cudaFree(ksi2d_y_d); + cudaFree(eta2d_x_d); + cudaFree(eta2d_y_d); + cudaFree(jac2d_d); + cudaFree(psi2d_d); + cudaFree(psi2d_ksi_d); + cudaFree(psi2d_eta_d); + cudaFree(psideh_d); + cudaFree(nx_d); + cudaFree(ny_d); + cudaFree(jac_side_d); + cudaFree(imapl_d); + cudaFree(imapr_d); + cudaFree(invm_d); + + //de-allocate any local memory + delete [] rhs; + delete [] invm; + delete [] q_print; + +} diff --git a/myf90/ESBGK/Stu_DG2D/dg_input.txt b/myf90/ESBGK/Stu_DG2D/dg_input.txt new file mode 100644 index 0000000..accb724 --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/dg_input.txt @@ -0,0 +1,7 @@ +2 +2 +F +1000 +1 +2 +3 diff --git a/myf90/ESBGK/Stu_DG2D/dg_solvers.f90 b/myf90/ESBGK/Stu_DG2D/dg_solvers.f90 new file mode 100644 index 0000000..b4b5d65 --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/dg_solvers.f90 @@ -0,0 +1,304 @@ +module dg_solvers + implicit none + +contains + + subroutine compute_rhs_dg_ntp(rhs,q,u,v,ksi_x,ksi_y,& + eta_x,eta_y,jac,psi,psi_ksi,psi_eta,nelem,npts,nqs) + implicit none + integer,intent(in)::nqs,npts,nelem + real(kind=4),dimension(nelem,npts),intent(in)::q,u,v + real(kind=4),dimension(nelem,nqs),intent(in)::ksi_x,ksi_y,eta_x,eta_y,jac + + real(kind=4),dimension(npts,nqs),intent(in)::psi,psi_ksi,psi_eta + real(kind=4),dimension(nelem,npts),intent(inout)::rhs + character(len=36):: fmt_string,fmt_string_nq + ! local variables and arrays + integer :: ie,k,i + + real(kind=4) :: wq,e_x,e_y,n_x,n_y,h_k,u_k,v_k,q_k,h_e,h_n,dhdx_k,dhdy_k + + ! initialize rhs + + rhs = 0. + do ie = 1,nelem + + ! loop integration points + do k = 1,nqs + wq = jac(ie,k) + e_x = ksi_x(ie,k) + e_y = ksi_y(ie,k) + n_x = eta_x(ie,k) + n_y = eta_y(ie,k) + + !interpolate at integration points + u_k = 0.; v_k = 0.; q_k = 0.; + do i = 1,npts + h_k = psi(i,k) + u_k = u_k + h_k*u(ie,i) + v_k = v_k + h_k*v(ie,i) + q_k = q_k + h_k*q(ie,i) + end do !i + + ! interpolate at integration points..more + do i = 1,npts + h_e = psi_ksi(i,k) + h_n = psi_eta(i,k) + dhdx_k = h_e*e_x+h_n*n_x + dhdy_k = h_e*e_y+h_n*n_y + rhs(ie,i) = rhs(ie,i)+wq*q_k*(dhdx_k*u_k+dhdy_k*v_k) + end do!i + end do !k + end do !ie + +!!$write(fmt_string,*) '(1X,', npts, 'ES12.3)' +!!$write(fmt_string_nq,*) '(1X,',nqs,'ES12.3)' +!!$ +!!$ +!!$write (*,*) 'output from compute RHS:' +!!$write(*,fmt_string) ((rhs(i,ie), ie = 1,npts),i = 1,nelem) + + + + + + end subroutine compute_rhs_dg_ntp + + subroutine compute_flux_dg_ntp(rhs,q,u,v,psideh,& + nx,ny,jac_side,psi,nside,ngl,nq,imapl,imapr,nelem,npts) + implicit none + integer,intent(in)::nside,ngl,nq,nelem,npts + real(kind=4),dimension(nelem,npts),intent(in)::q,u,v + real(kind=4),dimension(nside,nq),intent(in)::nx,ny,jac_side + real(kind=4),dimension(ngl,nq),intent(in)::psi + integer,dimension(4,2,ngl),intent(in)::imapl,imapr + integer,dimension(nside,4),intent(in)::psideh + real(kind=4),dimension(nelem,npts),intent(inout)::rhs + + ! local arrays and variables + real(kind=4),dimension(ngl)::ql,qr,ul,vl + integer::is,iel,ilocl,l,il,jl,kl,ier,ilocr,ir,jr,kr,i + real(kind=4) :: wq,nxl,nyl,nxr,nyr,qlq_k,qrq_k,u_k,v_k,ul_k,vl_k + real(kind=4) :: ur_k,vr_k,unl,unr,claml,clamr,clam,fxl,fyl,fxr,fyr + real(kind=4) :: flux_ql,flux_qr,flux_q,diss_q,h_i + character(len=36):: fmt_string,fmt_string_nq + integer :: ie, nqs + nqs = nq*nq + ! initialize local arrays + ql = 0.; qr = 0.; ul = 0.; vl = 0.; + + do is = 1,nside + !store left side variables + iel = psideh(is,3) + if(iel .ne. -6) then + ilocl = psideh(is,1) + do l = 1,ngl + !get pointers + il = imapl(ilocl,1,l) + jl = imapl(ilocl,2,l) + kl = (jl-1)*ngl+il + !left element + ql(l)=q(iel,kl) + ul(l)=u(iel,kl) + vl(l)=v(iel,kl) + end do!l + + ! store right side variables + ier=psideh(is,4) + if(ier .ne.0) then + ilocr = psideh(is,2) + do l = 1,ngl + ! get pointers + ir = imapr(ilocr,1,l) + jr = imapr(ilocr,2,l) + kr = (jr-1)*ngl + ir + ! right element + qr(l)=q(ier,kr) + end do !l + endif !if ier... + + ! do gauss-lobatto integration + do l = 1,nq + wq = jac_side(is,l) + !store normal vectors + nxl = nx(is,l) + nyl = ny(is,l) + nxr = -nxl + nyr = -nyl + + ! interpolate onto Quadrature points + qlq_k=0.; qrq_k=0.; u_k = 0.; v_k = 0.; + do i = 1,ngl + qlq_k=qlq_k+psi(i,l)*ql(i) + qrq_k=qrq_k+psi(i,l)*qr(i) + u_k=u_k+psi(i,l)*ul(i) + v_k=v_k+psi(i,l)*vl(i) + end do!i + ul_k=u_k; vl_k = v_k; + ur_k = u_k; vr_k = v_k; + + ! computer rusanov flux constant + unl=nxl*ul_k+nyl*vl_k + unr=nxl*ur_k+nyl*vr_k + claml=abs(unl) + clamr=abs(unr) + clam = max(claml,clamr) + + ! flux variables + fxl=qlq_k*ul_k + fyl=qlq_k*vl_k + fxr=qrq_k*ur_k + fyr=qrq_k*vr_k + + !normal flux component + flux_ql=(nxl*fxl+nyl*fyl) + flux_qr=(nxr*fxr+nyr*fyr) + flux_q=flux_ql-flux_qr + + !dissipation term + diss_q = clam*(qrq_k-qlq_k) + + !construct rusanov flux + flux_q=0.5*(flux_q-diss_q) + + !loop through side interpolation points + do i = 1,ngl + h_i=psi(i,l) + + !left side + il=imapl(ilocl,1,i) + jl=imapl(ilocl,2,i) + kl=(jl-1)*ngl+il + rhs(iel,kl)=rhs(iel,kl)-wq*h_i*flux_q + + !right side + if(ier .gt. 0) then + ir = imapr(ilocr,1,i) + jr = imapr(ilocr,2,i) + kr=(jr-1)*ngl+ir + rhs(ier,kr)=rhs(ier,kr)+wq*h_i*flux_q + endif !(ier... + end do !i + end do !l + end if !(iel... + + end do !is + +!!$write(fmt_string,*) '(1X,', npts, 'ES12.3)' +!!$write(fmt_string_nq,*) '(1X,',nqs,'ES12.3)' +!!$ +!!$ +!!$write (*,*) 'output from compute flux:' +!!$write(*,fmt_string) ((rhs(i,ie), ie = 1,npts),i = 1,nelem) + + end subroutine compute_flux_dg_ntp + + + subroutine dg_advection_2d_non_tensor_product(ntime,dt,kstages,q0,u0,v0, & + psi,ksi2d_x,ksi2d_y,eta2d_x,eta2d_y,jac2d,psi2d,psi2d_ksi,psi2d_eta, & + intma2d,psideh, jac_side, imapl,imapr, nx,ny, nelem,npts,nqs,ngl,nq, & + Mmatrix_inv,nside) + implicit none + + ! input arguments + integer, intent(in)::nside,nq,ngl,nqs,npts,nelem,ntime,kstages + real(kind=4),intent(in)::dt + integer,dimension(4,3,ngl),intent(in) :: imapl,imapr + integer,dimension(nside,4),intent(in) :: psideh + integer,dimension(nelem,npts),intent(in)::intma2d + real(kind=4),dimension(ngl,nq),intent(in) :: psi + real(kind=4),dimension(npts,nqs),intent(in)::psi2d,psi2d_ksi,psi2d_eta + real(kind=4),dimension(nelem,nqs),intent(in)::ksi2d_x,ksi2d_y,eta2d_x,eta2d_y + real(kind=4),dimension(nelem,nqs),intent(in)::jac2d + real(kind=4),dimension(nside,nq),intent(in)::jac_side,nx,ny + real(kind=4),dimension(nelem,npts,npts),intent(in)::Mmatrix_inv + ! in/out arguments + real(kind=4),dimension(nelem,npts),intent(inout)::q0,u0,v0 + + ! local variables and arrays + integer :: itime, ik,e,ie,i + real(kind=4),dimension(npts,npts) :: Mtemp + real(kind=4),dimension(npts) :: rtemp + real(kind=4),dimension(nelem,npts)::rhs + real(kind=4),dimension(nelem,npts)::qp,q1 + real(kind=4) :: a0,a1,beta,dtt +character(len=36):: fmt_string,fmt_string_nq +integer, parameter :: time_step_rep_freq = 100 + ! initialize the state vector + rhs = 0.; q1 = q0; qp = q0; + do itime = 1,ntime + + if(mod(itime,time_step_rep_freq) .eq. 0) then + write (*,100) itime +100 format(' ','Commencing time step', I6) + endif + + do ik = 1,kstages + + select case (kstages) + case (2) + select case (ik) + case (1) + a0 = 1.; a1 = 0.; beta = 1.; + case (2) + a0 = 0.5; a1 = 0.5; beta = 0.5; + end select !ik + case (3) + select case (ik) + case (1) + a0 = 1.; a1 = 0.; beta = 1.; + case (2) + a0 = 3./4.; a1 = 1./4.; beta = 1./4. + case (3) + a0 = 1./3.; a1 = 2./3.; beta = 2./3. + end select !ik + end select !kstages + dtt = dt*beta; + + ! compute RHS + call compute_rhs_dg_ntp(rhs,qp,u0,v0,ksi2d_x,ksi2d_y, & + eta2d_x,eta2d_y,jac2d,psi2d,psi2d_ksi,psi2d_eta,& + nelem,npts,nqs) + ! compute/apply flux + call compute_flux_dg_ntp(rhs,qp,u0,v0,psideh,nx,ny,& + jac_side,psi,nside,ngl,nq,imapl,imapr,nelem,npts) + ! apply inverse mass matrix + do e = 1,nelem + Mtemp = Mmatrix_inv(e,:,:) + rtemp = rhs(e,:) + rtemp = matmul(Mtemp,rtemp) + rhs(e,:) = rtemp + end do !e + +!!$write(fmt_string,*) '(1X,', npts, 'ES12.3)' +!!$write(fmt_string_nq,*) '(1X,',nqs,'ES12.3)' +!!$ +!!$ +!!$write (*,*) 'output after inverting mass matrix:' +!!$write(*,fmt_string) ((rhs(i,ie), ie = 1,npts),i = 1,nelem) + + +!!$write(*,*) 'dtt = ',dtt + ! update inter-stage + qp = a0*q0 + a1*q1+dtt*rhs + +!!$write (*,*) 'qp after interstage:' +!!$write(*,fmt_string) ((qp(i,ie), ie = 1,npts),i = 1,nelem) + + + ! update + q1 = qp + +!!$write (*,*) 'q1 after interstage:' +!!$write(*,fmt_string) ((q1(i,ie), ie = 1,npts),i = 1,nelem) + + end do !ik + + ! update q0 + q0 = qp + + end do!itime + + end subroutine dg_advection_2d_non_tensor_product + +end module dg_solvers diff --git a/myf90/ESBGK/Stu_DG2D/lgl_interp.f90 b/myf90/ESBGK/Stu_DG2D/lgl_interp.f90 new file mode 100644 index 0000000..20cd390 --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/lgl_interp.f90 @@ -0,0 +1,119 @@ +module lgl_interp + + implicit none + +contains + + subroutine legendre_poly(L0,L0_1,L0_2,p,x) + integer, intent(in) :: p + real(kind=4), intent(in) :: x + real(kind=4), intent(out) :: L0,L0_1,L0_2 + integer :: i + real(kind=4) :: a, b + real(kind=4) :: L1, L1_1, L1_2,L2,L2_1,L2_2 + + L1 = 0.; L1_1 = 0.; L1_2 = 0.; + L0 = 1.; L0_1 = 0.; L0_2 = 0.; + + do i= 1,p + L2 = L1; L2_1 = L1_1; L2_2 = L1_2; + L1 = L0; L1_1 = L0_1; L1_2 = L0_2; + a = (2.0*i-1.0)/i + b = (i-1.)/i + L0 = a*x*L1 - b*L2 + L0_1 = a*(L1 + x*L1_1) - b*L2_1 + L0_2 = a*(2*L1_1 + x*L1_2) - b*L2_2 + + end do + + end subroutine legendre_poly + + subroutine legendre_gauss_lobatto(xgl, wgl, ngl) + integer, intent(in) :: ngl + real(kind=4),dimension(ngl), intent(out) :: xgl,wgl + integer :: p, ph + integer i, k + real(kind=4) :: x, dx + real(kind=4) :: pi + real(kind=4) :: L0,L0_1,L0_2 + + pi = 4.*atan(1.0) + p = ngl - 1 + ph = floor((p+1.)/2.) + xgl = 0.0; wgl = 0.0; + + do i = 1,ph + x = cos((2.0*i-1)*pi/(2.0*p+1)) + do k = 1,20 + call legendre_poly(L0,L0_1,L0_2,p,x) +! get new Newton Iteration + dx = -(1-x**2)*L0_1/(-2*x*L0_1 + (1-x**2)*L0_2) + x = x+dx + if (abs(dx) .lt. 1.0e-20) then + exit ! end the loop + endif + end do + xgl(p+2-i)=x + wgl(p+2-i)=2./(p*(p+1)*(L0**2)) + end do + + + ! check for zero root + if ((p+1) .ne. 2*ph) then + x=0.0 + call legendre_poly(L0,L0_1,L0_2,p,x) + xgl(ph+1)=x + wgl(ph+1)=2./(p*(p+1)*(L0**2)) + + endif + + ! find remainder of roots via symmetry + + do i = 1,ph + xgl(i) = -xgl(p+2-i) + wgl(i) = wgl(p+2-i) + end do + + + + end subroutine legendre_gauss_lobatto + + subroutine lagrange_basis(psi, dpsi, xnq, wnq,ngl,nq,xgl) + integer, intent(in) :: ngl, nq + real(kind=4), dimension(ngl),intent(in) :: xgl + real(kind=4), dimension(nq), intent(inout) :: xnq, wnq + real(kind=4), dimension(ngl,nq), intent(inout) :: psi, dpsi + integer :: l, i, j, k + real(kind=4) :: xl, xi, xj, ddpsi, xk + + call legendre_gauss_lobatto(xnq,wnq,nq) + psi = 0.0; dpsi = 0.0; + !Perform Quadrature + do l = 1,nq + xl = xnq(l) + do i = 1,ngl + xi = xgl(i) + psi(i,l)=1.0 + dpsi(i,l) = 0.0 + do j = 1,ngl + xj = xgl(j) + if (j .ne. i) then + psi(i,l) = psi(i,l)*(xl-xj)/(xi-xj) + endif !(j .ne. i) + ddpsi = 1. + if (j .ne. i) then + do k = 1,ngl + xk = xgl(k) + if((k .ne. i) .and. (k .ne. j)) then + ddpsi = ddpsi*(xl-xk)/(xi-xk) + endif !((k .ne. i) ... + end do ! k + dpsi(i,l) = dpsi(i,l)+ddpsi/(xi-xj) + endif ! (j .ne. i)... + end do ! j + end do ! i + end do ! l + + end subroutine lagrange_basis + +end module lgl_interp diff --git a/myf90/ESBGK/Stu_DG2D/lin_alg.f90 b/myf90/ESBGK/Stu_DG2D/lin_alg.f90 new file mode 100644 index 0000000..760ba7e --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/lin_alg.f90 @@ -0,0 +1,166 @@ +module lin_alg + implicit none + + +contains + + subroutine diag_matrix_invert(A,A_inv,n) + ! of course this is dumb...but until I get something put together for + ! non-diagonal matrices, this will have to do. + integer, intent(in) :: n + real(kind=4),dimension(n,n),intent(in) :: A + real(kind=4),dimension(n,n),intent(inout) :: A_inv + + integer :: i + real(kind=4) :: a_tmp + + A_inv = 0. + + do i = 1,n + a_tmp = A(i,i) + A_inv(i,i) = 1./a_tmp + end do + + end subroutine diag_matrix_invert + + subroutine matrix_invert(A,A_inv,n) + integer, intent(in) :: n + real(kind=4),dimension(n,n),intent(inout) :: A + real(kind=4),dimension(n,n),intent(inout) :: A_inv + + ! local variables + integer, dimension(n) :: indx + + call MIGS(A,n,A_inv,indx) + + + + end subroutine matrix_invert + + ! Updated 10/24/2001. + ! +!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! ! + ! Please Note: ! + ! ! + ! (1) This computer program is written by Tao Pang in conjunction with ! + ! his book, "An Introduction to Computational Physics," published ! + ! by Cambridge University Press in 1997. ! + ! ! + ! (2) No warranties, express or implied, are made for this program. ! + ! ! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! + SUBROUTINE MIGS (A,N,X,INDX) + ! + ! Subroutine to invert matrix A(N,N) with the inverse stored + ! in X(N,N) in the output. Copyright (c) Tao Pang 2001. + ! + IMPLICIT NONE + INTEGER, INTENT (IN) :: N + INTEGER :: I,J,K + INTEGER, INTENT (OUT), DIMENSION (N) :: INDX + REAL(kind=4), INTENT (INOUT), DIMENSION (N,N):: A + REAL(kind=4), INTENT (OUT), DIMENSION (N,N):: X + REAL(kind=4), DIMENSION (N,N) :: B + ! + DO I = 1, N + DO J = 1, N + B(I,J) = 0.0 + END DO + END DO + DO I = 1, N + B(I,I) = 1.0 + END DO + ! + CALL ELGS (A,N,INDX) + ! + DO I = 1, N-1 + DO J = I+1, N + DO K = 1, N + B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K) + END DO + END DO + END DO + ! + DO I = 1, N + X(N,I) = B(INDX(N),I)/A(INDX(N),N) + DO J = N-1, 1, -1 + X(J,I) = B(INDX(J),I) + DO K = J+1, N + X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I) + END DO + X(J,I) = X(J,I)/A(INDX(J),J) + END DO + END DO + END SUBROUTINE MIGS + ! + SUBROUTINE ELGS (A,N,INDX) + ! + ! Subroutine to perform the partial-pivoting Gaussian elimination. + ! A(N,N) is the original matrix in the input and transformed matrix + ! plus the pivoting element ratios below the diagonal in the output. + ! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001. + ! + IMPLICIT NONE + INTEGER, INTENT (IN) :: N + INTEGER :: I,J,K,ITMP + INTEGER, INTENT (OUT), DIMENSION (N) :: INDX + REAL(kind=4) :: C1,PI,PI1,PJ + REAL(kind=4), INTENT (INOUT), DIMENSION (N,N) :: A + REAL(kind=4), DIMENSION (N) :: C + ! + ! Initialize the index + ! + DO I = 1, N + INDX(I) = I + END DO + ! + ! Find the rescaling factors, one from each row + ! + DO I = 1, N + C1= 0.0 + DO J = 1, N + !C1 = AMAX1(C1,ABS(A(I,J))) + C1 = max(C1,ABS(A(I,J))) + END DO + C(I) = C1 + END DO + ! + ! Search the pivoting (largest) element from each column + ! + DO J = 1, N-1 + PI1 = 0.0 + DO I = J, N + PI = ABS(A(INDX(I),J))/C(INDX(I)) + IF (PI.GT.PI1) THEN + PI1 = PI + K = I + ENDIF + END DO + ! + ! Interchange the rows via INDX(N) to record pivoting order + ! + ITMP = INDX(J) + INDX(J) = INDX(K) + INDX(K) = ITMP + DO I = J+1, N + PJ = A(INDX(I),J)/A(INDX(J),J) + ! + ! Record pivoting ratios below the diagonal + ! + A(INDX(I),J) = PJ + ! + ! Modify other elements accordingly + ! + DO K = J+1, N + A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K) + END DO + END DO + END DO + ! + END SUBROUTINE ELGS + +end module lin_alg diff --git a/myf90/ESBGK/Stu_DG2D/mesh_tools_2d_dg.f90 b/myf90/ESBGK/Stu_DG2D/mesh_tools_2d_dg.f90 new file mode 100644 index 0000000..7e42818 --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/mesh_tools_2d_dg.f90 @@ -0,0 +1,881 @@ +module meshtools_2d_dg + + implicit none + +contains + + subroutine create_grid_2d_dg(coord,intma,bsido, & + iperiodic,xgl, npoin, nelem, nboun, & + nelx, nely, ngl) + implicit none + integer, intent(in) :: npoin, nelem, nboun + integer, intent(in) :: nelx, nely, ngl + real(kind=4),dimension(ngl),intent(in) :: xgl + real(kind=4),dimension(npoin,2), intent(inout) :: coord + integer, dimension(nelem,ngl,ngl),intent(inout)::intma + integer, dimension(nboun,4), intent(inout)::bsido + integer, dimension(npoin), intent(inout)::iperiodic + + integer, dimension(npoin,npoin)::node + + real(kind=4) :: xmin, xmax, ymin, ymax, dx, dy + real(kind=4) :: x,y,x0,y0 + integer :: nop, nx, ny + integer :: ip,jj,j,k,l,l1,l2,ii,i,j1,ie,ib,i1,i2,ip1,ip2 + +!write(*,*) 'entering the subroutine' + node = 0 + + xmin = -1.; xmax = 1.; ymin = -1.; ymax = 1.; + dx = (xmax-xmin)/real(nelx) + dy = (ymax-ymin)/real(nely) + nop = ngl-1; + nx = nelx*nop+1 + ny = nely*nop+1 + +!write(*,*) 'generate coord' + ! generate coord + ip = 0; jj = 0; + do k = 1,nely + y0=ymin+real(k-1)*dy + if(k .eq. 1) then + l1 = 1 + else + l1 = 2 + endif !if(k .eq... + + do l = l1,ngl + y = (xgl(l)+1.)*dy/2. + y0 + jj=jj+1 + ii=0 + do i=1,nelx + x0 = xmin + real(i-1)*dx + if (i .eq. 1) then + j1 = 1 + else + j1 = 2 + endif !if(i .eq. ... + + do j = j1,ngl + ii = ii+1 + ip = ip+1 + x = (xgl(j)+1.)*dx/2. + x0 + coord(ip,1) = x + coord(ip,2) = y + node(ii,jj) = ip + end do !j + end do !i + end do !l + end do !k + +!write(*,*) 'generating intma' + + ! generate intma + ie = 0 + do k = 1,nely + do i = 1,nelx + ie = ie+1 + do l = 1,ngl + jj = (ngl-1)*(k-1)+l + do j = 1,ngl + ii = (ngl-1)*(i-1)+j + ip = node(ii,jj) + intma(ie,j,l) = ip + end do !j + end do !l + end do !i + end do !k + +!write(*,*) 'generating bsido' + ! generate bsido + ib = 0 + do i = 1,nelx + ie = i + ib = ib+1 + i1 = (i-1)*(ngl-1)+1 + i2 = (i-1)*(ngl-1) + ngl + ip1 = node(i1,1) + ip2 = node(i2,1) + bsido(ib,1)=ip1 + bsido(ib,2)=ip2 + bsido(ib,3)=ie + bsido(ib,4)=6 + end do + + ! right boundary + do i = 1,nely + ie = nelx*i + ib = ib+1 + i1=(i-1)*(ngl-1)+1 + i2=(i-1)*(ngl-1)+ngl + ip1=node(nx,i1) + ip2=node(nx,i2) + bsido(ib,1)=ip1 + bsido(ib,2)=ip2 + bsido(ib,3)=ie + bsido(ib,4)=6 + end do + + ! top boundary + do i = nelx,1,-1 + ie=nelem-(nelx-i) + ib=ib+1 + i1=(i-1)*(ngl-1)+ngl + i2=(i-1)*(ngl-1)+1 + ip1=node(i1,ny) + ip2=node(i2,ny) + bsido(ib,1)=ip1 + bsido(ib,2)=ip2 + bsido(ib,3)=ie + bsido(ib,4)=6 + end do + + ! left boundary + do i = nely,1,-1 + ie=nelx*(i-1)+1 + ib=ib+1 + i1=(i-1)*(ngl-1)+ngl + i2=(i-1)*(ngl-1)+1 + ip1=node(1,i1) + ip2=node(1,i2) + bsido(ib,1)=ip1 + bsido(ib,2)=ip2 + bsido(ib,3)=ie + bsido(ib,4)=6 + end do +!write(*,*) 'periodicity' + + ! periodicity + do i = 1,npoin + iperiodic(i) = i + end do + +!write(*,*) 'x-periodicity' + ! x-periodicity + do i = 1,ny + i1 = node(1,i) + i2 = node(nx,i) + iperiodic(i2) = i1 + end do +!write(*,*) 'y-periodicity' + ! y-periodicity + do i = 1,nx + i1 = node(i,1) + i2 = node(i,ny) + iperiodic(i2)=i1 + end do + + end subroutine create_grid_2d_dg + + subroutine map_deriv_2d(f_ksi,f_eta,psi,dpsi,f,ngl,nq) + implicit none + integer, intent(in) :: ngl,nq + real(kind=4), dimension(ngl,ngl),intent(in) :: f + real(kind=4), dimension(ngl,nq),intent(in)::psi,dpsi + real(kind=4), dimension(nq,nq),intent(inout) :: f_ksi,f_eta + + integer :: l,k,i,j + real(kind=4) :: sum_ksi, sum_eta + + f_ksi = 0.; f_eta = 0. + + do l = 1,nq + do k = 1,nq + sum_ksi = 0. + sum_eta = 0. + + do j = 1,ngl + do i = 1,ngl + sum_ksi = sum_ksi + dpsi(i,k)*psi(j,l)*f(i,j) + sum_eta = sum_eta + psi(i,k)*dpsi(j,l)*f(i,j) + end do ! i + end do ! j + f_ksi(k,l) = sum_ksi + f_eta(k,l) = sum_eta + + end do ! k + end do ! l + + + + end subroutine map_deriv_2d + + + subroutine metrics_2d(ksi_x,ksi_y,eta_x,eta_y,jac, & + coord,intma,psi,dpsi, & + wnq,npoin,nelem,ngl,nq) + + implicit none + integer, intent(in) :: ngl, nq, nelem, npoin + real(kind=4),dimension(nq), intent(in) :: wnq + integer, dimension(nelem,ngl,ngl), intent(in) :: intma + real(kind=4),dimension(ngl,nq),intent(in) :: psi, dpsi + real(kind=4),dimension(npoin,2),intent(in) :: coord + real(kind=4),dimension(nelem,nq,nq),intent(inout):: ksi_x,ksi_y,eta_x,eta_y,jac + + ! local helper arrays and index variables + real(kind=4),dimension(ngl,ngl) :: x, y + real(kind=4),dimension(nq,nq) :: x_ksi, x_eta, y_ksi,y_eta + integer :: ie, j, i, ip + real(kind=4) :: xjac + + ! initialize the variables + x_ksi = 0.; x_eta = 0.; + y_ksi = 0.; y_eta = 0.; + ksi_x = 0.; ksi_y = 0.; + eta_x = 0.; eta_y = 0.; + jac = 0.; + + x = 0. + y = 0. + + + + ! loop thru the elements + do ie = 1,nelem + + ! store element variables + do j = 1,ngl + do i = 1,ngl + ip = intma(ie,i,j) + x(i,j) = coord(ip,1) + y(i,j) = coord(ip,2) + end do !i + end do !j + + + ! construct mapping derivatives + call map_deriv_2d(x_ksi,x_eta,psi,dpsi,x,ngl,nq) + call map_deriv_2d(y_ksi,y_eta,psi,dpsi,y,ngl,nq) + + + ! construct inverse mapping + do j = 1,nq + do i = 1,nq + xjac = x_ksi(i,j)*y_eta(i,j) - x_eta(i,j)*y_ksi(i,j) + ksi_x(ie,i,j)=1./xjac*y_eta(i,j) + ksi_y(ie,i,j)=-1./xjac*x_eta(i,j) + eta_x(ie,i,j)=-1./xjac*y_ksi(i,j) + eta_y(ie,i,j)=1./xjac*x_ksi(i,j) + jac(ie,i,j)=wnq(i)*wnq(j)*abs(xjac) + end do !i + end do !j + end do !ie + + end subroutine metrics_2d + + subroutine create_side_2d(iside,jeside,intma,bsido, & + npoin,nelem,nboun,nside,ngl) + implicit none + ! arguments + integer, intent(in) :: npoin,nelem,nboun,nside,ngl + integer,dimension(nboun,4),intent(in) :: bsido + integer,dimension(nelem,ngl,ngl),intent(in)::intma + integer,dimension(nside,4),intent(inout) ::iside + integer,dimension(nelem,4),intent(inout) :: jeside + + ! local arrays + integer, dimension(npoin) :: lwher,lhowm + integer,dimension(5*npoin) :: icone + integer,dimension(4)::inode,jnode + ! a whole butt-ton of local integer variables and counters. + integer ::in,ie,ip,iloca,iloc1,iele,iwher,ip1,in1,ipt,j + integer :: iold,in2,ip2,is,jloca,il,ir,is1,is2,i1,i2,iel,jnod + integer :: ib,ibe,ibc,ilb,irb,ier + + ! initialize some local arrays + lwher = 0 + lhowm = 0 + icone = 0 + inode = 0 + jnode = 0 + + ! fix lnode + inode(1) = 1 + inode(2) = ngl + inode(3) = ngl + inode(4) = 1 + jnode(1) = 1 + jnode(2) = 1 + jnode(3) = ngl + jnode(4) = ngl + + + + ! count how many elements own each node + do in =1,4 + do ie = 1,nelem + ip = intma(ie,inode(in),jnode(in)) + lhowm(ip)=lhowm(ip)+1 + end do + end do + + + ! track elements owning each node + lwher(1) = 0 + do ip =2,npoin + lwher(ip)=lwher(ip-1)+lhowm(ip-1) + end do + + + + !another tracker array + lhowm = 0 + + do in = 1,4 + do ie = 1,nelem + + ip = intma(ie,inode(in),jnode(in)) + + lhowm(ip)=lhowm(ip)+1 + + jloca=lwher(ip)+lhowm(ip) + + icone(jloca)=ie + end do + end do + + + ! loop over the nodes + iloca = 0 + do ip = 1,npoin + iloc1 = iloca + iele=lhowm(ip) + if(iele .ne. 0) then + iwher=lwher(ip) + + + ! loop over those elements surrounding node ip + ip1 = ip + do iel = 1,iele + ie = icone(iwher+iel) + ! find out position of ip in intma + do in = 1,4 + in1 = in + ipt = intma(ie,inode(in),jnode(in)) + if(ipt .eq. ip) then + exit + endif + end do !in + !check edge element ie which claims ip + j = 0 + do jnod = 1,3,2 + iold = 0 + j=j+1 + in2 =in+jnod + if(in2 .gt. 4) then + in2 = in2-4 + endif + ip2 = intma(ie,inode(in2),jnode(in2)) + if(ip2 .ge. ip1) then + ! check whether side is old or new + if(iloca .ne. iloc1) then + do is = (iloc1+1),iloca + ! iside(is,2) + jloca=is + if(iside(is,2).eq.ip2) then + iold=1 + exit + endif + enddo !is + endif !iloca + + if(iold .eq. 0) then + iloca=iloca+1 + iside(iloca,1)=ip1 + iside(iloca,2)=ip2 + iside(iloca,2+j)=ie + elseif (iold .eq. 1) then + iside(jloca,2+j)=ie + endif!iold + endif!ip2 + enddo!jnod + enddo !iel + + !perform some shifting to order the nodes of a side in CCW + do is = (iloc1+1),iloca + if(iside(is,3).eq. 0) then + iside(is,3)=iside(is,4) + iside(is,4) = 0 + iside(is,1) = iside(is,2) + iside(is,2) = ip1 + endif!iside + enddo!is + endif !if iele + enddo !ip + + if(iloca .ne. nside) then + write (*,*) 'Error in SIDE. iloca nside = ', iloca, nside + endif + + !reset the boundary markers + do is = 1,nside + if(iside(is,4) .eq. 0) then + il=iside(is,1) + ir = iside(is,2) + ie = iside(is,3) + do ib = 1,nboun + ibe = bsido(ib,3) + ibc=bsido(ib,4) + if(ibe .eq. ie) then + ilb = bsido(ib,1) + irb = bsido(ib,2) + if((ilb .eq. il) .and. (irb .eq. ir)) then + iside(is,4) = -ibc + exit + endif !ilb + endif !ibe + enddo !ib + endif !iside + enddo ! is + + + ! form element/side connectivity array + do is = 1,nside + iel = iside(is,3) + ier = iside(is,4) + is1=iside(is,1) + is2=iside(is,2) + !left side + do in = 1,4 + i1 = intma(iel,inode(in),jnode(in)) + in1 = in+1 + if(in1 .gt. 4) then + in1 = 1 + endif + i2 = intma(iel,inode(in1),jnode(in1)) + if ((is1 .eq. i1) .and. (is2 .eq. i2)) then + jeside(iel,in) = is + endif !is1 + enddo !in + + ! right side + if(ier .gt. 0) then + do in = 1,4 + i1 = intma(ier,inode(in),jnode(in)) + in1 = in+1 + if(in1 .gt. 4) then + in1 = 1 + endif + i2 = intma(ier,inode(in1),jnode(in1)) + if((is1 .eq. i2) .and. (is2 .eq. i1)) then + jeside(ier,in) = is + endif !is1 + enddo !in + endif !ier + enddo !is + + end subroutine create_side_2d + + subroutine create_side_2d_dg(psideh,imapl,imapr,iside, & + intma,nside,nelem,ngl) + implicit none + integer, intent(in) :: ngl,nelem,nside + integer, dimension(nelem,ngl,ngl),intent(in)::intma + integer, dimension(nside,4),intent(in) :: iside + integer, dimension(4,2,ngl),intent(inout) :: imapr, imapl + integer, dimension(nside,4), intent(inout) :: psideh + + ! local arrays + integer, dimension(4):: inode, jnode + + ! local index variables and other temps + integer :: l,i,ip1,ip2,iel,ier,j,j1,j2,jp1,jp2 + + psideh = 0 + imapl = 0 + imapr = 0 + + inode = 0 + jnode = 0 + + ! construct boundary pointer + inode(1) = 1 + inode(2) = ngl + inode(3) = ngl + inode(4) = 1 + jnode(1) = 1 + jnode(2) = 1 + jnode(3) = ngl + jnode(4) = ngl + + ! construct IMAP arrays + do l = 1,ngl + imapl(1,1,l) = l + imapl(1,2,l) = 1 + imapr(1,1,l) = ngl+1-l + imapr(1,2,l) = 1 + + imapl(2,1,l) = ngl + imapl(2,2,l) = l + imapr(2,1,l)=ngl + imapr(2,2,l)=ngl+1-l + + imapl(3,1,l) = ngl+1-l + imapl(3,2,l)=ngl + imapr(3,1,l)=l + imapr(3,2,l)=ngl + + imapl(4,1,l)=1 + imapl(4,2,l)=ngl+1-l + imapr(4,1,l)=1 + imapr(4,2,l)=l + end do + + ! loop through the sides + do i = 1,nside + ip1=iside(i,1) + ip2=iside(i,2) + iel=iside(i,3) + ier=iside(i,4) + + ! check for position on left element + do j=1,4 + j1=j + j2=j+1 + if(j2 .gt. 4) then + j2 = 1 + end if + jp1=intma(iel,inode(j1),jnode(j1)) + jp2=intma(iel,inode(j2),jnode(j2)) + + if((ip1 .eq. jp1) .and. (ip2 .eq. jp2)) then + psideh(i,1)=j + exit + end if + end do + + !check for position on Right Element + if(ier .gt. 0) then + do j = 1,4 + j1 = j + j2 = j+1 + if(j2 .gt. 4) then + j2 = 1 + endif + jp1 = intma(ier,inode(j1),jnode(j1)) + jp2 = intma(ier,inode(j2),jnode(j2)) + + if((ip1 .eq. jp2) .and. (ip2 .eq. jp1)) then + psideh(i,2)=j + + exit + endif + end do + endif + + ! store elements into PSIDEH + psideh(i,3) = iel + psideh(i,4) = ier + + end do! i + + + end subroutine create_side_2d_dg + + subroutine compute_normals_2d(nx,ny,jac_side,psideh,& + intma,coord,nside,ngl,nq,wnq,psi,dpsi,nelem,npoin) + implicit none + + integer, intent(in) :: nq,ngl,nside,nelem,npoin + real(kind=4),dimension(ngl,nq),intent(in) :: psi,dpsi + real(kind=4),dimension(nq),intent(in):: wnq + real(kind=4),dimension(npoin,2),intent(in)::coord + integer,dimension(nelem,ngl,ngl),intent(in)::intma + integer,dimension(nside,4),intent(in)::psideh + real(kind=4),dimension(nside,nq)::nx,ny,jac_side + + ! local arrays + real(kind=4),dimension(ngl,ngl) :: x,y + real(kind=4),dimension(nq,nq) :: x_ksi,x_eta,y_ksi,y_eta + + ! local aux variables + integer :: is, ilocl,ilocr,iel,ier,i,j,ip,l + real(kind=4) :: wq, rnx + + ! initialize variables + nx = 0.; ny = 0.; jac_side = 0.; + x = 0.; y = 0.; x_ksi = 0.; x_eta = 0.; + y_ksi = 0.; y_eta = 0. + + + ! loop thru the sides + do is = 1,nside + + ! store left and right elements + ilocl=psideh(is,1) + ilocr=psideh(is,2) + iel=psideh(is,3) + ier=psideh(is,4) + + ! store element variables + do j = 1,ngl + do i = 1,ngl + ip = intma(iel,i,j) + x(i,j) = coord(ip,1) + y(i,j) = coord(ip,2) + end do + end do + + ! construct mapping derivatives + call map_deriv_2d(x_ksi,x_eta,psi,dpsi,x,ngl,nq) + call map_deriv_2d(y_ksi,y_eta,psi,dpsi,y,ngl,nq) + + ! compute normals + do l = 1,nq + wq = wnq(l) + select case (ilocl) + case (1) + i = l + j = 1 + nx(is,l)=y_ksi(i,j) + ny(is,l)=-x_ksi(i,j) + jac_side(is,l)=wq*sqrt(x_ksi(i,j)**2 + y_ksi(i,j)**2) + + case (2) + i = nq + j = l + nx(is,l)=y_eta(i,j) + ny(is,l)=-x_eta(i,j) + jac_side(is,l)=wq*sqrt(x_eta(i,j)**2+y_eta(i,j)**2) + + case (3) + i=nq+1-l + j=nq + nx(is,l)=-y_ksi(i,j) + ny(is,l)=x_ksi(i,j) + jac_side(is,l)=wq*sqrt(x_ksi(i,j)**2+y_ksi(i,j)**2) + + case (4) + i=1 + j=nq+1-l + nx(is,l)=-y_eta(i,j) + ny(is,l)=x_eta(i,j) + jac_side(is,l)=wq*sqrt(x_eta(i,j)**2+y_eta(i,j)**2) + + + end select + end do !l + + ! normalize normals + do l = 1,nq + rnx=sqrt(nx(is,l)**2 + ny(is,l)**2) + nx(is,l)=nx(is,l)/rnx + ny(is,l)=ny(is,l)/rnx + end do + end do ! is + + end subroutine compute_normals_2d + + subroutine create_periodicity_2d(psideh,iside,coord,nside,nboun,npoin) + implicit none + integer, intent(in) :: nside,nboun,npoin + real(kind=4),dimension(npoin,2),intent(in)::coord + integer,dimension(nside,4),intent(inout)::iside + integer,dimension(nside,4),intent(inout)::psideh + + ! local arrays + real(kind=4),parameter :: tol = 1e-6 + + integer, dimension(nboun/4) :: ileft,iright,itop,ibot + ! note the hack...I believe this will only really work on domains where nelx = nely + + integer :: nleft,nright,ntop,nbot + + real(kind=4) :: xmax,xmin,ymax,ymin,xm,ym + + integer :: j,is,ier,i1,i2,i,isr,isl +real(kind=4) :: yl1,yr2,xl1,xr2 + + ! initialize + nleft = 0; nright = 0; ntop = 0; nbot = 0; + ileft = 0; iright = 0; itop = 0; ibot = 0; + + ! find extrema of domain + xmax = maxval(coord(:,1)); xmin = minval(coord(:,1)); + ymax = maxval(coord(:,2)); ymin = minval(coord(:,2)); + + !loop thru sides and extract Left, Right, Bot and top + do is = 1,nside + ier = psideh(is,4) + if(ier .eq. -6) then + i1 = iside(is,1); i2 = iside(is,2) + xm = 0.5*(coord(i1,1)+coord(i2,1)) + ym = 0.5*(coord(i1,2)+coord(i2,2)) + + ! check grid point + if(abs(xm-xmin) .lt. tol) then + nleft = nleft + 1 + ileft(nleft) = is + elseif (abs(xm-xmax) .lt. tol) then + nright = nright + 1 + iright(nright)=is + elseif(abs(ym-ymin) .lt. tol) then + nbot=nbot+1 + ibot(nbot)=is + elseif (abs(ym-ymax) .lt. tol) then + ntop = ntop + 1 + itop(ntop) = is + else + write (*,*) 'No match in PERIODIC_BCS for is ier = ' + write (*,*) is, ier + endif + endif !ier + enddo + +!!$write(*,*) 'nbot = ', nbot +!!$write(*,*) 'ntop = ', ntop +!!$write(*,*) 'nleft = ', nleft +!!$write(*,*) 'nright = ', nright +!!$ +!!$write (*,*) 'itop = ', itop +!!$write (*,*) 'ibot = ', ibot +!!$write (*,*) 'ileft = ', ileft +!!$write (*,*) 'iright = ', iright + + ! loop through periodic BCs + ! first - do left and right + + do i = 1,nleft + isl=ileft(i) + i1=iside(isl,1) + yl1=coord(i1,2) + !search for corresponding right edge + do j=1,nright + isr=iright(j) + i2=iside(isr,2) + yr2=coord(i2,2) + if(abs(yl1-yr2) .lt. tol) then +!!$write(*,*) 'found edge. isl=',isl,'isr=',isr + psideh(isl,2)=psideh(isr,1) + psideh(isl,4)=psideh(isr,3) + psideh(isr,3)=-6 + iside(isl,4)=iside(isr,3) + iside(isr,3)=-6 + exit + endif !abs(yl2... + end do ! j + end do !i + + ! second, do top and bottom + do i = 1,nbot + isl=ibot(i) + i1=iside(isl,1) + xl1=coord(i1,1) + + ! search for corresponding top edge + do j=1,ntop + isr = itop(j) + i2=iside(isr,2) + xr2=coord(i2,1) + if(abs(xl1-xr2) .lt. tol) then +!!$write(*,*) 'found edge. isl=',isl,'isr=',isr + psideh(isl,2)=psideh(isr,1) + psideh(isl,4)=psideh(isr,3) + psideh(isr,3)=-6 + iside(isl,4)=iside(isr,3) + iside(isr,3)=-6 + exit + endif + enddo + enddo + + + end subroutine create_periodicity_2d + + subroutine reshape_1d_to_2d(psi2d,psi2d_ksi,psi2d_eta,& + ksi2d_x,ksi2d_y,eta2d_x,eta2d_y, & + jac2d,intma2d,npts,nqs,psi,dpsi,& + ksi_x,ksi_y,eta_x,eta_y,jac,intma,& + nelem,ngl,nq) + integer,intent(in) :: nq,ngl,nelem,npts,nqs + real(kind=4),dimension(nelem,nq,nq),intent(in)::ksi_x,ksi_y,eta_x,eta_y,jac + real(kind=4),dimension(ngl,nq),intent(in)::psi,dpsi + integer,dimension(nelem,ngl,ngl),intent(in)::intma + real(kind=4),dimension(npts,nqs),intent(inout)::psi2d,psi2d_ksi,psi2d_eta + real(kind=4),dimension(nelem,nqs),intent(inout)::ksi2d_x,ksi2d_y + real(kind=4),dimension(nelem,nqs),intent(inout)::eta2d_x,eta2d_y + real(kind=4),dimension(nelem,nqs),intent(inout)::jac2d + integer,dimension(nelem,npts),intent(inout)::intma2d + + ! local variables + integer :: e,k,j,i,n,m,l + + ! initialize arrays + psi2d=0.; psi2d_ksi =0.; psi2d_eta =0. + + do k=1,ngl + do j = 1,ngl + i=(k-1)*ngl+j + do n = 1,nq + do m = 1,nq + l=(n-1)*nq + m + psi2d(i,l)=psi(j,m)*psi(k,n) + psi2d_ksi(i,l)=dpsi(j,m)*psi(k,n) + psi2d_eta(i,l)=psi(j,m)*dpsi(k,n) + end do !m + end do !n + end do!j + end do !k + + !initialize the next set of arrays + ksi2d_x = 0.; ksi2d_y = 0.; + eta2d_x = 0.; eta2d_y = 0.; + jac2d = 0. + do e = 1,nelem + do k = 1,nq + do j = 1,nq + i=(k-1)*nq+j + ksi2d_x(e,i)=ksi_x(e,j,k) + ksi2d_y(e,i)=ksi_y(e,j,k) + eta2d_x(e,i)=eta_x(e,j,k) + eta2d_y(e,i)=eta_y(e,j,k) + jac2d(e,i)=jac(e,j,k) + end do !j + end do !k + end do !e + + ! intma + intma2d = 0 + do e=1,nelem + do k = 1,ngl + do j = 1,ngl + i=(k-1)*ngl+j + intma2d(e,i)=intma(e,j,k) + end do !j + end do !k + end do !e + + end subroutine reshape_1d_to_2d + + subroutine create_Mmatrix2d_dg_NonTensorProduct(Mmatrix, & + jac2d,psi2d,nelem,npts,nqs) + integer, intent(in) :: nqs,npts,nelem + real(kind=4),dimension(npts,nqs),intent(in) :: psi2d + real(kind=4),dimension(nelem,nqs),intent(in)::jac2d + real(kind=4),dimension(nelem,npts,npts),intent(inout) :: Mmatrix + + ! local variables + integer :: ie,k,i,j + real(kind=4) :: wq,h_i,h_j + + ! initialize Mmatrix + Mmatrix = 0. + + do ie = 1,nelem + do k = 1,nqs + wq = jac2d(ie,k) + do i = 1,npts + h_i = psi2d(i,k) + do j = 1,npts + h_j = psi2d(j,k) + Mmatrix(ie,i,j)=Mmatrix(ie,i,j)+wq*h_i*h_j + end do ! j + end do ! i + end do ! k + end do ! ie + + + end subroutine create_Mmatrix2d_dg_NonTensorProduct + +end module meshtools_2d_dg diff --git a/myf90/ESBGK/Stu_DG2D/problem_setup_2d_dg.f90 b/myf90/ESBGK/Stu_DG2D/problem_setup_2d_dg.f90 new file mode 100644 index 0000000..a67325a --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/problem_setup_2d_dg.f90 @@ -0,0 +1,96 @@ +module problem_setup_2d_dg + implicit none + +contains + + + subroutine exact_solution(qe,ue,ve,coord,npoin,& + time,icase) + implicit none + integer, intent(in) :: npoin, icase + real(kind=4),intent(in) :: time + real(kind=4),dimension(npoin,2),intent(in)::coord + real(kind=4),dimension(npoin),intent(inout)::qe,ue,ve + + ! local constants + real(kind=4)::w,visc,sigma0,sigma,rc,a,den,timec,x,y,r,xx,yy + real(kind=4)::xmin,xmax,ymin,ymax,xl,yl,xm,ym,xc,yc + integer :: ip + + ! initialize + qe = 0.; ue = 0.; ve = 0.; + + !set local constants + w = 1.0; visc = 0.0; + xmin = minval(coord(:,1)); + xmax = maxval(coord(:,1)); + ymin = minval(coord(:,2)); + ymax = maxval(coord(:,2)); + xl = xmax - xmin; + yl = ymax - ymin; + xm = 0.5*(xmax+xmin) + ym = 0.5*(ymax+ymin) + xc = xmin+0.25*xl + yc = ymin+0.5*yl; + sigma0 = 0.125 + rc =0.25 + sigma=sqrt(sigma0**2+2*visc*time) + a = 1./(1.+2.*visc*time/(sigma0**2)) + den = 2.*(sigma0**2+2*visc*time) + timec=time-floor(time) + + ! generate grid points + do ip=1,npoin + x = coord(ip,1) + y = coord(ip,2) + r = sqrt((x-xc)**2 + (y-yc)**2) + xx = x-xc*cos(time)-yc*sin(time) + yy = y+xc*sin(time)-yc*cos(time) + qe(ip) = a*exp(-(xx**2 + yy**2)/den) + select case (icase) + + case (1) ! gaussian in CCW direction + ue(ip)=w*(y-ym) + ve(ip)=-w*(x-xm) + + case (2) ! gaussian along x + ue(ip)=w*xl + ve(ip)=0. + qe(ip)=a*exp(-(x**2 + y**2)/den) + + case (3) ! gaussian along y + ue(ip) = 0. + ve(ip)=w*yl + qe(ip) = a*exp(-(x**2 +y**2)/den) + + case (4) ! gaussian along x-y diagonal + ue(ip)=w*xl + ve(ip)=w*yl + qe(ip)=a*exp(-(x**2+y**2)/den) + + case (5) !square in ccw diretion + qe(ip) = 0 + if(r .le. rc) then + qe(ip) = 1. + endif + ue(ip) = w*(y-ym) + ve(ip) = -w*(x-xm) + + case (6) ! square along x + qe(ip) = 0. + if (r .le. rc) then + qe(ip) = 1. + endif + ue(ip) = w*xl + ve(ip) = 0. + + end select + + end do + + + + end subroutine exact_solution + + +end module problem_setup_2d_dg diff --git a/myf90/ESBGK/Stu_DG2D/stuDG.f90 b/myf90/ESBGK/Stu_DG2D/stuDG.f90 new file mode 100644 index 0000000..124079d --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/stuDG.f90 @@ -0,0 +1,378 @@ +program stuDG + use lgl_interp + use meshtools_2d_dg + use problem_setup_2d_dg + use lin_alg + use dg_solvers + + + implicit none + ! Basic input variables + integer, parameter :: nel = 100 + integer, parameter :: nop = 2 + logical, parameter :: exact_integration = .false. + integer, parameter :: ntime = 1000 + real(kind=8),parameter :: time_final = 1.0e-0 + integer, parameter :: icase = 2 + integer, parameter :: kstages = 3 + + integer, parameter:: solver = 1 + ! solver = 1 --> use FORTRAN 90 solver + ! solver = 2 --> use C/C++ solver + ! sovler = 3 --> use C/C++/CUDA solver + + ! variables derived from basic input + + real(kind=4) :: dt = time_final/dble(ntime) + real(kind=4) :: c + real(kind=4) :: pi + integer :: noq + integer :: nelx, nely + integer :: ngl, nq + integer :: nelem, npoin, nboun, nside + integer :: array_stat + real(kind=4), allocatable, dimension(:) :: xgl, wgl + real(kind=4), allocatable, dimension(:) :: xnq, wnq + real(kind=4), allocatable, dimension(:,:) :: psi, dpsi + real(kind=4),allocatable, dimension(:,:) :: coord + integer, allocatable, dimension(:,:,:):: intma + integer, allocatable, dimension(:,:) :: bsido + integer, allocatable, dimension(:) :: iperiodic + real(kind=4),allocatable,dimension(:,:,:) :: ksi_x, ksi_y + real(kind=4),allocatable,dimension(:,:,:) :: eta_x,eta_y,jac + integer, allocatable, dimension(:,:) :: iside, jeside,psideh + integer, allocatable, dimension(:,:,:) :: imapl,imapr + real(kind=4), allocatable, dimension(:,:) :: nx,ny,jac_side + + ! variables for conversion from tensor-product form + ! to non-tensor-product form + + real(kind=4),allocatable,dimension(:,:) :: psi2d,psi2d_ksi,psi2d_eta + real(kind=4),allocatable,dimension(:,:) :: ksi2d_x,ksi2d_y,eta2d_x,eta2d_y + real(kind=4),allocatable,dimension(:,:) :: jac2d + integer,allocatable,dimension(:,:) :: intma2d + integer :: npts, nqs + + real(kind=4),allocatable,dimension(:,:,:) :: Mmatrix,Mmatrix_inv + real(kind=4),allocatable,dimension(:,:):: Mtemp,Mtemp_inv + real(kind=4),allocatable,dimension(:) :: rtemp + ! solution and problem parameter variables + real(kind=4),allocatable,dimension(:)::qa,ua,va + real(kind=4) :: time + + ! re-shaped arrays for solver + real(kind=4),allocatable,dimension(:,:) :: q0,qe,ue,ve + + + ! aux variables for stuDG + integer :: i, j, e, ip,ie + real(kind=4) :: top, bot, l2_norm; + +real(kind=4) :: time_start, time_end, time_per_timestep + +character(len=36) :: fmt_string + + ! Begin execution section + + + write (*,*) 'Beginning simulation' + write (*,*) 'Number of elements = ', nel*nel + write (*,*) 'Order of interpolation = ', nop + write (*,*) 'Exact integration = ', exact_integration + write (*,*) 'Time final = ', time_final + write (*,*) 'Time step size = ', dt + write (*,*) 'Using ',kstages,' stage RK integration' + + pi = 4.0*atan(1.0) + + select case (icase) + case (1) + c = 2. + case (2) + c = 1. + case (3) + c = 1. + case (4) + c = 1. + case (5) + c = 2. + case (6) + c = 1. + + end select + + if (exact_integration) then + noq = nop + 1 + else + noq = nop + endif + + !ntime = ceiling(time_final/dt) + nelx = nel + nely = nel + nelem = nelx*nely + ngl = nop + 1 + nq = noq + 1 + + npoin = (nop*nelx+1)*(nop*nely+1) + nboun = 2*nelx + 2*nely + nside = 2*nelem + nelx + nely + + ! get the interpolation points and weights + allocate(xgl(ngl),wgl(ngl),stat=array_stat) + if (array_stat .ne. 0) then + write (*,*) 'Failure to allocate xgl and wgl arrays' + endif + +write (*,*) 'Computing LGL points and generating Lagrange basis' + + !Compute LGL Points and weights (for interpolation) + call legendre_gauss_lobatto(xgl,wgl,ngl) + + + allocate(xnq(nq),wnq(nq),psi(ngl,nq),dpsi(ngl,nq),stat=array_stat) + if(array_stat .ne. 0) then + write (*,*) 'Failure to allocate xnq, wnq, psi, dpsi' + endif + + !Compute Legendre Cardinal functions and derivatives + !(for mapping computational domain and integrating) + call lagrange_basis(psi,dpsi,xnq,wnq,ngl,nq,xgl) + + allocate(coord(npoin,2),intma(nelem,ngl,ngl), & + bsido(nboun,4), iperiodic(npoin),stat=array_stat) + if(array_stat .ne. 0) then + write (*,*) 'Failure to allocate coord,intma,bsido,iperiodic' + endif + + + +write (*,*) 'Generating the mesh' + ! create the grid + call create_grid_2d_dg(coord,intma,bsido, & + iperiodic, xgl, npoin,nelem,nboun, & + nelx, nely, ngl) + + allocate(ksi_x(nelem,nq,nq),ksi_y(nelem,nq,nq), & + eta_x(nelem,nq,nq),eta_y(nelem,nq,nq), & + jac(nelem,nq,nq),stat=array_stat) + + +write (*,*) 'Computing Metric terms' + ! compute metric terms + call metrics_2d(ksi_x,ksi_y,eta_x,eta_y,jac, & + coord, intma, psi, dpsi, & + wnq, npoin,nelem, ngl, nq) + + allocate(iside(nside,4),jeside(nelem,4),stat=array_stat) + if(array_stat .ne. 0) then + write (*,*) 'Failure to allocate isde and jeside' + endif + +write (*,*) 'Computing side data' + ! create side data + call create_side_2d(iside,jeside,intma,bsido,npoin,nelem, & + nboun, nside, ngl) + + allocate(psideh(nside,4),imapl(4,2,ngl),imapr(4,2,ngl),& + stat=array_stat) + if(array_stat .ne. 0) then + write (*,*) 'Failure to allocate psideh,imapl,imapr' + endif + + ! create side data for DG + call create_side_2d_dg(psideh,imapl,imapr,iside,intma,nside, & + nelem,ngl) + +!!$write(*,*) 'after create_side2d_dg, psideh = ' +!!$write(*,'(1X,4I3)') ((psideh(ie,i),i=1,4),ie=1,nside) + + + + + + + + + allocate(nx(nside,nq),ny(nside,nq),jac_side(nside,nq),& + stat=array_stat) + if(array_stat .ne. 0) then + write (*,*) 'Failure to allocate nx,ny,jac_side' + endif + + ! compute normals for flux calculation + call compute_normals_2d(nx,ny,jac_side,psideh,intma,coord, & + nside,ngl,nq,wnq,psi,dpsi,nelem,npoin) + +write (*,*) 'Set-up for periodic boundary conditions' + ! periodicity + call create_periodicity_2d(psideh,iside,coord,nside,nboun,npoin) + +!!$write(*,*) 'after create_periodicity_2d, psideh = ' +!!$write(*,'(1X,4I3)') ((psideh(ie,i),i=1,4),ie=1,nside) + + npts = ngl*ngl + nqs = nq*nq + allocate(psi2d(npts,nqs),psi2d_ksi(npts,nqs),psi2d_eta(npts,nqs), & + ksi2d_x(nelem,nqs),ksi2d_y(nelem,nqs),eta2d_x(nelem,nqs), & + eta2d_y(nelem,nqs),jac2d(nelem,nqs),intma2d(nelem,npts),& + stat=array_stat) + if(array_stat .ne. 0) then + write (*,*) 'Failure to allocate non-tensor-product arrays' + endif + +write (*,*) 'Re-shaping arrays to Non-Tensor Product form' + ! for solvers that want the non-tensor-product form + call reshape_1d_to_2d(psi2d,psi2d_ksi, psi2d_eta, & + ksi2d_x,ksi2d_y,eta2d_x,eta2d_y, & + jac2d,intma2d,npts,nqs,psi,dpsi, & + ksi_x,ksi_y,eta_x,eta_y,jac,intma,& + nelem,ngl,nq) + + allocate(qa(npoin),ua(npoin),va(npoin), & + stat=array_stat) + if(array_stat .ne. 0) then + write (*,*) 'Failure to allocate solution/parameter arrays' + endif + + time = 0. + call exact_solution(qa,ua,va,coord,npoin,time,icase) + + + + + + + allocate(Mmatrix(nelem,npts,npts),stat=array_stat) + if(array_stat .ne. 0) then + write (*,*) 'Failure to allocate Mass Matrix' + endif +write(*,*) 'Computing the Mass Matrix' + ! generate the mass matrix + call create_Mmatrix2d_dg_NonTensorProduct(Mmatrix,jac2d, & + psi2d,nelem,npts,nqs) + + allocate(Mmatrix_inv(nelem,npts,npts),Mtemp(npts,npts), & + Mtemp_inv(npts,npts),rtemp(npts),stat=array_stat) + + ! invert the mass matrix + Mmatrix_inv = 0. + Mtemp = 0. + Mtemp_inv = 0. + rtemp = 0. +write (*,*) 'Inverting the Mass Matrix' +! have suspicion for the matrix invert routine, but it +! seems to be doing the right thing + do e = 1,nelem + Mtemp(:,:) = Mmatrix(e,:,:) + call diag_matrix_invert(Mtemp,Mtemp_inv,npts) + !call matrix_invert(Mtemp,Mtemp_inv,npts) + Mmatrix_inv(e,:,:) = Mtemp_inv + end do + + allocate(qe(nelem,npts),ue(nelem,npts),ve(nelem,npts),& + q0(nelem,npts),stat=array_stat) + q0 = 0.; qe = 0.; ue = 0.; ve = 0.; + ! re-shape arrays here... + do e = 1,nelem + do i = 1,npts + ip = intma2d(e,i) + qe(e,i) = qa(ip) + ue(e,i) = ua(ip) + ve(e,i) = va(ip) + end do + end do + +!!$write(fmt_string,*) '(1X,',npts,'ES12.3)' +!!$write (*,*) 'qe:' +!!$write(*,fmt_string) ((qe(i,ie),ie=1,npts),i=1,nelem) + + + + ! ----- Ready to send data to the solver --------------------- + + q0 = qe; + + write (*,*) 'Set-up complete, commencing time-stepping' + + call cpu_time(time_start) + + select case (solver) + + case (1) + + call dg_advection_2d_non_tensor_product(ntime,& + dt,kstages,q0,ue,ve,psi,ksi2d_x,ksi2d_y,& + eta2d_x,eta2d_y,jac2d,psi2d,psi2d_ksi,& + psi2d_eta,intma2d,psideh,jac_side,imapl,imapr,& + nx,ny,nelem,npts,nqs,ngl,nq,Mmatrix_inv,nside) + + case (2) + call dg_advection_2d_ntp_c(ntime,dt,kstages,q0,ue,ve,psi, & + ksi2d_x,ksi2d_y,eta2d_x,eta2d_y,jac2d,psi2d,psi2d_ksi, & + psi2d_eta,intma2d,psideh,jac_side,imapl,imapr, & + nx,ny,nelem,npts,nqs,ngl,nq,Mmatrix_inv,nside) + + case (3) + call dg_advection_2d_ntp_cuda(ntime,dt,kstages,q0,ue,ve,psi, & + ksi2d_x,ksi2d_y,eta2d_x,eta2d_y,jac2d,psi2d,psi2d_ksi, & + psi2d_eta,intma2d,psideh,jac_side,imapl,imapr, & + nx,ny,nelem,npts,nqs,ngl,nq,Mmatrix_inv,nside) + + + end select + + call cpu_time(time_end) +!!$write(fmt_string,*) '(1X,',npts,'ES12.3)' +!!$write(*,*) 'q0 after return = ' +!!$write(*,fmt_string) ((q0(i,ie),ie=1,npts),i=1,nelem) + + ! ------ Solver returned, upack the data and check the result --- + time = dt*ntime ! update the time + + call exact_solution(qa,ua,va,coord,npoin,time,icase) + + ! -- reshape arrays --- + do e = 1,nelem + do i = 1,npts + ip = intma2d(e,i) + qe(e,i)=qa(ip) + ue(e,i)=ua(ip) + ve(e,i)=va(ip) + end do!i + end do !e + + ! --- compute norm + top = 0.; bot = 0.; + do e = 1,nelem + do i = 1,npts + top = top + (q0(e,i)-qe(e,i))**2 + bot = bot + qe(e,i)**2 + end do !i + end do !e + + l2_norm = sqrt(top/bot) + +time_per_timestep = (time_end - time_start)/ntime +write(*,*) 'Time per timestep = ', time_per_timestep + +write (*,*) 'l2_norm = ', l2_norm + + ! --- all done + + + ! deallocate dynamic arrays + deallocate(xgl,wgl,xnq,wnq,psi,dpsi,coord,intma, & + bsido, iperiodic,ksi_x,ksi_y, eta_x,eta_y, & + jac,iside,jeside, psideh, imapl, imapr, & + nx,ny,jac_side,psi2d,psi2d_ksi,psi2d_eta, & + ksi2d_x,ksi2d_y,eta2d_x,eta2d_y, & + jac2d,intma2d,qa,ua,va,q0,qe,ue,ve, & + Mmatrix,Mmatrix_inv,Mtemp,Mtemp_inv, & + rtemp,stat=array_stat) + if(array_stat .ne. 0) then + write (*,*) 'Failure to deallocate dynamic arrays' + endif + + +end program stuDG diff --git a/myf90/ESBGK/Stu_DG2D/stuDG_F90 b/myf90/ESBGK/Stu_DG2D/stuDG_F90 new file mode 100755 index 0000000000000000000000000000000000000000..89e24abae29d2e3d73042fe63ed25287b0244f87 GIT binary patch literal 94586 zcmeFa3t&{$wKqNk14fIU!A6@Z*65wKi4U46+Ki7jgA+Xm&uG+GVlOqJ0~A^m3KONb zNTbOpV;q`RRBk_SOYgP!($-#RYkO0wIN_Cm4+4tdqZ;35KqZPoM3nFMTYI0G6UYd( z{_o?v|50+zeyzRs+H0@1_S%p0gBttdA%4Hl{SEPb#m7*6uv((jU!Sk>%6y(`UzP7j zU$yUZz7u`Jk!u+Kx%sN~oBQn7#oT8R^71zn{{r||HH0BPtNF-1M@h)%bDwVQe6qg$ zuUfyV=Y1EcobJ=-E0uJ{64$LUmA*ALQR3Za9TCE>nMld6RQ?h@mh0E$`gJ+?>C^bT zPq%#p!+&);y>3V~LJU{);fu=O{L`(r_SM8BzK5z5o}sFU@rbrRmeYdG44mO`mk;^l8^O%{g;U)wyS$J7)Ck8KcjVb_an7 z^Ime<)q;{s#}b5yb20vr#%%{^g?)k1U^m-JE|o>8p>H{$>vI{~i7jhH6*P zM>b^~e`PhB?c)}WxJZQ#{QHjpaO{E6zF2|yrwYVRL!E=c`6*xz7Jp6wIHX7EAo!eI zAbuL+2ZKMOK>VTtcz##_p8$wI82r8h`OhwZ=kW#LTwWmm%mVQLy+C|tfp$GqfDVro zfKyx`{?-EV?FHaJT_FF?0(ibqApVa9;Qw8L{Id(hI|cB$qd@+i0{K6s^Yo_eb!ff#YBk=A}l|Q|JjL2@zYP|m_OIo(D1G6W?bJeJ27!qqM^apFz$wiD<^+z z+U&&SSr<*8ID7Ww*}ev=F4iz<@+=))7rSWsjO!=YP5kEc$trE|l&@bu@j73_*Dt&D zs>C_*#P?=Sw&FFFXE($eZk#r0^7ReVuAi2uxq%t3nsLsQ=}ohzO6aOY<+%+FiK(+@ z+}QBViIW=Q6K5xwYtqD8K2-9I*@;<`Cte2_vnS6=SaHNw#-~o4)sUDqaav+FDyq2> zac6&X;_PYh2A$igYpAPesJZIghT2Q7I=f-IRGp}rHe)s_nmu`p1%#rN8l^;7l||~H z;LnLqpFD9^4JtXS!2$&QOH53F9>8!+d`443pwE~-yJnht3!~v0-ATl5MjeRFwx2z- zDNz%b_x+H$h{&M3zBOgWEHvzTcG%gz1~Ca!Cm z-Y{#7CoWfF`iyT)iw{yLmtorViOJsr=>4J2PN=AcNnnA&AhR1L#=$)Y7M+}Em@zXk zZN~LrsT`UvKu=VBW>Z6A+I5o~yr|Dc>4~$xHJb#vPOH5_$&-DT*w)yK8qOMh*64H8 zdwqSw*`vpd4i5-bR(gZJOD?|r%DO8fmo-$59^(yEV6Y8=rVPO98K#YlVbF-fd`HNC zKKx}4<|#t1p?DjXdk#a~5fWx;{`~mwdB^})A2ppY?MT8n5prCm4;7uu^lJP|{h%DV zi+o>1I^!nx--GaQ-#I#*DE~1`iDABIuDsvZc>w%?2mTP>GRP3Z;;(h``^Y`iSF7PQ zOL&Oy8ya5!9B;-`AzM|k3Adg71v#Lw}> zALEH{_Qe04Cw{&s-t99cw|nBfcElo2{HHx|7JK4p`?;Py9)q{C%GI&w1keJ@K^t zazCwK$s?t?NQ^vB{3)LJpeO!RPrTXkMv>Vz?4+@$_{?M`G0gcs63v#?MXOY%edipD zqN#m5JHjefWL;D26$7vON+W!Z3$SHY4FZ5rjEZvdbj=4#ONG*~Jpx z#xRFQwq3$6G0Y*6ZI50ZJgYS1fOq8JM%wWcjZK%b0u;uF9jQ} zVDMoYj?A6?HsYtWSH5YbHePF{UNT!=>8q`?x-zwAMSNygXQq;mE?=#g+ShBG+|0NRhTQy}{X>!8OszCK zw|&+0_n4`*=EjejP<}b{n#oPM@-OM~%@* zKNt#Pc$i(Sp%OeNfed&p#p`9j39n3TR6Hmy@|xG^3nz?cRK}^CM4S z5a^tN6@Wh83!SXGVj-dXJOGyi6W@ya8U0yLZfBGy&mfi77OwOqPs!!v?-k1kz3qd1 zlT@XTGWtGGZs*@Tc?PNUio%s%=K-VW*hA>YdjP6RUt#nNPj2TdPo6<4y|Qqn{hpG` zDT!7%g#PnAo|bZ;t>|a;?VjAuL!LZ?RC-n6N?SZ&mQ#GK@SP3lXLtarmX=|_ILCQ% zJEkYkAeGh?uJjlWnB`>L6(d>cf4uK$sW4IrqkBB{EEl#L#pskLxAPzZ4^x_hb!ih$ zQm>8q*J`{5zbF}I;uZfh@8bMxGhR3Tlf+iwwH<;C88UcvFyfqFvASf&Od$?x%~+T* z{+L1F6$V<%zUuR}EesKIsx!ICu+~>*lAGooVa7Yn&fXy?_%Pu4OrxgXY0N|FN1V~R z{VN|k-hZG_9c;X>@z7xTS;)7dk)*l4KTk;K=M1zx5Mo_s+k=9=nO3+u(QAPB9YTi6 zP04kJwGnu)Bi;~(`{BKWfD{GF&={lj{tKniq+#W{gEaa(`NBm~`1>eq##@DCW_+PS zVVxO&8#NeKUvd-vt_KN=k{RQNfv{*5FBMP71MCwjo_06@S5+MVaNX`)^MZoms~3=9 zs=G}PiFn&-MJyAYO(Kx6%x$N2OF&?vQ&ahhnF@rQ;%nK;g&{JR6@OHiKs60Yvh&5E ze80StR3&FzoJXNVP*`<7i$4&mMsZ+5_Ojx8le<8hk803UnTYXI`D?B=sqOiyx%>S})U@<{Gu&_7yV6W`n}1R{_I;J88cVy(aCgMGXLZE?WKA;D zRFvFh{P5quOd{1JUon#B0g)9ziOSSW?KV>ngtjt7tlM|x2r5v|SCP_cTd{)!>$vyw(T*trHIPj7eO!SvDjF(biCT#n5h6FS)#Ng^fWwpwz>GmVH7Qi&{5?z&4<{@}N%UxZ z&jGw(m@8dY=!LR5ELPBbVBFJb{BfOseV!4{AtPXfC+{XmJq%L|^E`0JyEz`H1>$Ir zVu2~lpJ!mAZO2z8ciD!uwsK!GQ*HdT3j;il^9==YwvDBob~v*?`&!1*40zYYIk`{q zZXk#9=mID|y)%b$y%*(=U>6(~$``Z1pd;?^Q2zQDm&zZ(9vY0wU&y08)r<1u%zrpk zzWifRwm{{d4HOausJ!fK7v&=eaS+Cs{Z5X`jb4(+`O1 zqky2`pje=wpa2mpAQDJ8Gc`9?BoK#Bp(_a7A5WTq~F1bQ&Uf$OaLAyT2RGQgp7of7%5VN$zOcZ7PRq>}ui z#8yHFK^g~=n`{uR8_AY&S;iIVbJ%8OFU$dGU?>296lDV&$Pz~D4Ddaxl;%ra3b=yo zh4_U$#QA_|w2om>P3=MY>3+Z;T+lmho(@?q0{(pV)0iEee)<*y>=)Hf$`+)Et|CoA zcvssI7a9QUHxiPqW$y%Yor(@9vQwdu6SzsWmE`RJcD}qtK!XdER8vmlvYckDBua{m z&C=-AXDJv&OJGF)bReW*$xWbK&a!yptsIa=0^|q#0~v@2iXVW;l01xM0^<_G$k9hx zDWnTi*}aNB4~B>~Nx?>&pi7aN8Uv2ILkLACB}&fsZe+7Gl@DeRT3@2zaRrWy-6Qkp zbqIQ=R2~Mvq5}bTyqTj}MgaCy90tIy(Fb7otvtXU0oX=>U|%v$j#prs3!%VJ?BP5i zVL_l2hi(Y9ht`I|L#u+#)Fqc=#^qbFm_zI7?fGbD#2W~JRZ55>@Sek`%Ob$XWf@W?mg;rx z#lkivB7!Kz)>f`dJr31xEZvmkG^`y_GHrlWvP|AFc`)O<%=kJrg}CI*1NvhNF;XL| z7xf#hw_XA)z=8xui9|q^h2|u~4+MB?9$=A73!RG(1K^r74#4o7Jiwp;{1*+d`cN=z zJrLk4+j33DiVvDxRiiLGB!Cwk2=LB4z&@ejPY(y+>!%&i`0NoFp6tK4%zwx%(+KyX8EJ zOY4e7k36uhG5NYmq`x{Y0>HwL0uwuy`JN(YSdt~o`o)ifRe%obBWaF3!75-z4AjZd}6+E ztrT93!X$SI+XECB%P#9-mx02pg6$=qIb)?O5Te@@HWIV(=ez@x!ktTShKb1EU~I2J{77CO6``MoX(Oo-Om;cgqoWGV4Tb#pq0iCx|u!D3AF%<05b zVASr7cKNy_FoD8DSk3f)^aaWgPEOw+i#%Q*(XLl>-P4TxM(Y;jh_>A)4S*xXOg$Ld ziUd9vlWv6v)*jHNfGJu1T7Cp65B2BI!NAJCl8<&q5fk(ZIgGs_EQsEbfrefodhP&Y z?~uaI_*!I+t4uI<;yv1ZvI{O5Oq|!^Aumq9ij`3fQ_nf*R({Y2AB;CXk@;XNeHKBUC0`gb4z1 z+r5dcS%N5x&dfXvXq^B{0B$9?S~t0_);-R9FH3b~-Ki*Hv`+KX#qu_2rbTsk9{QpI z&?R9o==bEIHwyId*xc}dLlnUK469qvghO@X0oaYtm!B`?O_Yza6J*T7!-!ZLDq@*D z9}IQiX~jFEAo>zw3qq92vo%zYB(vo~3Hyu%YZ=4I)+3uOc_pRu5ee^0IkuRLok_%{|!Ul!;kryq`_|0_fJ~O<2)*EK{kwk-)de$m` z#7@0ncJ3Hrrvnc$TMR3K=Zs6&jr3dT5eq9fuES$(!Y6AQ0O zl(;LH-g?5TZtAmII{mhPm*w9U3ulewg%~Vmiz0+y=nl%*L)Abr{jG{Bpucoi7;*sE z7|ARaGD%~S+{UmH`a(L8zU&hA4m*|o4q}}pou`=J9;(>Brt(=1$f@j!HjEu+N!_P= z;eM%4_pc;D_7T{)^q6~W^fl2iJ6}fkd>Di1Snvym(_hm)@Dh4p8+zcGmY0Uw;XRgd z$rh_+-*C&gcnjA=M;RykymU#Mr2ma?Sm7;+hFI!-t9+xK+AtXXj|KhhyD)&U_}Pjo z*~&ra-#IA#t#B9Vzshg=pBMVSY$U%5`je2o0A-_%!h{G5xj{er%h}P2Jn0Y zTYA9-4!B@zOYcxB4d%gM%!?dGa9HXXxN}>7cOd~hAOI>ibp_~g2Ytp3A>L-V=SB=t zl|3e=iTQQ}=4dZWUQWk}YAS zfI;zKrUSPTF&PZdnTnjPp+!0mtCu`Tll&-Os?(&YfFq3n4yrNti+yF6=bDaJ6Qx~- zTMUSx(~t1i#^__TeY#75o22c%jKhfJScVRvJ`_B_Uc`^XoltCqiHq7#VJ6?aG+n$% zG|9ZrOljX~1VOAlJB^WnTSZxf`h) zgbMJ#*NeYG94#1P`k#SCrtmlYk6ZrjP(n?Q5(1$QP()?y9*nhKAYdd@TH%mEn4O)3 z`YRZOJ;rBc2F+apFi|oGp0xkX9%u(kg{kDpa9ADE#u&T7Tx2YmjqOAbB7}lIos1AH zNyVMC0bF}fymf#7j96w?s>^6yfYn_ZfMN^-4HDUPyTtd{aR(+E7ldJfRIw^>fN&ia zh%u990%oX4P_yWNKaq_XP}y5gOJ!s{lGted7pYNqXRhY;Hf&CGt^-)K->?`>E2p1z zY79^_03ejFnax!-TL;#h$k$BXHd^l}(AGXyw}&n9q4lVlt?e67GrMp=&16zl^D8GD z(B?f`HG$U}tzW{5hmDB@2}>hEYS_-Ub+ZQa8K3C0lL3+JO2|}@tWp%DVjvE~#Yln} zA!mqTp*WdiM_{3!FHMiRn^K5eZKc*=D1MxNG{f7VO*We?A5w>GE~Fnf(yrT_7;mMr z2hoprAEtgZ!)ug1Vf){*{LjR~uqTv$R84RN3<$?YlUn3LQ0o)~Ni7k-pcco0rWT_> zB@&B+-a{>>aOepqD{8?MyH8VVHROk};HLmTxHVygH^Yk9Y=K&Ql3FLxem=mS0KM$+ zHrf*hqt}7<#FiVC{o7r@{>stoHPFinZ#R+`YkPuS4mwG*(MFeE98jRwM5K~l9INCm z4q|P7aad9tF^UqNfIz=81WCU+c!Jx+>H>Tj{SZ_AOh!F~PpK_XKvU9_b(0H;3#Xej zZ~wI*UCrPui5?_eZL~~Rb-)^%SMqe*+|es$b@K;iTp{q;h{Grl+YWw zHN8dulHLS#aC##b?E(=^sE6=LZ-Mf^oZg(cz9jTnD6U)+jNK0Uv}~7&pc?L4VZksU z!QCl7Egohin=ZlAEbpNuR?pY;gcF2~mob3!4n8|y8)CXrg=JZ-Ow;cZ_Z~T;7WwS- zm3|YtPO*iW{97>pI8>^UPW!Z4%F-mn`ZLQw-H@9iw3Zbf9A#E&mxVbTfC+DbV=55q zW*`hCEkYEUMv(9VE%*sX8FH%Xv6F=*FziJdARkAm?eDoP9k>UU{w>U^fGg(#lFam& zn;A{O5*%^Mauo9A$CDNJVI~ptZ#7{{!?=qkeXxv^Q|pN)MYu@7K1uny{LfBCuG_fv3kLd%zvt>?rm?Wrm&|D?Mf@q@tC+;HNYL zNDR?D^v4`b^8{SY^M!nWpbG-G(z?D)xnh(nAPaSsxG|P`%PN106I=LmE=!L}qYl*e zY+{_5woilKrtSoPEM5FV+Rfx4cx!%Cxrt_QHT;@ofn(v7@kfkL2Z9mfGtpoXXMT2g z-K>{w|JGP|XX1F+%LD7Tpwi#8T2}kbq2O*3aH#4+5>5!|672l@qyvrQM`9Vf681jD zyoKGtWdK@W!dTgppBXNPQ{K0to)gAr9?wTRM{iI8oVY@VIRCEG7Sk|;_jj2*v6_n~ zy}~mktj$%|Ym5QA7|Nnp0GUUhGsH-nWXF0okcp($dn!(kd>Nd9jIkW2{{}9WThsk1Jz6^(g9hi>D6T_ zn6gy|WTmEG2TGaZ=s*breL4V_FH{j(8w^7gL=QrRPLwsBK16a&=iydTe&z~<(3QQQ zB8VM%Y~Mbx%dqdghiE`3LCKxQHN&n5o`Ueu`*+zMQN3@_c&T*Hm{#O(nK&U!;k=_K-3-Kk<>dJ*i51#Hu+3{ z-H#~BBdhaOneP*7H3h7sjJ(kZrNfHj|a* zPgbG^hEynGRgOt0N`V8C5+^LsEWq9(fut(SYX%l$SFjT?sG%8t=4u92xvLqrx|)H2 zk<(0T&EV<-NA!WK0*2fiFj2mT08*KHk$hi_c@N5Qa1wPgInAMIYh)0of+o?U18yQASA zyk&M8^-tqpPZMTRYyEck`KF~t*~*>9wX1d-S9BYr*BQ6ZFC(J|p+d-TTDAcaHtIkO zfhe#}m;$-B--(mSKxJ@ghYE&Oyd$Ntn2XZRS&!u9oH{x&pM&t(b&P&M^0;-&!2+jq zeZGM6u4KdJ7GN|8BHt>Ah_l=SX8BPd-HLlzZ(P?aEz@mQ^;(_+&Rk>zi*qwS#%}US z928>@W)ODz2m%VZA3!mK@6c%sI#24fGU=?=cAdtc)1lLfRNB}7PQhf*xmTxAQUPYQ zF2|tLtkZ(3+?%JWatu0?bQ*;g%C+k>1|3tU`Bd5$bsB@t89J>>rTy*{1(QK1sMAc9 zR;kk%bl&Y&ZLd{ny``!ggU(i+Rn4)ataLN;7pDgU%S8=2L0+lqg&nbV_tupQ`UQ94;iD z3_5$)sI*mc6PN&tXv?FvHgH9FF#Ad&5B(>9S=Zg}774j=3VxM!o zE?uFp9gb+Wl0j$B8chLIids8^DvkDr^Sn-LR<*uns5Azhl}L->$VBT%#M-IXo!{!r zQ&r~kai)}&FzDQ+(;8J;bSWSM@?c*=Ho-ld^B&r}bb zX3u-QCx5Nyz0dRR^W>K?uVMCh@(2&=_?e#gIi7dZ^Dca;^ZPvU%RKLkJnxO3_o<%u ze$P7x41c6y(DPp6d9QNcT^*A+nhYZPC9%rGg(_Y+aDtAP{#Eh9^HCjNL6|C@ug)1d zUiw|d^VKQV@$*Sf70*{^?@Ds;CG1G2mz0TB}ubHpT z?{s`oz6QQJKhg2RO$sPqom+H#ZypX`op0-S;WCw9wB?mLK9kSSSLZw(-;s~!t8=Q3 zZ_MMvSLY}lKQ$lESLdBh8P9Cw>&~kPkP^-LOnh}7)A9BB%J}O1LC4qTUdusK3|>dbiA;-?k48HTF3X~;qcYDP{(i0$Me&lA2|)0{nX}xVSB%Dnqkh_ z(D%7Quk4rUEh;Y$&@WM)v$jiV)5r@gT)w>LzWj{FIb*Cx-6Xzvp!|~F3r~34e&N$M z4ov?w99Q1_Z|oO7^NE4!*Y6iz$}pN*>TT**xl4?Zq_0b~`*dV4x*!*M)Be?*I4|;G} zj*aO!B+kvflOOb;xQGFir=s95`xrS25gC_Kg{jkx`~|9xc&k1Y`Buzf)!Y8;srn>c zb&swoJDd(M09JZ$&Mup5d+fpBX?+lsEik$HD*eF+~0)@)T(4 zcyD!>&G}Z8usWZ&y0dk4^WEyc!665LKDVWK7cqkJAA<&WL=-@y$5Y?(FA7w!GY84@T$f}W;GuR6E`)!5Mu)Jt|k!UO;8t>p>vD+@^^BN72v&fRY zpIWwXb#n_>_l&1`s;{e9U5mH6KQe%Zxx+pCuaxKe3+PqMQs(InmnJZ)y`}l?+R`w$ z8FYsZ-zz*2a)TG~)$TyXd$T7`ktd#(4u6Ex>&cVr5ya<4CEhzcc^W-sX`%4PGu*7s zlT`01QtNr|^St{!`R9A$dpz%3J?}F;?{hrwrsw^?+KSF95sOoLXq*{VTO5q!6|X_o z{sMnUr~jFqpFIj|6#s+Km+NEx=TiUw>@7@v_KGWX`0k1gz8_3xXsP3h4J~OtX~~gw zJ~%Y=vnGAu1-fbfcDilLp8i)9A~08QX5ar9-Y!(D7XaT~5@!E=7{FW%1_1bn!vK)0 zwdN2DWk(%iMxyjd^8sQ0SEz>U8}A=FTm5%#^ucsscIf|-hLV$f9>u@R^S;RQ-spLs z>UnoO??KOdiRZmazh^JvwF;1Zj*Zo_2 z&`K`(^w+lhu#vj%5Rp2p@!(~SLVa~ur0d#q=v;cZm75;! z^|{01YoBn1gI4nPVb~gn(yGBKAwx{{wD6!R=ALbzEU612bx=f(`703Fely!!=4o${ z=iO^&xb)1{ABOaG3$YTHRrr3;wv5gGbwoW_CE2X=*8v_3{}lkTUl^p2 z+f&&inaH&oxP$DJWA4Vi95|PSPQ@!MtU|c^>Y1fo(_DK-E8H$X(o7GMi#+dcKj)6uBkvMV zdbJ0}d{2C%=Y6Vv&(4JBIfoCBEJFhE8mvOb4^~kQ5pKdn<;*u&#o3{Ba~EbVjs+AQ zjJaC>C0}wmHe&f!^pLsE{O5eM(};k3X$Dc4<7uzk6O3Q%iEno4@Jp_lAkpoB%1zlc zy2|-ADmKcxv)|c*1w33wk)8VMec05J!CTM2t0FT!^ltRLPu1_)uP@VOzx*3jqSlkE z-t&G?NI^hSl`qzy4y*XarEc-AkHJ0rD`bhpapS!M&&xPoZ)t_jpPho`Y!4c4T^fyE1@&*LE<5yZsv^t1ZtnhE z-rtV1GQ7gwO6`?f+~)(uwi9+njo{J8;#fucJNxYT_MOIQj~Qj1*vDJX4ZOqluEpNj zsBuO=Vouut3ko}%UV~l^L z#wg3OpU$~q4R&|rdS_%blt6!F(ii#7Mc9v??1z>^7vZo<-~!y}IM?>?#;)mTXbz7P zjAPmIXsA)>gL0j?0IQxyA!==8Hc`PyJ5<5Pq)?ELGeZ2E2$ce(4EkTTbPs!LHBr14 zbvsA=l($--tl0)AvYEOKFi|SVQg{f6pn?|!^xSR5pvB){WjY%AN{u1ximDWlxhZ+*l6eviyTB=GP+TS|mJ<`mfmx z9`5GLUyQPsk>Q#hfZLDAUKHDa=j-_Q2ur?%8hZHh2wtAwSz}z&W0bv(hLH@fqQ)Cs`g+Gp@wWsIKy+Si@5 zX%!B*8KZaJcx1~qzmePtV**)UK}{olI7h^;x)WKf7AJuCN$k99SECz`Fk5y6&bLEz zahnLrT$B!cHM|OA3bZ3JL*N`^{F^Qs8sm?fZ2#-gBYVv72IKBdEBtP?apziG0*x!5b{+ZlDczOp zA}#xd;gaQKx6yh!Fu*kzM(e3~F~+YmTegib%8r`~TydbR7!@8*MP9ise?APO3dA># z4IE({A3dVTSerS~7{8K-x{7(ItGEamkUAWx!*RgN2{H$hk7x4na^@=%?!NIqYQh^L z>4-lu6W3uy%Dc??Dx6j!rPg<%_T4YwM2JWXUbZwYeZfPd66_2}%k-ELmj6lH|F-3S ziZY}bJ&TbRj1J#xj28k&i@}7!sDCdI2$15@SU6iN(w1>Q0CZ zY&KkPw&7Ab+@u#YQ;$=ERPaW<5?qwTvp~&21$St_?6klGAi81*7}RV_go?(c1HV8j z&IKVrncfA3PWuG=w5H?gP=G2kTiy!rGNoAR322Y5h)6?Q-A%q`zx$X5yA8+;KZ8U2 zPk113 z#r!YF!jBp03(;V+jjaU{g1UKab)izwgLDB+r-w?JLr8)9bn^F4{cY`)^Wo3v`=_2L zaR1cDo)6Er_aOIANqb)!&|WkNSIbA+szg#0V|X&6qNeiA%4fma=cI5HF4|TRN#B6O z$$?N!au;GlxVqzMVVJs55stF|hTLPNwhOn${Oj$ZxXRcaFF1g7Vx$Ya*8J5-7vM9e zqm(J_j_~^!H$h+G+ZYpOYMW8EUAP|-58`Zs!^9zwZIp~_IKqY)WoyZAVM@EB5+~py zpW`v4kar`pL84M6;Pu*-G-5Co!+lMNlEU>o4mHvbxni|s0^zrtZfWVl>#6{b_~Kt@ zAiUd1UP?}&B=o-;3BO<@he>#~x!*W$JJJ}xj~eypDGxy^Krb7~b7pTIv4~JbM+V8ZCsL?0G-Oz^b0Lg25m;r(9qWBPa4vGTW ztXnileK+zzw~F2{u7ysy770*F^{WwA4|}j|GnMR4{@G4=&mrn{>QCbu(uvaRnm(iK zB@l~xG`JKzs7YLC_ znhpxRUZFLIY#Z$RZ2?QJ6Q~1typX)8sy_;uxCDlXm5v66m_yM4k(M>MO(04dYtrM+ zDbSyB9f(?j19J;dQ`jV-;D8rgDbe0JSBqQ)AKIM*`FBw z_MZYUYg)j*_;wD$(VL;7aO6>>MLnbjFFtxF5?v&V+n;2AHu zdi4E&t1(`7TmGGP_*GZm&kglZqOirJL66n7ibPRuWb53`1(Kd1Va%~=wKIBDeEW5F!3ck4p!xG4++p`BM_VPu)JUGyr@ouHwG*i4OA5tsAGt=pinnk6f$3yJah2}%=+`@qfCDzT6 ziOcrjRiihR-T*v6bhU`i6Z0bd${o4owY&$RM$_NV;>$1*1*p5Itju(=dLq7>qu4a# z9N0Km&ZQyH?!17Ro*NFBr{AGsjHPQS!J+sPzz8L3S<_bPS_FAnHiSY@DNd1Dq9>UH zq=M%lfdVk#dd;B;rBUlyK1dOsF^5u{fQ8UP8?qWgZxJzG%zCKDl;2QEL}T zKjMVSO%V8_K5P45P%{z#^B#LMaI9EeJmya01_sZ_z;IYLBa40b8!R4a`(gM8=^IhSyfmbshJnQTA_{WE1+`ap$`K(r@@* z)ZVH>?XCKMvAtLe1+NX>WL#u{w0x_Toh5v|7zuwfU$4J$yIg2jkgw11@U{1T+`o#S zqr%U)+Xbj|1h_$d$93K{aSeZ zQo&;U+FHS?@eA37p*o44@^PYGMSj_^MB{Z(QpVCc$g_0x7^p|wQwf=mY0s?~^^-7R zph79(pQ*@kYmu1>gm?=lKjXMkbfTWx^+FTCkmUpdB};e{4N{bqL&Lvh+cars(Uc~g^fb>45OE)4t$n=r{=}=af4$55O0Q;hahikoV&*gP_;wc92U4rWxibjC!fjOAZ%`FR73 z<^Ra?e`xvNHTm}g)4$X5zhmK6WPe|Neg#DtI0D3`-W2T_3@8`pF#13PL!nfocxARB zl4q#AGiGE!I=B6N^ik6vp)?mEoIw8#h2a2X=R;jCGLn^|Qac0Y`6EI`a*ZfY5)n98 zi~0G1jHL-wi4VwyxY{#0#n^$?f){un-d1GN(?{M|2S26Q4>G>g{A95E;q7o1^ZCd0 zir8ZxKc7ct-cSC9RsJZ>baOrrqwyisi(Ba5%-;`xz`+AO@wnzF& zyqUDY8Zp9X(l%PDN8}bKbbIG3L&Q6GhY*e5txiwBwTO2=@v^@NreyH%&5Pm=c02Wy z$glkTr!u21)N}QF+lb*{j?@U)Pg|pyV($?9`!(CYBN|@Sw2|I4@v32p{VRCF+tCK^ zoVs1jr{QPBC#=oXdsmq$_&&GcTh^SM#^RNcHXnUVN=Zc1bB4k{geB<3fP=~JQ^^dz z@YyeZ8eT8F2+@(weYkXdQ!Ea*HSP{X444Q$s2yw|*ugfUj-A@!9D56=y*o_ok!VIH zXNVgHLpkr=EQ#=OIot3G2)K||^V%Jf-&oKB4DD2Az&+=7xF7t~1uIuyc;AJ|@+#aZ zY_zVy-9nX{U=NW75P-RU^f$Pq$@gKuV2Ew8QyRc)>nH6MKSp-Y{D)571LNh~5C0{znaf4U5Ex^Byk6dsG zGkRw~vcx6J#$gykVR))~Zh8?($6bL~Jcq@YI(n$lG8W|p2K>$2k5hIT;TWR{rlk_e zx;z$CgGY} zpoM)WeRxm_Q?hQP48=9qOzOZtpsxfD{Qxupke{vLFNz7>VmuM`9K-n0{XgtGzo0a* zp~dB6JDoN31!CJ4(n$v3(4I(JU0WcMdOZDIf6Jadkc7QRYU<%V zfP`Yk=cd`)7e0HR^!hn*;BUq2Qdry$PWtm++b^pc#ptm5Np~W+hpilOC)D~iy++v^ zSTM*MbvumF?-{pJ1+^#&`FR^WNfk*Y+V^I!g=z=$h^Gja!fTDP zM{e4dU5yn`HNM;Id-n1wToMec&R^j@co_jW;t+iavxd zHr(aoJ8!Uygm)Ur=YbJKed1IIc=}&Uc%Nl6pg0rea=xhoBF}N<60^NNFxxo?X3Hc1 zr7H*M)1Dk$G-%y};Q$Rsy#>vgS8~}bTnIV-1l?XT>RvZSziZsey>1``7o<~(g?qCV zFnut7WZXc{W455UkGBRc2J+kr^IJS!ej+a9PSEtE;(sA)>Ot3buNbzKysq-N+nIv? zmtId7TO_WqAmRvB^-C;X+%Ya*cNcp)nL)V`cP746xvAZXk64J6X)J;tW5yGIpsxl^ z2>jltTV;&yF>ZU3_`m`P_cT38?|LNse$(%H-65+Zp^0D=2yw8~V6n7y5?V<%X9v?u z(v2ruaOFRPx4Eo%+*Z7S?W_3#3u4NzPcXL&hJFAme!*|CA^<1YLtJ~eTlQj`5SR!u zRQL$2oaB}qEaMz+#%=QYglOAd7wx-{l}yN00pKlc)AZ|XH=4G8~?}X?ylKtn@29u$;#t2(+Sdp_H5d8@-Z$HS&Kxm!I@& z71+7_QFpvGtds&h7;nYERp^9u#rpMJT*OD*%b*4(^eXPWu^Qu@dv5ABM2yjI-*CCx zm%MS(if;nP{%CJ@Uanrv`FXZq1&P2u$CMAN4HOFur4sk-&6+4>rXQ4vAOdl{f`|@r#F7_Td z$j9|ynpfDx*oAG4STwmBmU5&n@?ZchzWS;fTJM;g(&O9_i*Jk~;BdN%_1;ae%D366 zEI!Z2JUwh=bz8Z#)r=o$#{YtO^c$j2O-z%;g*+FNS!3xV=I$+o*vmvD@QZJqR!cTu zVrLInYAg1#{08+yY0l-kG~BYN&RDP%QC8a0j#BT3!mRLmEa{_B%F>b>S}Rr)P0KXoYuh!K#@6y!Bb(l{f9+CW=f9Uzlil1H01BrIt%R zi`m^9^Txq-)XE#=Nq=0WB;S7_{GG-RFeZLb9fVvd#TU?KZBsdVD%VN<2c{?AOxO6O zujh?dawpZ5Iwo~laruVSxT4g!QGnWYF~^%Pbzaw4^>Te6_3c3FABt1+ic;5?q`qC! zb@2$vF>h4YSb2kBEroI@4R@+c zI@Z4u9(j7o;kQ^trwDHME&e6vwO|FVJp%?l!MAi3T}OSEe~sl|i(C#pO6lk=0d{%^ z=nbCAT-;{=B+zzVTlDsT6<=qiM~whmEURgoJAd!kwz)qJMA}9z!bZH-ZD=2QvOLmu z#*(mOw7!in)`uf)Bj$&@jn;kyqv3v|br{-^+&3>`#y3RTCd}V!rOy~4WdLIzGR4w! zkM@O~OWH2zFj_x_%m9Ldgk!WGjWD5c(f#6-{ZlMfuMCV$&Al5hYvBZ`<1#}H`q?8( zrT@|w`~|ROBWyJdNPdr~2wTk1o$T^Ruu1iEdd$FlR4>Ta=gF779x5ps1`&+pr%;u2 zed?I(7x5mRz*uf~`4a@553wT|i#0=a+4-M_s@M+zB8)LAO;y zDDsW@$kp&2z{-RB;3{|xssJIG@f`Z`X%X5GxxohuzO$% zIjAq%7R(vt@yBBES0DyraR>$z8#rwG_0$c|*vi7KAIwW*%}p;Ovz3Z8yF-~+r^_Y) zn@ex%#__fzu|EdI#MOwn42pA?3~GQO3XjHK8QJ9ZZ$O-K9+pgUnF(V3Bzz*j`ujn)@%wdlvr%T0U=5j*&{4*rR@ zv&q(2<>tKHwQpz&9q8Za_xd;T^K*Q#k^YUoLGk4@tVX$B{?3c%f!~zdjwUAR3!>@z zUv4ctKOcbKg1sO<>_mP&4*2m-*9F)wbinhdIw*d0atYEz!I8`on9Q-W0Lf@=>tqPq;>p!vdS3Yw*A~ zu1-gz7*n~a1z#e>4qr+Yp+jV7@$g8hZWy-;VV`by(_6BW6Sfef!N3Cvde5rner{Sd zTAwWf9DZ*N`g#G%KpscnlVg`-9fliM6+dA;;D9d^YJQS^tI)?Q$?w`f{1pAe{roDS z$Dfx!kAUF>e#clvZ%qi5rSY9RD9l>7qp?*7gd9IE-Wnhw$Hf;0NJuJ45EC1hZT(~KDin88ZC5t5vCAcf+tXrH5ety$UY_9T2X-{mJRaBgGP(xdB0kg zfnKk*{41&7A$jXC8PL;8s`fcF2}&+ckGMUSzA9k)GnC;zw$|UNSCo{ygo3-CB}{yY zQtft7w|)?SX_j|d<()3?Q=oIAL1`^$WB#WmpxRYVUV(J&1Db?8mG(tKeKnk@qg@qk zf&Ty-wDiLbb0YK``=%R7R=oJP#)6A67hx+10EP<6VyUNxetN7%Jt*VP!rl>XDv7^l z`ajK9aKU+-qS$!@uXE z)R-0$HW*2raVNi&jE?+5;hKPG#bhEFfq`|YC@@kNmMm=iW%P`U4vaN zCOeRPRsIosh}?lLxD(CFEk>n_f5c3*4yUxkgg_-tCq&x<%c3cKUTqnoVORt@E>1;@ z_NGR3D1js2Tfhjw;>+_)xi2MUGaj;?RSK%wbLA)BRzec>7 zR-dw`6~H-bJyJMN#%m*97vUdXYw>7CZoJmvRc6F1r$PJc(#uam;)=&n&iT$%Mn536 z<&r@$pK~({KrJ~43{3f!2bA+|0suXjjV49gkr-7(nBRruQ+5Sy!PwxZ8;cP^v#De` zvLip<++ohSe-#uW2RU{PuQa)6xt@EHu$FuVK9He3nJFxqVBlfR$4=YB01YG3#wgO& zVgIP*#|MS6l4_@(;jo35Rr#eb1g-Rh7DzuxKqQn*X>XhRyFhYZ6J!mySb;Bs(Y8~^ zw%I=p5UmVwjJ8c!jv2#wm;j7pu9B7=@NuXWf1F}jzWr$WbyP3+C%!l~Y|jr|+*a2d zXxbhPPgoLYdLEFn0TBWYPPsoTyR*~b_4`Qv_L+{ffE1jB&Xs-yg zedUt2;(IP>8?_W49BDc-8h)WE0H@g!8}qD7+Ri|%(=;*~-qJLJP^!?wNNu;%#h+z4 zmO<(creZnR@;}D1(eT!$V_4=$azj}pRXiVC`{sQomcC}+<*59kwh=#LX{dEnxJ4>l zgdTv!F0)i;CvfrD`MMHhVxq|?&_??fg!LV(So>D=$%{XPrPxZ2YMrPO*SU$-9 z2u?bz(W6r+Aumdu@%!xAaGi-gNQ=FS{lgxC?Anrz&==RqA#9T<#HyS|?Y{6zS6DZd$Otz|OuF4KFWaSVW!TR1|>{m8rY z?lA73fk*3~&>##S_-#@jVrH^7mU>Qh6#|LADeYF;9K_p$?Z%kXtj97p<_O~K8GK*d zO!#{qp9OnvJwG72$tquO^Rx6kh&19M){D7Y0{e&TQy@qZX?}bb%HJp-4&Dfw;rF>8 z8Q*JMKlpwc9FT=w4P6(5&F6m(A2c=6JncnQmccmhlXl47AviqPT()P9ZqMuZyf$_q z?0%u;b$sR?U!lBYk7~~zOd1R9^J9B<(>~q8_H4GwKjf)d`MR3g!^+hVv?p;g+H=A0 zDtEy?U5oZ;ch%!;Pj3v{7Y5xQG_XC-=Gucn@h-PH{eVIz9gq$K`@XAOKHL$ax01qr zLoVv!;KFd~rm$NEj_RX8#QvA9a`FHMV~6xPw*SQZyKMmIhr5VEFC+Ut#_y7y4q!iV z6>ZJ!;zh|Gjt?>yH=w)=1pq;flOtd@x}rKLgF4_)0LmHK)*0FUU=PKHXmLy7z(CCZ z4#p>!4XLt-S@*c~+;iYKN?Zsh>+`U}aT3!2{s6$VC6hOew6k)tT{g&G=R%iiQ*)$Om?+-%jm; zCS9YC9L<4U0k_R+vwS5x9}|Q=Op)L_UZ~c|P%EW1dv`!qv|pm3nx)ExG6%+K4i_Fz zFX>g^-94z`T?V^#8375wq&6Qg7XKA&FM;*)#@&|c%$C%#B%D9Zxexdf1XEUU`VPILmu0>#B7 z8nVv2*D@?weN2BPAEpiLJc%&BREhjSgpCDf!fC;MSj_~?{>2kYjQioLO;~bMt*+Hz zzG)m4kmuQh)N9>wbDW^JhK#r*h{=4*tDoo->eb zkTU-^%JuYn59J$9ICpx_kK76jY0hD4FZ;7FH}gMl`Tt^LcX+xEEn;g)PTex0tlKhm zx<(SBh|qBRT!$#MDw9q`n{WOK^ICu|5^_8#Zr*c5+}fE-8EuI zkC;I`lR(pNCjZpx&p=X1GZ%)x+8hE04EDI>&R)W~01)ok%MX;jDqc8jB*HAa_}mU+ zG~jX&?TgDFA1`(Nm4W;?mGfjmXE{IqoK^lLy`;F|VE^xJALKll++PCw?a>1AK+ccD zR)GbAnZgHlc(ai{pCbl}*c3%9^sA>5ZsZBfTJ!wrA=5~IhzSFgX%!mCFLy2`sYq$z zx%*YT+a**`jm*7-xrD$UOGud?jb`(F-|eit1I{f3KJ z(nbe1W+?h83TYCGJ(4G39Z5)5&XZ0^Nvm%^M*LqxzgA_W4Z?pE1+eyUcfj~n z{HN?scl>I;tKFY;;I2W%-%N_X+MM`v{Xp)05`XKs#d%r?ORl~Z7`KOTp7w8*f7qtZ zMne8k^g#&D)bejH$&|>!v<+UiGz@{|OVfOvv<~I+`{&=c3h-{Q03fbQAAY>l}K$Kv^S$oo#^s0j=C1 zt=zy3E!YV>(4MU2x$7IjC+}0fe`J+E9>Y9H^ki;El3@Gr*<*E2fmgm;;hq~|FFUdD zdxhI4pGo$ked}$u?@9=6UJ1Aiy;#GA(O<%~%m#N$*Xzy0>Ef>`zQH7CD#sJLmsShz z7Wgu>tzOsxXU&1v1u;Bsi!3rt!bC{HQqvs@uc;V-tdOTD#pN440?PqvDraGxo$nQ_ z26l=(kzcq;P9EvOi#VBJA+Myl(>UN>#InpwCxqIC#_8E2-6YP02wI9`PG0Y>LM@gI zt)XT&dxcq4=y-tkA69N z4PIoe6o8~EG-)dUM%$E2{!YAe21oL7G?Ma>^i-@d-mh`LS0>un3BZ9khsoj!7JoO` zaS;k)tK?5FO^+$ZIVArE@S7PwihGCQ-Dtzs(igeQkM56PppCD$!u8Hgn3<+6^ks6x zXNI-IwZ8<<2-bT3l|9?bNI30@E;skV-A`Y&gCKqRr6e3sQe-N|v0Qh8xlsJsMmc=5 zQC_kZbS~YlW;9-qgvAk;ho#C9j-vAzz*HUW93=phZ=PHBfQwNE_-#Jit%4eo*y{q@ zdcm(WQY1C~OlSuCrg9x`B1{KvrTKwH;?3kxGiMB*-m2_zIJ2fQMYKmXLG7Oc{h<~3 z{wa-2sU57$oOX^Y2xV{~=H_}Ja@jyQOV3)#1RbP!y10@qxLP_SQQchHP#APE|GVhH zDjWf{{X21x7k6aIc{0<#!Sp|_?0ER6f|FX=(^y}@CRS|o#{2Xn&yEG_TIlKGm1gom za_l_Cl{wd92WnjgaoDWSo7&`pnn&s9I^$Sl{PV{6r@%2IKEv{MQYPfc!>_cD95s@E zi-QQ`HlUFsjt9rVWCHhNNq(qwy%UsaF7Mf5`I9w zMo*zbxa3yew8<6mZd6S9C8G=Fg1TIpE~nz}P?wLQJez|vTwDNnP^5$%{{qpj((sl@ zIkjMF#vMS;um(V`W)6D2Fe73nl&|2^ur|XnZN@)KnA#7x8%qRg9Se0~kj8Wf&>wy0DzJgWD_es=!{op zxfxkdJFDzP4bFEZEH1f5WN{|u!$k~VL6ErB4f)9t#)3H3cGz)SiJz9xRu``LCEN!^ zO4_OWL^42zi@mBl2!cQg;g$OOFJ*^Vsi)XEJ*-ZzThg)5`FbHs zOMAAs8H<;Be*ArJn`y7AYeomX?=9!o8tvL^+!Anr|E|z|Z{ok>{UEz=tCt;q*1i8l zjR7IZ18OxCJ=F{vTl#h!O&>)BQYdFHf?Yqs_TxH&|NQpTsaHxKL@TK^bSuF?^CedI zH@$!N3!6~;FGH7;smKZhz^o(;8H~Vk{qxDw2j_2}q-$E1;flro9e)^DYhtNKloiEw3{qA5!3K2{+hJUzcFB`9y!ol-Ne-ui{x7*t=}Defe$gt& z_J<9MUI9(7Ozu7<@q>T-cpBUvZiQF7_mAze{M$iRBVEJxYkE-xQE=(zL2$Bp7-c!C zfPWap0?eqUGJ`qD53k+G&C{SS;~IoP?os4?=BSVN{||nD;=uN1rM=L((q1-0TyEVW zgO!}8#oCJmj|}cXP|F!=0=AqvU=eb^FAp%lwTg)=RkWJ4|t@IUxT=U|( zMY)_Uw`!j1jz8~t0(XH#_9tS(_pIC#fkjn1MvJu)+{E=?xrm93;hpd4`7FlI!22O? zMgzniZ&&C}BG5>pJc)Povt+r6@*Ki*QparZnmCeWc_9xOlaajImwof!LMWnI` z*s+ZIf24h`094x}^%99@Oi!adn5(}ZHRSar!_RNkL%A{bI^5!fkv^Np z7F(#KF!WFdqmf-{H$55N{fy)!6iHr5RN#z2L^mR|-hvq~gaAT)f*4Zq#O%1uuT4jq z`fL$40-EvZw*-s@KViVRM9=C{Q`YEIkhyJZr))H{HkDw<5qg z`kAQPM0Wv+O`=!GQ%&v31ij3kQtko@IOo}Cl;OWWfW|=qxXX>W!C6ZX=O+R{Ez7Ij zrN)A_Lsac00D{Aq%*z2sM1y$JlAD7kX51%jpTc-JnI!&RkxIPN(O8PTO!w;@6bWpp zM9ZfV26(I{^_1KS)5oJV%IGmMf1ZzITt$lu$+7+EDDn&I=Z)cfXJJ}|Gx6m98}`{U z#&HzkiLG>`t`JB#ZVOA1Mx5u7&b(-Pp`OGoR1M;OPo+GVR84THm*_yGBIV3SJ>Bgl zon~q)du<6aZa<5DsYuISztQ>{b7P|w*%kZ4{cCOi7OQ1DjIeAV=653Dw~W?zkO#SO zUubZ|uhlNSd7X$Dn)rMC4KzJWnNAG2lc(+oYdm0R)2n?GNS8JFynh zh6%-AWvGfM7$c8mOyH|VS)VcfoZ*utu5$Zr%t6H?`Ij=(%aEIYspN3EfNTb9V?<=&Ml+%Q6bPOXp~_q5!xho?BU7pQ(Vz) zj9zElzDKIQ_GwTa{+*4$Y2>I7mZnjC5WBg;o13s3TFDC-K1%9g%wzuyaml$}mSATi zR}0uK?Hb{$(|&=9FG6k|&sS%Tj&IDz^VMn8@x92d^Yhg)b-XX1pRZ1Zj<1(^alP;r zPL{2Bcw8?>FysEQtLVNfqoV=MRZEn3nGmXQ9ApRO@QS+;@S+lxcL7DR>*$YSVR2KS zkUE9VZIxJHPrwU@O*OSfIRJ2HZ2biXbA0-*E1XDnHL+ z#zCM2k*)1msvsRvfrx0^c}y!Ii1@Lxw^ZFRjn!W95;aR^LNk`qJO*^^Fj6;hcAyk$ z#(nhBe&?ewC{LFPcIF!kK9pk9Lyb|;yVOjr<-J}(9;cFv3AstIQ7=-N@b+-}qHds~ z{uy-srtNgkQKEwWko^7ynz4h-W4L}N_rV`9-_4&V+B0j1Y-n0F*nLwp*Y$l2tSB0Byyc$k1Ss%6OV&sAMyg!=4cR#Mh?c}Ga z8&h~MdwR@F=!f~pfd%R5>@C_yWkjBm0w3f85=QTjg7SUnD2QT7QjaBz^IaNcw9-Bk6hA!k!3?KpU%3 zF}`P!nPOZI`Ghftv35P8*^V;2)-mFnosvAlOwYY#FH=ga_-aJMtVS~!Q8C4s&TH^V zOk5*!3qL6W>Sd^?5CtKKe!-fmv0TnL%4z#h=g}FJ;rO^f;H($XZR!jw<2KGC=NK`M!lyY ztGTLCi&YKEK#}KIkuNcE^12rHob$-h12P_9I+4%QiO=~VvpL96Zn_`JpW_dQqoy0l zU+o3BxRSuU=PbV4j^8iBfvTag1{XaEThoTXDxK7*78N0RF@!0?Ig~?qz$0%yfzk91 zHvVJZ>hlqd9*fQ~?$k|>ktv?Km+qw@90gT2SOsUJvDS(Z^ka`)uM7u2*aXMnWSb`5 ziOE!6K72Fku;Yo4#8h)fH#%hcaI7>YpfX_}*5DXVwGAIx%8}1j@Fsa4QNXSs(Qj=N z)N#|)337*gSJPW^h5T0AzcCu#!0xPVdPMEv$9>-KH7#Y^Ly(|VHR)&-SQ{bGFbK{B z!k`8WA8-ZEaFL=$@+z2th@Fb9vb>-|bRTlgp;T-}4t3TB6AcE~b#p*)cPAnbf#7T( zmfteuwK>pfZOl_4m$bmp#e9`g8>1kJ2g<9(?6wKN-P=}f`=3sgjZIBx#_a=j^8*+S zUD1xfaOML`MZ=q#wnzoE#}G+!ZV@1?2}XT00xp^~$(c!NmqkJz$#l{7uMzYjxWWQ+ zwBK6scbK#U`H8tk(k$NrP5=j4xC`4w7-`fWX1ea?z0uT*7~|)_h|l8FLm&f0@y#fb z2sOsi6H-NJ_*4tW8FoL-E3X@I_f^{&zl*k=X8NDGI5mC#UN&|_b621v*F^u@Hck>f zeM!3RRKLqHZ%8XVO^n6S%@2l}C7wITbgQCZCK8`KP=Uizx(CMmh%{(~JwnP$ld)$6 zM^DS&GX41MP^;+0h_UoZ)!=TBmpO5q61)6?Jm*fxzR25Yiw{R+k5-?sH-#YeYBg=W^{rGX$61OF(v`x9g;Me6l}4?Cb6b*D6Xz`yWX#x-`3rA zd)vFzx~_L?{Utp6jL)<_5bLTFL5&YZYN>O-Kj(YCnMs(?-TU|2{~RUf;d>shbDis4 z=Q`K9&JV+VXK_HisJH(1qrG#(y@4CKQtv`zcL3TQD}SY9KNV{R3e``AOlhdM5$?*i z6jSbUsLle{WnZE!6P!hESf;@GOc zHGA=q0X^bdWJBPZ?W~F{oh?&-a0MN5FP*PuA%mw0mcQTjw*N>NE0kw>UFNnHlkS$I z6sxwen)K$;cwJNER^10-W`y`VkQ0}ur8MTpHHaukia?2hb|vI`Az=JKi!k6SY4hUDmGc^0`_UH*yQTBx}dKe?7f5wREqnxvg!1FYYW+{R$Cps*L?Y8zL zumjHKGRzo8djB{IyWpl!+S-dKr0dAUh_h&~QLPqIokD8eEpqvT1$)AVq$BuR%E5u& zQXE+S7KRKGS4DgGWDFB=dHTQ84xMf%&XE8t7=}t`HnQ4&*B(*K|9*A( zM>-(uR6b1W9hEEOjg;@h9#QYBxf*OVviBn4bh^LCtfHZ=JaD0XPcplioIec6KEgQ7 zoI^9KlYyH@>>1^1UzUA@+eS4Idmux9p@A<<(0y~aeU-8uXgZkuHaj?V@+rY_KpJs7 zqO%{>onhGE81FDBnfIV%wR_zUd9#7v#Aw%UF8n>Rp8XaNO72K}mghH({Anxz-+{B6 z@YZtHewo<3h|P#Gvm+Cm)6Kr|64}@6lYLEXRLZ%+_Q>_JUpeI$WCDep&D#QoKa{c0 zzf(b7HxOY;Id}?oTj2}QwvwIC@0hQRe6{7gO2_zm(rP|zI@QQ=NPDXut0Pw&UVb0! z@%79Y5L6AVBq37lwkT{mwl{z6bk$HAYX85qt=@ldhT)S`FYVrp-OWjv`i9*4CLCTy zeV@f-=PxOA)-7&z-`H=e+i^XQn+9k7T7_c!a}P;P*==7KY6AkEo`Ku1Ra7wdw4ad1 z>8j1Z=T%BL?H0ur|DIbkaL7FNu!`7eFi8BZFKRW8 zZ()*frQU<6Q5Zo@lL~#7LWx@?Q){FojRbnZlQvL>c{gjtqs|pPk;b+z!k!=q81?WZ z8fsyZP5C=FJGxIRO#5&pA|u09Y8@BcO|da7#F#b1>cS+^vSNc%w2>Ftv>huU)r_Ll zfnp;#V%?=|UYrzr#wWIRs;r(Cn95haTFGLMpGGW6*4(O{gPN!MR!A2NHwd?R#wu>8AYa9lt`8pM6ewq z>0|indAZbCg&}M6jha`Y@nyuL1Aqd?-nn8z(_-+{w6N7P5W|Q~!{~Ui zo(K7N{^4qGWgE2ZG8*(919LUV)Q!F_JFZO@Gh=X{5yzC>(C-E`_OO61ysy<3u#ghVv>yHc=ZgdP(5T!6ylw74@C^U#Z`7w z^-usKi0(s~3r6C8MBJG5T11^H=6MlpWlPf8F_KV!j@7+VMW{MMX?Z-#e6dWS z_=i!9;d>@e_v0iHpF1-rOL(oz3l_kT6ec5ihT|(FFbhu6FoUzi7TODBVgfVlAkvHs zBX~nKrb2RZf|W9IcOi{nA3GCH5D~@gN7?t%-V7IbqyvpOcp}Nrn3utmZiz9XSGx6;rwVfi#{g@KOUI$L`OcR zD#h6Nx6_>By&cbTzA`QFLrk|I=6bZ_X9BC=*?cXaN7OLWw4*xtb?pOZ@;Eb1I_;B7 zVLM%|9QR^h*E`s_=hX-VP836b!d%G0;Am>oMVBUP{?}1npLFL>q;4KEnW+@QqF1wn znQ~S5-&lJfSMhDdP}j3xZU z$&s#|NMf=slMhY!f_UkjWm&|JvS%j0?H_tU;w6S8d7XdgQT#4reHuQ)9~u7Hef-?o z|L1hQ>(GjDA9ufFdK3fT?s7S3cp)%bx|R|YgcJ&zFBpf(CC>lKU-5h`x>CKnzq?+W z6ohkCu=Caa(lnx7e_jls2}WU_(3z*@GdM?n$RcPb10nm&dbAV{%KttapJAE*jg{}| zctPmKpuKMr!@HPqTfN1m@4VPc4;bhJT$XdCusxbovqC%9v%U_<_*w(wfx~}mET=wj@FoE&6!WkEDFRiQbA2o`r z$V?>%DUSJnr+Hco3I}K0CTV#w^Lg`8vr!YhF1(kiz~TpEqBmCc`m|n5{5V|xRL5Wx zQBXDJzy8_07?Nj)ELmFr7Wk!c^V!bJABIby*`FT>q^|0annKH(bi;mjo)!ZvoZh5)w!Q^<> zl||E>Z|(_nsyHKAPP!MSu>bq86bv7W&VH=A_xH!>BZ1{u6fhUKbL%6$G)f<7qH?P&RLUHv%;SuPg?bp(9+ndW*WRLC$rpw*PO!}QkA`+>TFvap9g{Tn8*P1{=X71dZL(D?R9O6Lp!PPV}*e0T0q5LjO#q47b zq?A{aQj-|6?@HFBa4&{C+&9w@;aKmCqG;WVI=?li$K0%8!3{%Dp`9Gh>LOeiG=DdB z$}G^+wSijbd3UA!o_A;J@aGKv!u!SgW)z8|DR4G_PsDm0W1~5Yqdtnf#=Q4(J8eWh zmG$0>%%ZyXA^UULhXjjRhA=k@ejJ>>()O*v3`Dk^*G?v)%?u$d9Ru_&{}0DgWC4*btP z$@)cPl7&eTAy*Z}P1&z73lavv>{};V2*vSML|E_3_zZ@=6EMnR%wLIOmA+Whz16Jp zfP}+LNX($VC1ZY0okgxFFr}DY_Xd(D4}v0?=21mSU-)~0t|K2qQ;MDLM#MrqG!gF0 zb2f`#5G@gzKh)AtIHl{Wxh2r^G(;+t zXz_FmX?b7b7Xd(Y7tQ+a6}(*`{)wktEzz?L2U<7MP?!wI;Ndv!%x+{>W1dnwv=Qz@ zC&Af#9+Zt{D$+m-pbYX6ctVL?u*3c6fU8VA8szvscT6}%CE#kPd2F3_%`1SjMmE9_ z>XhC+n-Ib*)#nX5A5j1Je8W+VY$pNJ4HA^ceZBNUuJ^@BwS9t!k#)?{b1 zmJ_G=NoIKQ3r_I?y=iahAw!}Iri{Rx8)3Fm2y|XosHunXFu|A2$8I;YG6NJ+KF7bk z9S@=`a8TSm&)kT!EoMTk;=!Ahuk!CL6bE+n^jA7pNP4pJHT~h+ibF=eXznfkUU84= z55hAMK?Xx1m}tS{Q1d;Qk#2$v z*f+J;c1?)w7H)ODBkcSjzmx?aM18w3+Y5$%z(XxD|89zhzc;KVXj*u0Ut7@2FwaHj zOyVa%Q_`;p7&Bvmdpm9wE5Ft@XwKZhSZO~pC?*8`gamIuFqYf2{Yc4p z8Rq@iJ`uV}$+s#ww%(m6XYM62SFeGOTdM2c0B$^#aTG^3%-qwXJDkl@uGDN9l={rp zCCNGFe1f$_`}o0#eZt)EzM{My-P9YM|E#RzD5x3B2oYqj96p%`+duEFb2;tnypV`< z-W8es9O@He?IN5DiIROCV;5nBAC>6e-q-W&BgD#&V;{k_AMj4h|4g;_p11T|w({>L zJtuJ=q)Px4F^4<;4G!TNy%~L9Q6~YEz&LM<>4{Ciza&{?z{CjnKi0}5;%_wIOW$VE z1L7|+^Ka64^$WOnv;FdfGNN2ZM4Gnn$h^yKOO-g%TZndvU1P{sK(?2xKD}9?H0-K; z7-b&5-yw{6Ky5_CuT{>Vg1h)hJ%kVFKbW;pK7wmUf-@Ct9SMdMm#^-xoolCFGb!kEkXvoQRSI4ZBHO z-(ub;!OgR^*1)f=WO01;8Pv+wD@=Ydm3#vl!>{T?tf)u-;70BVU3C@LO!u_5MuhQ z2}x=(h;Lq>b`91!wa@Z8SbNJkWH~Sz0=wpagE_d_F?^zf?n2Wm8sRT zRqbBaFf&$}D1po-Nb$<|qi9)C(5;5*G^qBEK43x9R`wfRg~ZYgoQ0JzWO%UJKxeI_l3 zEV^lXjojMr0d~%2vr^_hH)J7WKM}eXAxOfA@7b>Y zsjpj!U;t66&;^rs`S(e{xxGKu`ywFS#s|~wNAN_F%LA-kiW_C+&UU<~{{d5uQ)T*o z8sH9eL42kKjDm}Ios*l&?j%3Fl4vD00CJxqPy1sfya?V{4VNOm;c#G zIu-fKQ4TdL+aW#(^hsX zjy?q%;;&-8vN(-meRki^Qz_~m!g&^?p@4J6+Y!B$4aKQY77!<5+&j(fF%{8olDkEO z_=fyT<0i_~N>I-n3%Pex_kL2{`!2CUk!Y_-gxb24H0~C&Zxa1Fl1EF{5=VoBIQ)-2 z^uq=El<#z(K^nEOlXz9G4uNcU6@5rL?bPq!gT`))N#L$}EXpBcp3{<7^0`F+HnTWd zk_j_tdJb6oDNGs4RE6lcGrP%XcidIfX+M?d(_U3Vx?%$d8we?&S!vH{?lf226bS0e z-S$g?0$WXBR#9sI2QKcbUFycLak*)w1ca|PuG?O*t)`_E8Br-}p=ffyZCiyRQ;PKZ zFlYm&)vj?C+iJ7Ak<7H(EL{ZP2rYj=z_fb5Ex@$;hzZc@_lMk8s{xF{2)Q>SB1rXd zKnvuAy?3zysMr$HY9Wv-&bjv^ zDXK%Gtq-w$m&G@hsRgqF)zrIZ3%u}ed=C(?=U=oAmh63zh8sj0cv}1VgEC?c>lg9! z;Os}&KPDXOCkMLTohB=MvzHI*8tyUg(<Xxol?o}@0 z9i%kqu~tN(g*kYhSp>h}0wu=N#t4M0<4hhuF9f25qI(vBvB#)*>eo zk1?zBW88~?v`;@~k>mWQjYV!SR_>OPIL(u{jIB3X-VZ);2p%0h;dsBGl$>#c#nDcJ=3b{oP%0CuyLig_tL(_|+!++% z_M{@5o^6uDF!2TnQI)g#3MMNzPDP1AAc@AYRmj5lMR0XTMpURm-@GjagTIyh9)g&v zR9^2#x8fS0Vps!)B-{`ifd+uiJqDrfy9Y_;2y4`HXp!<1>aV+8 zOV~Q~Pb?bQ2L^5jUWGpzvl9yza)!N<{h8U9O$!+1Bu-l}Z8*^OFgCk+I*67=cxIHh zo_ZhCDz}Y@=)fVY2+aRSyoMp`{;^VluUm0FOEnC*yFjz>1fMg#FrU3*zr*Xn1sE%%jeM_kqusbCh^vI*%Byrv3FpCa zNe%bawg(v4Sh*|mO`PY_r~gd<%b&ZRc37M+qQEqdOjIMarVC;$M)1w=_E`8bs6Hf@ zu{%}&HK;7>9$I);iG>6W*Kb&S%6p;ONKuK&*Tn+X1kPps(C_`*|d)a zWDYn-=D&a*6wrf)%@qZJb2G*wPPXvQY6qs@-!==C|8UIzQ3N}M0%CRfn@(~vlS1&) zoB+6Mfk!zTb{wZ$vL=rYVl))#@5_S^a(X87xC(a%OjM%A`iTD%Z|#ET6FzsAfHAGx zFh5T5e$?NF4)ARxWp_IRWJ54d;Y^5kV*k~a^3<8xrCjl%O8Gm!8A9-Rnxo;>Oh4|h znn8VBeKCXD{g(z6ZFkq$*J;en;4$O*$ol`lc%C}mcaj!^y?Zh0xn=`pkg)?t{; z3@Sul%jpq9x1qyI42uB~J}9%{Ht_3Npu72Km;=HUghNz>SEcxZ#`_^&U1CXnh)HoV zC2F3S6s2)rpnk-Babr@CV|mxC!Hmnu0t9BNnW(IM=;Ge#qp+rUKZ#Frh#R`fCV1$J z*=Xe^c%QhTi%1ZI{pT5EBlmxs`gC=_pElYn_wuHmv+*LEHeTM;Z*06Kt&NvA^^Z3G zh>8CzVXeFoE_EOw%LvkjHj+?u_?Ry{T>64y3vBKQk^Bv*zV)H5C;X{X1)r{+1!$Z= ze_VYNN#DZA?ph_<+%Gx)Ui8Q$xE+*KG$koMiZU_e!5rGGbv8e&#?-wd>3g!s=%4Uf z5n@-a^+kLCC^3glv17dhqcudIa27m)bi-Nngfv4t7_(UKc5`kG>uUARbmL}V{k!E) zBx)*q4@)ETiE!^nZ0?h@jIWVDfb7i;-ji_P+)JExGPiqB6ADrb)(JfkN=f zlG_dfNOmWPO(jbQ)LSt+!WW>g&--;3j@-aIFb$O?Xm+P*n*_Z|wUq&rL5}&S=-5k_pEo_C&aF`6`;GJ&^}tgVb+uZ2M`=IBoA4d0MIQ@UXM#mH{Qr_d zNBnP9mham13*=Rj$|Hq-04ekar{|D8q}zsoKZY!tqq}X~IwdwUX?B$M5BY@->GlLs5Y<(R0*;)W&3cbEnp_`@Ca(|1k`Xj^ zQ5X~*K%@5+b}G^Gv~%k~OnSq~fSDLHD#R7|NdwXHqqx=ZoVRYDwC@b8$NcF;=STf# z3{;mt&WC9EC-y_^0yN=F(mghYM1U52Y7CI0&t=(rdang_;s(#bu68&tra!1rC}M%OWs+IUwcXJEhinN82r;$jhnt7ww^Mh`O*a z6}sJ|3{oGFqGZUEO;jeoX*O7zA`*V|E7cv5WFfctxAa{eL={%dK>jqH(7S=Bc|NpG z))s<+Ep}524L}o10kh8y25vL{8x1VWcWk;}XGv&JQn9)|s7e6mhQ9gRj5hI{L_+Z? z;??l^+n5fb*5CaGB0Ji?6+fhI4q)&a?luyA0CO46g-2gykQB6cN50>**95jBCBDNn zx6`!hHVM!*y8A59%>c7}sR@FT1(+&TZJoFp)tR^kHkiMrF|C@)lpGog7C9Nu@jNS6dTBR$p#@C#Ia{bBbljsu4g(WJq zLO3^LUX`wK;VUe`_gfacLH_M`HZP<*?AD}X@|;urx>NkJkwa+A1cYG1Xp{^?^j}eM zzBEb{dc%aiZ$fPT=&5D{suj)(oloxF^q^DxA^T%C%v4*<<+qjZNR<|?8n=0 zWPMsok~|O7cDb)`pvq^`yrDnvWoVV4vGFqTD_%pltvv17Q|ZO# ztw4$%&q3zTxhNpjTEqZ|lx&^`!IM%Hik3}9bL=|Ww;%|hCQ`&6DHi(Q@vuUwh#{hP zxDngK{=LGZaMw}N7Z^uoE zIhL-Gq%y%WiaJTi={d(PWUb~bbpr|ASuY#jtxnH;q#PVgxjH51fAC30+vH0okk*%U z>NFEL>~xDx!6P2Q8Le^hC(dnkj&AzT{g0JT=I;NaBV9tBHq#|@*J8kQ%hfY6sI$u{ z{sUwEHvj$rPJvNFpHqL_24tO5T>v_F1ic-#dcA49q3dmcBLTvALX$!taKP!7*PA*8pMg1)SoVD&x*e}}C9LS>kCw}p#7BH>Qdhm=p$Q40u#{)t(yDCom#v?x&K_=;v| zaslkf75uvj2@z~Vf2Iu@wJJP+Pj@^Z2u zW*W&;HECbC8!loZ1$S;N`iL#+pQu~H?0#=tUao|0345o{I>DJa8$@LQp z7@Pgh$~WP5&F3ONv3}%RZz$$7v+|r%R=JX_s`6dy@0S&to!Bt|1UswierI@7Eo1w& zpbbj@KP9Al4SCH=F`m`)zo;bG=;;(VK!o5dNNcd(kXuNe)bPY-2a7lW)P-FFVAm{j z5?IB~qN;pg(_U%Fi^%VlxW*4`G5QVff>8Q7c1wILtuX5grUDFEk_UT!FL+n~8{4HKv$&#*ov<{MGQ?K>?BEVmQf25c#ng zh9?<~f$*fgU4J~0Q5pz$?{&IwC7m4&t)b>ICPM@Cim% z7%qQ&{YzoqZ#loJ_hX=2+F$KIC}3Waob`FZG3CLLWy?J~ws0!U3zt8F!1nMql*1JE42%bTB>{jsiy812`%_nvzOww?>tFsA3%Kv0XOW?p zZ#GBj$83qn*70y9cd2y!Z&ux`=J;=buje%>eGJj$k9N<1r3H@6m52clM ztG?}4+P(<~z8uoROws(M%+UzwbblEo10kW){mzM+a0^E!oZBGVlCrndl0B(yl!~BD zAjmFNXGr;Uql$ZBb$i5!tG$(+r?Ax)lip8okv&@RPE>RPrEVB%`8w1RyM1=&g+x;n zQa|Um?ji|gh}vQzw3(s0Re;G|0G9+K$xHdi-n^W>GK(OB8?y-VQn!4VF~k{J1kG7t_A;LDO1R&)W&H9e)^d9ob`SXlX@vx7LtvsCOR3J5R(rAL7kv5ITpM zNSzg(!sZ=z277Xkt$&~j6D-$O^$u6p?Um77vx3+YL7hkm;5zQuimSyFYrFo3F|ey= zwn#2OV-TPVMF$SGabANVsPyMYDf!6p|IJ?yK5 zy#f)RBs(w`6%h;l2J#CE-bYXcBh=QHJg56)Bq=t;Im~WkF@Q?*od(`%ND_Ps+APHA z%(d8j(frrWI7CQiy(!9aryopVg@v^y(;y!-$26uAdu!-rkF6;ls)SS`pK2*ykxR7`uP7p3)5@7PwQ1#A8!)ZB$pmO+>>#bQ z0zja&g?NmEcz{a%=KYaQn9(pWDh5U}#hD?L&gMf9G&bex(N^H-L0moWCk)G?b8!K{ zoqKn=dRQo#Wznh#I(V!*!zuvr%Mo0d{(h%;moj7DbmV}XV{v!CW60ZIt2Z#|(d0K> z9>4q%6b2*Z??H8TPa1|K=}CxK`BNfO;YrV*F^qAPfXpSFF&s^{LJGI*&1{?g{qB8E z_oJX;hCKa-kfepHw{Y;!Cd`Ptrw5=*g(rxsgnK_Zc7hoDe(#tX$~e7HEN2XsbES@{ zpYC7?s*-Ei%5v*_Z!zcR(DC4NY?KmSy|`jLT@8q zae`=gFg6%Dz#J$49*4u1IK++d6tJOZlB*{xp z<;W~UtRz|**mk6m_zsAhUjV(R>i&19ryP%28+5)T_?z{B{S@m1FlD=B5B|gD@3kGM z>VCUp0rT7iaz(!YhYilT&{@ zAkq(K0q{iK;HtCdE`u{+m@>KD=yXT(z12igExVc^`FCsZDigc4{8e{Bd6eoXjwrWK zQ^!z%mO|478Yui-I5w1@xX$`GWuHMRa!vSe2-m3r23eI}@iwMGC-Nu>PZ6jBcTi_O+ukX===Ao0ZnaYHxR1YgLJ;0h&J zLjsu7V<9p0iWchZT?U1Yfye#K8qvmWG)(iXAD0QD6Hf7G+P|TM%Y(&I_MJAy4z-CORAVA;tK^$-=2s?z}*w8}el41+J#GYzWE? zU*IzMP7g+hGz+rPzGo6AWID=zfRx!s){h8$#H%qjUf$HbHh!y#|Du*8d6D_uO~12t zd%YwqJv>ybjpNvlJ4(|-hnbRmzYAm=89zEEB%b^7@rG{bUNk)FO z4V~504@r|DHRjw$#~=H3{4p3wo_)%!5N!DueK^ua-6ttBXX6Q-jc+tYB+6p5U$5-N zQW6gw{W%{qiKQ>59)6LU5!D*|b6Ac*x)|_dF81DzwC>34C04|x5)JYRdFBfNO->hcmne#~z{!ZEA+%vlQ>eDNKD1tLTU+gQ?LObamUXQi z2^ASz91FF!(6V-uXiYVU=A6)a{b_A)oD-_&Z_ade zEiG;JbEs)u?VM}c8|z!bj=&)G+#I8+S<{hBh;X zV?Ow7yXzSKCQ6z^Ti1cgb!uP&Ts&Jir{>a_6v$@qkRFsXP zFoJRQ*EQF!W5_&g%MI782Is_=X26s%U^=_if@!2cSR`LslfnKVmF{vl;(nVXofi|l z_(dmLQSkbO92)y7_~-pPJ&#<7>Za^4{?6=w!pL}{rKPDoKCeDrzs4gQOIkwfY7=#9 znTnyMRm-B6g`*4O3l=Unn9 z-?@oqs41I%bN-FzQvb*=fG_)>+c)*k`{%=-`RBpIbNk<@Fl9}D({(#-c!dp@S$EjR zueD*y|7-p&;NOMyYvQ0HnAPs}-sdYmD|FVJvp~dbB!WOSw4OC56g_YD|NQBfF8-fC z{V8QsTs$Y_i#9VaYMUCrEqsHFL7l$UP?u%-|Xi|151{K>wW1~&sA$7+8x%>-4 zbnpxKAA9bTv~=%Ydf}7wQQX4)pQMWrUd_k7kno(}ev-Zo_t1Ii^b82{tc%j=dfa;| z(rMidzOyo&eoStd@8H!T46fJX-V-MtcUC>|48T3#Os6;EZtO^>cj9tHm7X(!_`gr5 zZI(-y()ixB|*zh0vCFSEb{aZRc3kI$50k##m8~=AYeKYPoTnGCQ_vmioVZT=W zE}gyycO&jjxi6>Fg$O(fUr(nKxEtR{r(crGuDb|f(w*GCXu_St{`NtI;odH!-*eys z>TzcXQ2hMgUgq;}2>L%)ctXJz{|Q0GOGr5U{3q!zDS`dvPVe&-@m=XxJo`yH2%LOV zgNvrlyu>+ieL<)1;!`j9(m6$cr6TJ27v|p<#^Hf);8`#=_`|#fQ)h0@UobV43@n&h zeDj3LsUyA7~Y-y+c!U*E&=XY^;eD5AI_@3-=F_CbP?&($6B@)LaPDQ_738#XcSAHi%R@q7O|o&H$xHE;`$hFhkd{8#gi zftM-B^q+D%!RIB8SI%+Zzt^8X&l$hH^%SA`(z!I9{*3@KW<0A##xtz(te(1cV*UdY z88zBEW;~ain13K|yzw;sr=I=)PdeQY@al=E|51&bsb}%XxaD7%KlzvnKnsI^I|aY1 z)9LOhqsob>oXALfE&dHp%>QL!&Q?WE%xlc~LG`G;bHIn@$e#p=lTD1oc1Qj zzpXW$-uD@=y=J^jd(C{w^n3Bt{`~yE`Ap7cRh_u__~ZBz^~~%>qi?!bPt^#1WyaCM zWlJD`;k5DU$r&%OxbvU!3*OH5wa}|4gU83#lb8SF-0-hD@v`yj5uV>eJ$L^@I{g#q zXHL9aBz!S=dC64f`3Je@dDv_BQOYa$Y4$u{Ecj&TMrJ%0Pi@UFVMueJdDV&G@tbYn zQ_uQP^0Rchb3D8;^W3yMI(3UbfBtwCRgYg$5g>VtdLpbte}~_9N!NrM1D8)Nxj8?X zmlvNrwP4YSE1;KQK7B~Ke$q8SH*&%;)9?Hx{$r+{N;mzw5PF;XMLPWj)+e`|41OI$ z_xtnm*A?Uh;3X&Ktr)js8uz=XXWOsR>BWqfTaTd=8G4*)Z+Pn3{Nlj))55g-LzVYX zI(_BH{LA3$G40+mIsdP6&%fyL=AY>PD)98~r_$-q;n#ZgbqoIl@!N=(#ow{;A5%Yc ze|)<4AmyxmHl6-{j`=XAoJ*!|%gaBIn@&Ubn#P5HJxwfvFWK^9HojakLPw&aubFi! zI(2FO>HeI<%se^HxCk$Ark*=qN~hVX|0R6Nf0atcE;t^YewK2UzLHMQ9Kr8#$2afK zxr;G;z%1~(_&_>+!$>*fuB+jxSLa8^>o>hIc-KUEYX{Tm+Noau#w~C0)WQ7xf1Z+a zSgKB3k~7}WW8vMi)HC?lC+Ra)F?jb4@U95&A>voOlTP<4{*oIf-aKLT)QYYXl7Y?n zKg`P$b#VEW8Ss_;q|-2b`uh|DpF-eM2z&~GPa*Ir1U`kprx5rQ0{haNWYqSmB#(7`HP&Xy4Lv2UJPH* z#&ZqE{8IM2WS;to7uN3&pQKyl-jP=!$H&&${K_U)7?-W2`Ss^j7#G&r{I1nrpPyt6 z`ptA5pR^pZe>NE#7PoyPJh%Pe$!E32)9)15@$H*jVSHZvRvY%(?VWML@|&MmKj-zP zKf{wMOgx7VCLC~+`G#%yAMAG<{a@}JvG5;T(Zs10ChBDCo@w3lth>m%U$-th8k7E8 z*1gfXw_Eq0tb4z8cU$*)>%L*#!`7YX*!rz|rghJ=?jq}c-MXu-`z`C6CGcPahBZH5}H5v{JH0xUD{zH zOJ}=bpKosab?Xwf*Wf1FTz9SKa}Tb4bDLWd4Rb3OM9=1s*@o6Mcg($}qp_*}?8bWE zT!q%wwy*Wg<(QB1TsP6?e$nv|XO*J@alW-RG}S7D^|dx7d~+K)4f4%RG;H86XW)%| zYca>Qa~syiSL>`ZzP6scUdYW6H-~_34$psW9m(veBUQIhAb;kXTHF$^X>a$h zDW2^7iM$H)DPi+#f1mvun?KA)v_4wA`aODfhm9z)`c#Qao7_U|eJqE%jgk3Ri_C2*J6#9sTo zbe{X)_%7P#g+2N=Z1Zok`PCf3&TGGo_I=yNd-*+jyVVkWkKTIqd-=WocaO>M(eHiU zf??sO{Du4*n}0KLf{C}@cw3=JcvO>@zX(M`{fv^c+|kW9IA9XkKlkHi=Z}=T`6K0S##P=p+VtM|d1dtDk=V}; z^9Vn(zuL|3!P!ecHveOMaI*4y{7TB^ADez`ct3enlUM(+U2li~gZW=3uW--H@2%I^ zjv|(y1+0gkUj1Ib-XKYK{YNgVHQA0_=H~a(n@Ik6Vb47{Ccnoo`Ifr*jjY$poFo5- zWAe}BHkp1iFL(3LFrmyJuU;?aV?2uU`mft)`pvQV$JQS*Uo(GD_%Qu1$wV2Sm;Xy+ z;@wucdGu)-|Gc`hlYZUSFHQOE-~idWXmtEc=4)+sl=0{LKI`-5WA?}VV^&3XJ_UU5X*!PjxxCF3-%`6j zHM%UM(9yY{-CMD(*ISpg-ui)a=sXvmA;S87nSZ&;zmWJm-xOcxbT_@X5B-jf_x8tU z64pKGq)C4){9hto^==u1Kf1iE`2I2Rx_qPfvHSE>iC24; zj%m+q6Mrl|h%BodZ~t$b?^|WkAB$JtAwE~XelSk_{}?C!-^PjGF;4tW8kVcxr%CT; zUTj@h={4QBgC_o1IO}7sat@jF$KrFISudv_Gmk!F<2`wWSNWO5`;i`uokw3GemZdO zOjVeqMrR4%0^&p9?b!KJZPI(CW?WrJRXO`MRG7$N`{{=^zT`9)4n}8yZ~r(>InR$1 ze{h`mJQfwzyVUY?Bs0J1Hh#qz{QRtmw~4ZxzaA%EI(;fX*ivCq8yzuxf>X%EjQ;cY z>elmppD;WnD4e4XQ@{~Mb=ly=h_T_Mu%9w+@Dyz+nKrZ+l` zq(2SD*|gKH`$p%7_=t(eGv;@#O~3RDZoiDK2x)E@C;i`!6aNbFg4+taUVGKNHct8z z*}tlt{TnMxQlm?XZ=Q+wa%bGdHht?rg^4sehWJ{E*L)m1UhBsxXNxUo>(AVBj7}uw z{KCfjEd4JswfKDZoA_h-XX!%c!rLcpIsMnU+R9zWe5|8`p6LU+;&OgW3 z(rRvU%q=~)bnd+Jb4q-u!E+PDD#F*bw))~%E>ASIN0%o`OY7Fww#5@|wT;*tudAwB zM%+2NJyRENLC)GxyDnN+6R%o+ZhYybb#uASK^?yM;@D*iDr50@d&f2LXgq!;#hBYU z4e_?M?MqiLkLUWYv9Zxq94m{*8@axNB(<(BK`$`+QRj=tFIl)O9*ZuoQBNA$xXpx4 zK22K}uWM>)ZfN)Ah>FMSTjFb)TCS;Wir1r)-yW~+*x*xJbXBB&?s?~zmQiw!L~%)P zU0H0Ldr@uI`CRF3eaE_W*O4VF6jxJ6C2d;MM5Jy_nQH|(Z(4Bw;~MT<#oHSa9j#e6 zFq#^h{nwrJ7vs+xFB<$_q1QjBU@{pbkp zDzqhh@$kwemDSOOe6uyjt1k=38=BU%`y^(L+hl9F&I1-S$Lm@<;<{E7pI=e2c+q9c zYL-hTQ2`^}Hf1z5aPz9Zts&0cF+c_q)V8)ZT^Cr1UOb*y+t#u^eobwCysoxA@vSj=M;g|ia2o=aG>>l0@%p@`wxhj0 zF3d|LT4>3dbuG&}jWGnc*MOo+yBfkgtvrOOlN)+MfMZD78Zo)fQ*ujh_ev%!U`(Y>5o%3L+a6oO!) ztg)pXc(gankAl$)OUI-hr}{dGK%&81BvS*{nJZIc#s}PDAPg$S*VQ&QAN$!Lz3q8Grk1I9{O+u2Ypi!?r7s@6 zjHcH!FFV>B>Y)P-LW6AdQzhe6k#i0)=i`|JeEk zy)0n`JC7A?`8n|>JGO48=a%~#Hnff2i#YUK`ovgCQkm<2w$$V)ntD0>k29>M`N?kr4a4kY)c$7>5en9>s%UW?y(AsTqN>v=>Gt>#I4}~ literal 0 HcmV?d00001 diff --git a/myf90/ESBGK/Stu_DG2D/stuDG_driver.m b/myf90/ESBGK/Stu_DG2D/stuDG_driver.m new file mode 100644 index 0000000..2f221ee --- /dev/null +++ b/myf90/ESBGK/Stu_DG2D/stuDG_driver.m @@ -0,0 +1,46 @@ +% stuDG_driver.m + +% script to prepare input file, invoke DG executable, read the output file +% and process results for DG code + +clear +clc +close('all') + +nel = 2; +nop = 2; +exact_integration = 'F'; +ntime = 1000; +time_final = 1.0; +icase = 2; +kstages = 3; + +% --- delete old input and output files ----- +s = system('rm *.txt'); + +% --- open input file --------------------------- +input_file_name = 'dg_input.txt'; +infile_fid=fopen(input_file_name,'w'); + +% ---- write input file data --------------------- +fprintf(infile_fid,'%d \n',nel); +fprintf(infile_fid,'%d \n',nop); +fprintf(infile_fid,'%s \n',exact_integration); +fprintf(infile_fid,'%d \n',ntime); +fprintf(infile_fid,'%g \n',time_final); +fprintf(infile_fid,'%d \n',icase); +fprintf(infile_fid,'%d \n',kstages); + +% ---- close the input file when done --------------- + +fclose(infile_fid); + +% ----- invoke the DG Program ---------------------- +tic; +s = system('./stuDG_F90'); + +program_run_time = toc; + +q0 = load('dg_output.txt'); + +% ---- visualization/analysis code --------------------------- \ No newline at end of file diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.depend b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.depend index 3ae4c43..7535aec 100644 --- a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.depend +++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.depend @@ -15,3 +15,5 @@ 1374818856 source:/home/manuel/git/aero-shock/myf90/ESBGK/TestingDGWENO/TestingDGWENO/quadraturemodule.f90 +1374992936 source:/home/manuel/git/aero-shock/myf90/ESBGK/TestingDGWENO/TestingDGWENO/main.f90 + diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.layout b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.layout index 5442396..4c18c75 100644 --- a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.layout +++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/TestingDGWENO.layout @@ -8,17 +8,12 @@ - + - - - - - - + - + @@ -31,4 +26,9 @@ + + + + + diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/configuration.in b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/configuration.in index 5d4fc25..72f76d6 100644 --- a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/configuration.in +++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/configuration.in @@ -2,7 +2,7 @@ $Parameters_list name="SBBGK1d", theta=0, nx=10, - quad=1, + quad=2, P_deg=3, RK_stages=3, IC_case=1, diff --git a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/main.f90 b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/main.f90 index 4b39eff..ae6e6d5 100644 --- a/myf90/ESBGK/TestingDGWENO/TestingDGWENO/main.f90 +++ b/myf90/ESBGK/TestingDGWENO/TestingDGWENO/main.f90 @@ -1,14 +1,15 @@ -! For Testing the DG&WENO Subroutines +! DG&WENO for Classical Botlzmann Equation ! By Manuel Diaz, NTU, 2013.06.16 - +! manuel.ade'at'gmail.com +! program main ! Load modules - use IDmodule ! create ID names for output files - use mathmodule ! linear Algebra functions - use dispmodule ! matlab display function for arrays - use tecplotmodule ! write to tecplot functions + use IDmodule ! Create ID names for output files + use mathmodule ! Linear Algebra functions + use dispmodule ! Matlab 'disp' function for displaying arrays + use tecplotmodule ! Write to tecplot functions use quadraturemodule ! Gauss Hermite and Newton cotes abscissas and weights - use initmodule ! load 1d initial condition to algorithm + use initmodule ! Load 1d initial condition to algorithm ! No implicit definitions allowed implicit none @@ -61,7 +62,7 @@ program main real :: rl,ul,pl,rr,ur,pr ! left and right IC values. real, allocatable :: r(:,:),u(:,:) ! classical macroscopic quantities real, allocatable :: p(:,:),t(:,:) ! r:density, u:x-velocity, p:pressure, t: temperature - real, allocatable :: rki(:,:),uki(:,:) ! Arrays for dealing with DOM + real, allocatable :: rki(:,:),uki(:,:) ! classical macro-properties as degrees of freedom: (P_deg,nx) real, allocatable :: pki(:,:),tki(:,:) ! r:density, u:x-velocity, p:pressure, t: temperature real, allocatable :: f(:,:,:),feq(:,:) ! PDF arrays @@ -69,6 +70,10 @@ program main namelist /parameters_list/ name,theta,nx,quad,P_deg,RK_stages,IC_case,fmodel,f_case,method,kflux,tau,tEnd,MM ! this list need not to be in order + !********************* + ! Main program start: + !********************* + ! Read file with parameters open(10, file="configuration.in", form="FORMATTED", action="READ") read(10, nml=parameters_list) @@ -79,11 +84,11 @@ program main ! Create IDs for simulation and result file call ID_name(name,theta,nx,P_deg,RK_stages,tau,IC_case,fmodel,f_case,method,IDn,IDf) - ! Checking IDs - print *, 'IDn: ',IDn - print *, 'IDf: ',IDf + ! Checking IDs (uncomment to diplay variables) + !print *, 'IDn: ',IDn + !print *, 'IDf: ',IDf - Print *, 'Build Nv points using Selected Aquadrature method'; print *, ' '; + Print *, 'Build Nv points using Selected a quadrature method'; print *, ' '; select case (quad) case (1) ! allocate w and nv diff --git a/myf90/ESBGK/TestingDGWENO/advection1d/advection1d.cbp b/myf90/ESBGK/TestingDGWENO/advection1d/advection1d.cbp new file mode 100644 index 0000000..a72d87f --- /dev/null +++ b/myf90/ESBGK/TestingDGWENO/advection1d/advection1d.cbp @@ -0,0 +1,40 @@ + + + + + + diff --git a/myf90/ESBGK/TestingDGWENO/advection1d/advection1d.layout b/myf90/ESBGK/TestingDGWENO/advection1d/advection1d.layout new file mode 100644 index 0000000..1053280 --- /dev/null +++ b/myf90/ESBGK/TestingDGWENO/advection1d/advection1d.layout @@ -0,0 +1,4 @@ + + + + diff --git a/myf90/ESBGK/TestingDGWENO/advection1d/main.f95 b/myf90/ESBGK/TestingDGWENO/advection1d/main.f95 new file mode 100644 index 0000000..58f583a --- /dev/null +++ b/myf90/ESBGK/TestingDGWENO/advection1d/main.f95 @@ -0,0 +1,8 @@ +! A fortran95 program for G95 +! By WQY +program main + implicit none + integer re_i + write(*,*) "Hello World!" + re_i = system("pause") +end \ No newline at end of file diff --git a/myf90/advection1d/extramodule.f90 b/myf90/advection1d/extramodule.f90 new file mode 100644 index 0000000..bee87bf --- /dev/null +++ b/myf90/advection1d/extramodule.f90 @@ -0,0 +1,20 @@ +module extramodule +implicit none + ! + ! here goes any global variables + ! +contains + + function linspace(a,b,n) + integer :: n,i + real :: a,b,dx + real :: linspace(n) + + dx = (b-a)/(n-1) + do i = 0,n-1 + linspace(i+1) = a + i*dx + end do + + end function linspace + +end module extramodule diff --git a/myf90/advection1d/main.f90 b/myf90/advection1d/main.f90 new file mode 100644 index 0000000..ea9a53c --- /dev/null +++ b/myf90/advection1d/main.f90 @@ -0,0 +1,51 @@ +program main +use extramodule + +implicit none + ! + ! we will try to solve u_t +a u_x = 0 + ! with IC u(0) = sin( x/(2 pi ) ) + ! + real :: u0(10), u(10), u_new(10), x(10) + integer :: i,j,k + real, parameter :: pi = 4*atan(1.0) + real :: dt,dx,tEnd,c + real, parameter :: cfl = 0.9 + integer :: time,steps + + ! create x domain + x = linspace(0.0,1.0,10) + dx = 0.1 + + ! load initial condition into domain + u0 = sin(2*pi*x) + + ! using upwind scheeme to evolve IC + u = u0 ! load IC into u + + c = 2.0 + tEnd = 1.2 + dt = real(cfl*dx/c) + steps = ceiling(tEnd/dt) + +do time = 0,steps + do i = 2,10 + u_new(i) = u(i) - c*dt/dx*(u(i)-u(i-1)) + end do + u_new(1) = u(10) ! periodic BC + u = u_new ! update info in domain + print *, 'step: ',time +end do + +! write answer +open(unit=12,file='output.plt') +write(12,*) 'variables =x, u' +do i = 1,10 +write(12,*) x(i), u(i) +end do + + print *, u + + + +end program main diff --git a/myf90/advection1d/output.plt b/myf90/advection1d/output.plt new file mode 100644 index 0000000..6bdc423 --- /dev/null +++ b/myf90/advection1d/output.plt @@ -0,0 +1,11 @@ + variables =x, u + 0.00000000 5.84559441E-02 + 0.111111112 -0.332437694 + 0.222222224 -0.577446103 + 0.333333343 -0.575789034 + 0.444444448 -0.347420245 + 0.555555582 -7.44608045E-03 + 0.666666687 0.313223124 + 0.777777791 0.525778532 + 0.888888896 0.566101491 + 1.00000000 0.382825762 diff --git a/myf90/advection1d/runme b/myf90/advection1d/runme new file mode 100755 index 0000000000000000000000000000000000000000..c7a519a4cf31946a9885eb510287a1b8c3ad9c58 GIT binary patch literal 13422 zcmeHOe{37qeSf4RTXt-b($ZA1Q=9AL!7*Y+N>oQq;P_;T_Q{}Z*O6*xlYOE^>c=9H z0(o@grk>SN%zU7gs9Rb%9R>&(1_W5qbwih_Y2Zk9R@s9$$gSIIb+!eD-AT3!cvF(9 zsT%I{efQq+cp`^jz<~XsC)|C1yzlq>zVG{f+&$h`X9B^!ZkLP6;%3h>qQ=V=&NAvsR7Di7B+3o5L#Y0`_?Q+!!d87BmqUf64MW+i12!`T42@2={MREsA^hLXg>@rDB zBU}V;Ak-IS`+_ZdnzSeUPx!yU_U%1*xb1M`0j`chlH8xp2I$&W;S1=A(W5IV!0cC&8?dyAVp zNsRBQt2t)gs(qatA)fS&7(N@H_6u3I+W2zsq`F!g-~OCvu<`Btq+;XS{ZVavxmS|? zLpDAFWtF2g{#uL6Sf`EufQ{d4fQ27 zHX)$dFFXN}?3Q}mJ)0F^ME8=={QFsyE!&8u(U_m%^yi7Dkk3zYdIQlVM1P0VRYX&W z=f^p{mS_s?{1~U*L{muThdKT6XF*da=Lb3cUqn*~=Q}xlhiD4j{2@-?B$`4tuWDj|1g7@_%pHN?2i;@T--%Zdz;<=} zSpYS)SfLuXUvwY4Gy!`&a336dhpMvgdo~XfZ-OW(Kebo_dDXb6PTkq9x@J`42kNy? zlF)yPWYzQ!b*X<-=nwB|gug6V*$&OQKZGO4eo*EI%e8Qc{83Z8&OZnNeWT?UrREs_|3xm;VmVX8( zJ@XQcRNx-w07f^{eoyWCp_;mbarOqkK2G-0((`2hcwjECy6x-Zpdp&m=pc=bA3@kb zf2OKLPKQVxq#01&n@s!l)J2!kemy1z5)55eg) z`AWOF&DR?r7cjQX(i`Sj0lfv%6MghR=t5fb@%M5S)RDgeiMc=%oe$2d-v@^Iu&e`* z@s?;E4RfqOn`dCloJXUFsBymbrf#~`(RE(wi!G$yeDT1zT3q(ludMEn*B^p(&Lz*t zKAIEeBG1i-#|!xFm=Durd^+t#-^8}!1IpL_bb2dLv8_}~?kIdU(hMZj$^odTqkqRa zX3)QfX4zY2Z>Q(}VUg;sO=br*_|k&UE0?aLFCuj&hdb`ocwiq9&Q@<%`F@9WkWA557;avI``UZOb z1dTgE8cxa2wuDBS7!2{_dT%4CjNvr+u-cKA6<~gZ=u+P`r=M z2?^h>p~rkly8(oRWf5jsXi({*LhzriUfp%Wb43BFeDn3!cXxMp?{;l0d4?GD60-DB zHhT{6DuU7hIDzng7jO)tLc3}lYn`5WZv#F8gF#G$WW zMq`EdWp~s1RqktO*|ZY|`@8V}W~r@y|GAP>X~~aDxs+*#`WwZheg@F8bJ7$bf4+t5 zJn9Xq-~aNeb7kEXN^13J=_@6#y4^=BXr5Aku2VhY@@YnxL|4j61Xd!j5`mQntVCcX z0xJ<%iNHz(RwD5KD*~Ns6z-5Bp!CjurXrrC>9>D7P2UQrtTO?7W3|Gu^@6U!L6^$I z{7&CnJ2|#d&~&z-vcUvwQPk7)4pi#{~XBg*v^K2F2?mE0aulqemx}M?+Um>z`X)~Nx*IaPYU?5 zfZq_%+1|dUrl);d+fO8;dh+pLUo;et`Re^o`Jdcgm*mX4M@^bBf8x}D9y$T2$4%HP zG5XFIA7uV$Opo{*_q1%+Lp_4(i6;Fgl70AK-4|wlqIyG#Ugi&Kb)o4P4B|$=kATfBo=7B^zbiH{5Q*x}0RAs&T9-8C<@|4OE1+_D z8FF26CT?PRNRs^yFj@;zzs;?H%H=)yIF-AV=1`8+;cnO0-3sUH!e9-jayO@bGpg{@ zN=&wY-mP$9na?5pFXPAzsJ zF9RW0rM}$H2E_&|_b+GsB!9%AFZai@Lf`w z7ZFxTGWre6u0OL{;dL`Ypa14&N$0*()Aw#Fui>)m%e<#?;pm94!-a8^dXoO8Ltj4s zYJ{MC4oG>Y|Kmbm#$OT7PetgHfA-?D6TsiVMgB{DnI|OkhngMo&i3DQ=vP-NT%)?u zRHWZT?M05C+btFPYaR=cm+{^Wf*uv9~Yt8x}W{wc`FuJ9i;t*Jr*?Qe(`= z?`+5f)?MmJ_yKC{`a?p0$Y-iZKcxYwC+UxY5KqRBK8I5|>*HWq{GIbB3o&ZH))cqy zQh%GpF%gobA|XvP**32bx4-q7NVyagdI8YunsY7D@(sUvpW&h>=L+i6p z`(vU*|K_N4wr{pl_+A3<8r;gVO3J}$E66WldO?1v^&V_DMp-=*1?5*U`5elZU(Mw6 zC7+L|tWxYOj+Y?IQ$7#ul_hxFmiv4@{{d_N&gZXVa$nBpueb6{^Z6A_?xXp950m?2 zKHqEQN#^q_ncUCv`Bm1rB%l9V)_EkK{~)vW67q9{xsUY4@|q~(}b=$Y%HuQ_sG6t#3dbF=kq#_>Mx7ksux z5t1fo6ny8o@d)@Xoc5%ja(#4j{_=Q^3Ok(+`>%8P<>N8|exZI{F2cVi?1$GXz!>A- z$^d6^x$vZ+53y8y{BG_;h3x+a>Rqgo)mAA&(gYuJ{_=QIZ6P}@9+0YK^QsDbSAqGn zQOL`4v)L*RJzReI{Go4^RTxK~=$G6J+ralpWg8si{N?lQCGfpCEz5C~Dt`oi6`t_U z^>~)^rBn`ns|f#Jitz6i;Xi_oROQTC*fZdhpU(B7g734Wa>x*s+enMeL7| zJm!fa;whhZf6C?8vvXqJ@^7?g=G#T={HzFn4W1`7UX;E@#o%}RUIM-s>uAj3=O!U9 z-%Eu!+alzh<5kDyrBn_c5q2nFjtc)a4A@ab{xtYL5_FWWflu*tw)^Kr>|7{fXO_zs zcwThH^@N`6?)G;fYj}UtVJ+D5TpMzRbD71OKA?3`ZgGOid|?|)L$a>c(;qv54B@aI zizl>Da){y3JlG%6BVqrJr|KGDxqy&H$9MBEuHk?kKgG=BdpJ2Ta0)6mN~2SMo+!@$ zIP1s4$^HmWaj~HYb?W}tW*7hE5Yz=6E1I?5V%{wabJlkQWH9vo# zv9+ZMQnnhNO3k!=!Gn7mgWAEpd!GxmX>E;rf&nte6NbAIN#3RZF*BQfSTIwInHJ%{ zG6@qMZ=Vhp6i>LGM?8C`qgp}Aqn?IwPYfK0nluNtF<0%MtGx%!BCgolGTM=YI>g4 zdRwcOY+hQrQ#?0_H0&6iC|nW4uHly;iQ7Aq%T~{`R?AawY#zv@bB9#HynGFr{BoWy TYs+bX{s%xV%KZNT;GOty2p8)m literal 0 HcmV?d00001