From 3541d363338f7297a010733bedc69638d5693893 Mon Sep 17 00:00:00 2001 From: tduretz Date: Wed, 4 Dec 2024 23:20:27 +0100 Subject: [PATCH] Stress BC back on track Recovered old branch and improved, periodic BC with normal load working --- .../Main_Visualisation_Makie_MD7.jl | 73 +- .../Main_Visualisation_Makie_MD7_Residuals.jl | 2 +- MDLIB/GridRoutines.c | 57 +- MDLIB/HDF5Output.c | 18 +- MDLIB/Main_DOODZ.c | 14 +- MDLIB/MemoryAllocFree.c | 16 + MDLIB/RheologyDensity.c | 4 +- MDLIB/Setup.c | 44 + MDLIB/StokesAssemblyDecoupled.c | 526 +++++----- MDLIB/include/mdoodz.h | 4 +- MDLIB/mdoodz-private.h | 1 + SETS/ShearBoxLaeti.c | 166 +++- SETS/ShearTemplate.c | 2 +- .../LinearPureshearAnisotropic.txt | 136 +++ .../LinearPureshearIsotropic.txt | 136 +++ .../LinearSimpleshearAnisotropic.txt | 136 +++ .../LinearSimpleshearIsotropic.txt | 136 +++ .../NonLinearPureshearAnisotropic.txt | 136 +++ .../NonLinearPureshearIsotropic.txt | 136 +++ .../NonLinearSimpleshearAnisotropi.txt | 136 +++ .../NonLinearSimpleshearIsotropic.txt | 136 +++ SETS/StressBC_FreeSurf.c | 138 +++ SETS/StressBC_FreeSurf.txt | 102 ++ SETS/StressBC_ShearTemplate.c | 216 +++++ SETS/StressBC_ShearTemplate.txt | 102 ++ SETS/ViscousInclusion.c | 50 +- ...eralStiffness_MDOODZ_7.0_v4_StressBC.ipynb | 7 +- ...lStiffness_MDOODZ_7.0_v4_StressBC_v2.ipynb | 897 ++++++++++++++++++ 28 files changed, 3216 insertions(+), 311 deletions(-) create mode 100755 SETS/ShearTemplateCI/LinearPureshearAnisotropic.txt create mode 100755 SETS/ShearTemplateCI/LinearPureshearIsotropic.txt create mode 100755 SETS/ShearTemplateCI/LinearSimpleshearAnisotropic.txt create mode 100755 SETS/ShearTemplateCI/LinearSimpleshearIsotropic.txt create mode 100755 SETS/ShearTemplateCI/NonLinearPureshearAnisotropic.txt create mode 100755 SETS/ShearTemplateCI/NonLinearPureshearIsotropic.txt create mode 100755 SETS/ShearTemplateCI/NonLinearSimpleshearAnisotropi.txt create mode 100755 SETS/ShearTemplateCI/NonLinearSimpleshearIsotropic.txt create mode 100755 SETS/StressBC_FreeSurf.c create mode 100755 SETS/StressBC_FreeSurf.txt create mode 100755 SETS/StressBC_ShearTemplate.c create mode 100755 SETS/StressBC_ShearTemplate.txt create mode 100644 misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC_v2.ipynb diff --git a/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl b/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl index ba3f2b80..14e3cf09 100644 --- a/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl +++ b/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl @@ -65,9 +65,9 @@ end # path ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiftingAnisotropy/ref_d4_HR/" # File numbers - file_start = 1500 - file_step = 10 - file_end = 1500 + file_start = 1 + file_step = 1 + file_end = 1 # Select field to visualise field = :Phases @@ -76,10 +76,12 @@ end # field = :Viscosity # field = :PlasticStrainrate # field = :Stress + # field = :σxx + field = :σzz # field = :StrainRate # field = :Pressure # field = :Divergence - field = :Temperature + # field = :Temperature # field = :Velocity_x # field = :Velocity_z # field = :Velocity @@ -93,32 +95,32 @@ end # field = :ChristmasTree # Define Tuple for enlargment window - zoom = ( - xmin = -500, - xmax = 500, - zmin = -300, - zmax = 5, - ) + # zoom = ( + # xmin = -500, + # xmax = 500, + # zmin = -300, + # zmax = 5, + # ) # Switches - printfig = true # print figures to disk + printfig = false # print figures to disk printvid = false framerate = 12 PlotOnTop = ( - ph_contours = true, # add phase contours + ph_contours = false, # add phase contours fabric = false, # add fabric quiver (normal to director) - T_contours = true, # add temperature contours + T_contours = false, # add temperature contours topo = false, - quiver_origin = false, - σ1_axis = false, - ε̇1_axis = false, + quiver_origin = true, + σ1_axis = true, + ε̇1_axis = true, vel_vec = false, ϕ_contours = false, PT_window = true, ) α_heatmap = 1.0 # transparency of heatmap vel_arrow = 5 - vel_scale = 50 + vel_scale = 0.00001 vel_step = 10 nap = 0.1 # pause for animation resol = 500 @@ -132,7 +134,7 @@ end # tc = My # Vc = 1e-9 - Lc = 1e3 + Lc = 1 tc = My Vc = 1.0 τc = 1 @@ -213,6 +215,8 @@ end τxx = Float64.(reshape(ExtractData( filename, "/Centers/sxxd"), ncx, ncz)) τzz = Float64.(reshape(ExtractData( filename, "/Centers/szzd"), ncx, ncz)) τyy = -(τzz .+ τxx) + σzz = -P + τzz + σxx = -P + τxx τxz = Float64.(reshape(ExtractData( filename, "/Vertices/sxz"), nvx, nvz)) ε̇xx = Float64.(reshape(ExtractData( filename, "/Centers/exxd"), ncx, ncz)) ε̇zz = Float64.(reshape(ExtractData( filename, "/Centers/ezzd"), ncx, ncz)) @@ -226,10 +230,6 @@ end ϕ = ExtractField(filename, "/Centers/phi", centroids, false, 0) divu = ExtractField(filename, "/Centers/divu", centroids, false, 0) - @show mean(τxx) - @show mean(τxzc) - - T_hr = zeros(size(ph_hr)); T_hr[1:2:end-1,1:2:end-1] .= T; T_hr[2:2:end-0,2:2:end-0] .= T if LAB_color ph_hr[(ph_hr.==2 .|| ph_hr.==3) .&& T_hr.BCp.type[c] != 30 && mesh->BCp.type[c] != 31 ) { mesh->p_lith[c] += model->bkg_pressure; - mesh->p_lith[c] += 0.5*model->gz * mesh->dz * rho_eff; + mesh->p_lith[c] += 0.5*model->gz * mesh->dz * rho_eff; } } } + + // --------------------------------------------- + + // Stress boundary conditions: Lithostatic pressure on lateral sides + for( l=ncz-2; l>=0; l--) { + + int kW = 0; + int kE = nx-2; + int c1W = kW + l*(nx-1); + int c1E = kE + l*(nx-1); + + // Initialise + mesh->sxx_W[l] = 0.0; + mesh->sxx_E[l] = 0.0; + + // Density + double rhoW, rhoE; + if ( mode == 0 ) rhoW = rhoE = RHO_REF; + if ( mode == 1 ) { + rhoW = 0.5*(mesh->rho_s[kW + l*(nx)] + mesh->rho_s[kW + (l+1)*(nx)]); + rhoE = 0.5*(mesh->rho_s[kE + l*(nx)] + mesh->rho_s[kE + (l+1)*(nx)]); + } + + // Initialise pressure variables : Compute lithostatic pressure + if ( mesh->BCp.type[c1W] != 30 ) mesh->sxx_W[l] = mesh->sxx_W[l+1] - model->gz * mesh->dz * rhoW; + if ( mesh->BCp.type[c1E] != 30 ) mesh->sxx_E[l] = mesh->sxx_E[l+1] - model->gz * mesh->dz * rhoE; + } + + // Add confining pressure + for( l=0; lrho_s[kW + l*(nx)] + mesh->rho_s[kW + (l+1)*(nx)]); + rhoE = 0.5*(mesh->rho_s[kE + l*(nx)] + mesh->rho_s[kE + (l+1)*(nx)]); + } + + // Correction + if ( mesh->BCp.type[c1W] != 30 ) mesh->sxx_W[l] += model->bkg_pressure; + if ( mesh->BCp.type[c1E] != 30 ) mesh->sxx_E[l] += model->bkg_pressure; + if ( mesh->BCp.type[c1W] != 30 ) mesh->sxx_W[l] += 0.5*model->gz * mesh->dz * rhoW; + if ( mesh->BCp.type[c1E] != 30 ) mesh->sxx_E[l] += 0.5*model->gz * mesh->dz * rhoE; + if ( mesh->BCp.type[c1W] != 30 ) mesh->sxx_W[l] *= -1.; + if ( mesh->BCp.type[c1E] != 30 ) mesh->sxx_E[l] *= -1.; + } + // // + eps diff --git a/MDLIB/HDF5Output.c b/MDLIB/HDF5Output.c index a93db3cf..ad469ac8 100755 --- a/MDLIB/HDF5Output.c +++ b/MDLIB/HDF5Output.c @@ -1011,8 +1011,10 @@ void WriteOutputHDF5Particles( grid *mesh, markers *particles, surface *topo, ma part_sxxd = DoodzMalloc( sizeof(float) *Nb_part_viz ); part_sxz = DoodzMalloc( sizeof(float) *Nb_part_viz ); part_index = DoodzMalloc( sizeof(int) *Nb_part_viz ); // Tracer: allocate - part_nx = DoodzMalloc( sizeof(float) *Nb_part_viz ); - part_nz = DoodzMalloc( sizeof(float) *Nb_part_viz ); + if (model.anisotropy == 1 ) { + part_nx = DoodzMalloc( sizeof(float) *Nb_part_viz ); + part_nz = DoodzMalloc( sizeof(float) *Nb_part_viz ); + } // // Tracer: store selected particle // for (k=0; kNb_part; k++) { @@ -1043,8 +1045,8 @@ void WriteOutputHDF5Particles( grid *mesh, markers *particles, surface *topo, ma part_sxz[k] = (float)particles->sxz[ind]; part_ph[k] = (char)particles->phase[ind]; part_gen[k] = (char)particles->generation[ind]; - part_nx[k] = (float)particles->nx[ind]; - part_nz[k] = (float)particles->nz[ind]; + if (model.anisotropy == 1) part_nx[k] = (float)particles->nx[ind]; + if (model.anisotropy == 1) part_nz[k] = (float)particles->nz[ind]; part_index[k] = k; ind += d_part; } @@ -1181,8 +1183,8 @@ void WriteOutputHDF5Particles( grid *mesh, markers *particles, surface *topo, ma AddFieldToGroup( FileName, "Particles", "sxxd", 'f', Nb_part_viz, part_sxxd, 1 ); AddFieldToGroup( FileName, "Particles", "sxz", 'f', Nb_part_viz, part_sxz, 1 ); - AddFieldToGroup( FileName, "Particles", "nx", 'f', Nb_part_viz, part_nx, 1 ); - AddFieldToGroup( FileName, "Particles", "nz", 'f', Nb_part_viz, part_nz, 1 ); + if (model.anisotropy == 1) AddFieldToGroup( FileName, "Particles", "nx", 'f', Nb_part_viz, part_nx, 1 ); + if (model.anisotropy == 1) AddFieldToGroup( FileName, "Particles", "nz", 'f', Nb_part_viz, part_nz, 1 ); if ( model.free_surface == 1 ) { AddFieldToGroup( FileName, "Topo", "height" , 'f', (model.Nx), Cheight, 1 ); @@ -1211,8 +1213,8 @@ void WriteOutputHDF5Particles( grid *mesh, markers *particles, surface *topo, ma DoodzFree( part_Vx ); DoodzFree( part_Vz ); DoodzFree( part_ph ); - DoodzFree( part_nx ); - DoodzFree( part_nz ); + if (model.anisotropy == 1) DoodzFree( part_nx ); + if (model.anisotropy == 1) DoodzFree( part_nz ); DoodzFree( part_gen ); DoodzFree( part_T ); DoodzFree( part_P ); diff --git a/MDLIB/Main_DOODZ.c b/MDLIB/Main_DOODZ.c index 1d198d2a..b3b3aa0e 100755 --- a/MDLIB/Main_DOODZ.c +++ b/MDLIB/Main_DOODZ.c @@ -57,14 +57,18 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) { .scaling = inputFile.scaling, .materials = inputFile.materials, .crazyConductivity = NULL, - .flux = NULL, + .flux = NULL, }; if (input.model.free_surface) { input.flux = malloc(sizeof(LateralFlux)); if (input.flux != NULL) { *input.flux = (LateralFlux){.east = -0.0e3, .west = -0.0e3}; - } + } + input.stress = malloc(sizeof(LateralFlux)); + if (input.stress != NULL) { + *input.stress = (LateralFlux){.east = -0.0e3, .west = -0.0e3}; + } } if (setup->MutateInput) { @@ -622,6 +626,8 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) { printf("Running with normal conductivity for the asthenosphere...\n"); } // Allocate and initialise solution and RHS vectors + RheologicalOperators( &mesh, &input.model, &input.materials, &input.scaling, 0, input.model.anisotropy ); // SURE? + ApplyBC( &mesh, &input.model ); // SURE? SetBCs(*setup->SetBCs, &input, &mesh, &topo); // Reset fields and BC values if needed @@ -743,12 +749,12 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) { WriteOutputHDF5Particles( &mesh, &particles, &topo, &topo_chain, &topo_ini, &topo_chain_ini, input.model, "Particles_BeforeSolve", input.materials, input.scaling ); } - // Build discrete system of equations - Linearised Picard withou Newton linearisatio nor anisotripic contributions + // Build discrete system of equations - Linearised Picard without Newton linearisation nor anisotropic contributions int kill_aniso = 0, kill_Newton = 0; RheologicalOperators( &mesh, &input.model, &input.materials, &input.scaling, kill_Newton, kill_aniso ); BuildStokesOperatorDecoupled ( &mesh, input.model, 0, mesh.p_corr, mesh.p_in, mesh.u_in, mesh.v_in, &Stokes, &StokesA, &StokesB, &StokesC, &StokesD, 1 ); - // Rheological operators including physical contrbutions from anisotropy + // Rheological operators including physical contributions from anisotropy RheologicalOperators( &mesh, &input.model, &input.materials, &input.scaling, 0, input.model.anisotropy ); // Build discrete system of equations - Jacobian (do it also for densification since drhodp is needed) diff --git a/MDLIB/MemoryAllocFree.c b/MDLIB/MemoryAllocFree.c index ab6b0a2d..e444b91a 100755 --- a/MDLIB/MemoryAllocFree.c +++ b/MDLIB/MemoryAllocFree.c @@ -516,6 +516,12 @@ grid GridAlloc(params *model) { mesh.bet_n = DoodzCalloc((Nx - 1) * (Nz - 1), sizeof(double)); mesh.bet_s = DoodzCalloc((Nx - 0) * (Nz - 0), sizeof(double)); + // Stress boundary conditions + mesh.sxx_W = DoodzCalloc((Nz - 1), sizeof(double)); + mesh.sxx_E = DoodzCalloc((Nz - 1), sizeof(double)); + mesh.szz_S = DoodzCalloc((Nx - 1), sizeof(double)); + mesh.szz_N = DoodzCalloc((Nx - 1), sizeof(double)); + // Grid indices mesh.kvx = DoodzCalloc(Nx * NzVx, sizeof(int)); mesh.lvx = DoodzCalloc(Nx * NzVx, sizeof(int)); @@ -635,6 +641,10 @@ grid GridAlloc(params *model) { if (model->anisotropy == 1) Initialise1DArrayDouble(mesh.FS_AR_s, (Nx - 0) * (Nz - 0), 1.0); if (model->anisotropy == 1) mesh.aniso_factor_n = DoodzCalloc((Nx - 1) * (Nz - 1), sizeof(double)); // To delete if (model->anisotropy == 1) mesh.aniso_factor_s = DoodzCalloc((Nx - 0) * (Nz - 0), sizeof(double)); // To delete + if (model->anisotropy == 1) { + for (int k=0; k<(Nx - 1) * (Nz - 1); k++) mesh.aniso_factor_n[k] = 1.0; + for (int k=0; k<(Nx - 0) * (Nz - 0); k++) mesh.aniso_factor_s[k] = 1.0; + } // Compressibility mesh.p0_n = DoodzCalloc((Nx - 1) * (Nz - 1), sizeof(double)); mesh.p0_s = DoodzCalloc((Nx) * (Nz), sizeof(double)); @@ -821,6 +831,12 @@ void GridFree(grid *mesh, params *model) { DoodzFree(mesh->bet_n); DoodzFree(mesh->bet_s); + // Stress boundary conditions + DoodzFree(mesh->sxx_W); + DoodzFree(mesh->sxx_E); + DoodzFree(mesh->szz_S); + DoodzFree(mesh->szz_N); + // Grid indices DoodzFree( mesh->kvx ); DoodzFree( mesh->lvx ); diff --git a/MDLIB/RheologyDensity.c b/MDLIB/RheologyDensity.c index 82a55466..b0cc6db4 100644 --- a/MDLIB/RheologyDensity.c +++ b/MDLIB/RheologyDensity.c @@ -1747,8 +1747,8 @@ void StrainRateComponents( grid* mesh, scale scaling, params* model ) { mesh->div_u[c0] = dvxdx + dvzdz + dvydy; // Normal strain rates - mesh->exxd[c0] = dvxdx - 1.0/3.0*mesh->div_u[c0]; - mesh->ezzd[c0] = dvzdz - 1.0/3.0*mesh->div_u[c0]; + mesh->exxd[c0] = dvxdx - model->compressible*1.0/3.0*mesh->div_u[c0]; // SURE? + mesh->ezzd[c0] = dvzdz - model->compressible*1.0/3.0*mesh->div_u[c0]; // SURE? } } diff --git a/MDLIB/Setup.c b/MDLIB/Setup.c index 89bab263..0af46449 100644 --- a/MDLIB/Setup.c +++ b/MDLIB/Setup.c @@ -1039,3 +1039,47 @@ bool IsRectangleCoordinates(Coordinates coordinates, Rectangle rectangle, double /*--------------------------------------------------------------------------------------------------------------------*/ /*------------------------------------------------------ M-Doodz -----------------------------------------------------*/ /*--------------------------------------------------------------------------------------------------------------------*/ + +char* DefaultTextFilename( char filename[] ) { + + int n_basename = 0; + bool start=false, end=false; + + // Go once reverse through the full file name + for (int k=strlen(filename)-1; k>=0; k--){ + // Look for slashes + if (filename[k] == '\\' || filename[k] == '/') end = true; + if (start==true && end==false) { + // printf("%c", filename[k]); + n_basename++; + } + // Look for start of file extension + if (filename[k] == '.') start = true; + } + + int len_basename = n_basename; + // char basename[len_basename]; + char *basename = calloc(len_basename, sizeof(char)); + start=false, end=false; + n_basename = 0; + + for (int k=strlen(filename)-1; k>=0; k--){ + // Look for slashes + if (filename[k] == '\\' || filename[k] == '/') end = true; + if (start==true && end==false) { + basename[len_basename-1-n_basename] = filename[k]; + n_basename++; + } + // Look for start of file extension + if (filename[k] == '.') start = true; + } + + strcat(basename, ".txt"); + printf("%s\n", basename); + + return basename; +} + +/*--------------------------------------------------------------------------------------------------------------------*/ +/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/ +/*--------------------------------------------------------------------------------------------------------------------*/ \ No newline at end of file diff --git a/MDLIB/StokesAssemblyDecoupled.c b/MDLIB/StokesAssemblyDecoupled.c index 914ce844..bfe65784 100755 --- a/MDLIB/StokesAssemblyDecoupled.c +++ b/MDLIB/StokesAssemblyDecoupled.c @@ -826,14 +826,14 @@ void Continuity_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesC, Spar double pc=0.0, uW=0.0, uE=0.0, vN=0.0, vS=0.0, oop_fact = 1.0, comp_fact = 0.0; - const double rhoW = 0.5*(mesh->rho_n[c1] + mesh->rho_n[c1-nx]); - const double rhoE = 0.5*(mesh->rho_n[c1+1] + mesh->rho_n[c1-nx+1]); - const double rhoS = 0.5*(mesh->rho_n[c1-nx]+ mesh->rho_n[c1-nx+1]); - const double rhoN = 0.5*(mesh->rho_n[c1] + mesh->rho_n[c1+1]); - const double drhodx = (rhoE-rhoW)*one_dx; - const double drhody = (rhoN-rhoS)*one_dz; - const double VxC = 0.5*(mesh->u_in[c1] + mesh->u_in[c1+1]); - const double VzC = 0.5*(mesh->v_in[c3] + mesh->v_in[c3+nxvz]); + // const double rhoW = 0.5*(mesh->rho_n[c1] + mesh->rho_n[c1-nx]); + // const double rhoE = 0.5*(mesh->rho_n[c1+1] + mesh->rho_n[c1-nx+1]); + // const double rhoS = 0.5*(mesh->rho_n[c1-nx]+ mesh->rho_n[c1-nx+1]); + // const double rhoN = 0.5*(mesh->rho_n[c1] + mesh->rho_n[c1+1]); + // const double drhodx = (rhoE-rhoW)*one_dx; + // const double drhody = (rhoN-rhoS)*one_dz; + // const double VxC = 0.5*(mesh->u_in[c1] + mesh->u_in[c1+1]); + // const double VzC = 0.5*(mesh->v_in[c3] + mesh->v_in[c3+nxvz]); const double adv_rho = 1.0; // + model.dt*VxC*drhodx/mesh->rho_n[c2] + model.dt*VzC*drhodx/mesh->rho_n[c2]; // if (fabs(adv_rho)>1+1e-10) printf("%2.2e %2.2e\n ", model.dt*VxC*drhodx/mesh->rho_n[c2] + model.dt*VzC*drhodx/mesh->rho_n[c2], drhodx); @@ -883,14 +883,20 @@ void Continuity_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesC, Spar /*------------------------------------------------------ M-Doodz -----------------------------------------------------*/ /*--------------------------------------------------------------------------------------------------------------------*/ -void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, SparseMat *StokesB, int Assemble, int lev, int stab, int comp, double om, int sign, params model, double one_dx, double one_dz, double one_dx_dx, double one_dz_dz, double one_dx_dz, double celvol, grid* mesh, int ith, int c1, int c2, int c3, int nx, int ncx, int nxvz, int eqn, double* u, double* v, double* p, int **JtempA, double **AtempA, int *nnzc2A, int **JtempB, double **AtempB, int *nnzc2B ) { +void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, SparseMat *StokesB, int Assemble, int ix, int iz, int stab, int comp, double om, int sign, params model, double one_dx, double one_dz, double one_dx_dx, double one_dz_dz, double one_dx_dz, double celvol, grid* mesh, int ith, int c1, int c2, int c3, int nx, int nz, int eqn, double* u, double* v, double* p, int **JtempA, double **AtempA, int *nnzc2A, int **JtempB, double **AtempB, int *nnzc2B ) { - double dx = mesh->dx; - double dz = mesh->dz; - double out_of_plane = 1.0; - if ( model.out_of_plane==1 ) out_of_plane = 3.0/2.0; - double uC_corr = 0.0; + const int nxvz = nx+1, ncx = nx-1; + const double dx = mesh->dx, dz = mesh->dz; + double OOP = 1.0, uC_corr = 0.0; + if ( model.out_of_plane==1 ) OOP = 3.0/2.0; + char *BCv_t = mesh->BCv.type, *BCu_t = mesh->BCu.type, *BCg_t = mesh->BCg.type, *BCp_t = mesh->BCp.type; + + // Regular or stress BC nodes? + double SxxBCW = 0., SxxBCE = 0.; + if (ix==0 && BCu_t[c1]==2) SxxBCW = 1.0; + if (ix==nx-1 && BCu_t[c1]==2) SxxBCE = 1.0; + // Stencil int iVxC = c1; int iVxW = iVxC-1; int iVxE = iVxC+1; @@ -904,48 +910,56 @@ void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars int iPrE = c2+1; int ixyN = c1; int ixyS = c1-nx; + + // Linear equation coefficient + double uS=0.0, uN=0.0, uW=0.0, uE=0.0, uC=0.0, vSW=0.0, vSE=0.0, vNW=0.0, vNE=0.0, pE=0.0, pW=0.0; // Periodic stencil - if ( mesh->BCu.type[iVxC]==-2 ) { + if ( BCu_t[iVxC]==-2 ) { iVxW = c1+nx-2; iPrW = c2+ncx; iVzSW = c3-2; iVzNW = c3+nxvz-2; } - - // double D11E = mesh->D11_n[iPrE]; - // double D11W = mesh->D11_n[iPrW]; - // double D33N = mesh->D33_s[ixyN]; - // double D33S = mesh->D33_s[ixyS]; - double D11E = 2*mesh->eta_n[iPrE]; - double D11W = 2*mesh->eta_n[iPrW]; - double D33N = mesh->eta_s[ixyN]; - double D33S = mesh->eta_s[ixyS]; - comp = 1.0; - - double uS=0.0, uN=0.0, uW=0.0, uE=0.0, uC=0.0, vSW=0.0, vSE=0.0, vNW=0.0, vNE=0.0, pE=0.0, pW=0.0; + // Variable PDE coefficient + const double D11E = mesh->D11_n[iPrE]; + const double D11W = mesh->D11_n[iPrW]; + const double D33N = mesh->D33_s[ixyN]; + const double D33S = mesh->D33_s[ixyS]; + // Flags double inE=0.0, inW=0.0; - if (mesh->BCp.type[iPrW] == -1) inW = 1.0; - if (mesh->BCp.type[iPrE] == -1) inE = 1.0; - - double inS=0.0, inN = 0.0; - if (mesh->BCg.type[ixyS] != 30 && mesh->BCu.type[iVxS] != 13) inS = 1.0; - if (mesh->BCg.type[ixyN] != 30 && mesh->BCu.type[iVxN] != 13) inN = 1.0; - - // Simpler - uW = (1.0/3.0)*D11W*inW*(comp*out_of_plane - 3)/pow(dx, 2); - uC = (1.0/3.0)*(3*pow(dx, 2)*(D33N*inN + D33S*inS) - pow(dz, 2)*(D11E*inE + D11W*inW)*(comp*out_of_plane - 3))/(pow(dx, 2)*pow(dz, 2)); - uE = (1.0/3.0)*D11E*inE*(comp*out_of_plane - 3)/pow(dx, 2); - uS = -D33S*inS/pow(dz, 2); - uN = -D33N*inN/pow(dz, 2); - vSW = ((1.0/3.0)*D11W*comp*inW*out_of_plane - D33S*inS)/(dx*dz); - vSE = (-1.0/3.0*D11E*comp*inE*out_of_plane + D33S*inS)/(dx*dz); - vNW = (-1.0/3.0*D11W*comp*inW*out_of_plane + D33N*inN)/(dx*dz); - vNE = ((1.0/3.0)*D11E*comp*inE*out_of_plane - D33N*inN)/(dx*dz); - pW = -inW*one_dx; - pE = inE*one_dx; + if (BCp_t[iPrW] == -1) inW = 1.0; + if (BCp_t[iPrE] == -1) inE = 1.0; + + double inS=1.0, inN = 1.0, RegS = 1.0, RegN = 1.0, NeuS = 0.0, NeuN = 0.0, DirS = 0.0, DirN = 0.0; + if (BCg_t[ixyS] == 30 ) inS=0.0; + if (BCg_t[ixyN] == 30 ) inN=0.0; + if (BCu_t[iVxS] == 13) {NeuS = 1.0; RegS = 0.0;} + if (BCu_t[iVxN] == 13) {NeuN = 1.0; RegN = 0.0;} + if (BCu_t[iVxS] == 11) {DirS = 1.0; RegS = 0.0;} + if (BCu_t[iVxN] == 11) {DirN = 1.0; RegN = 0.0;} + + double DirSW = 0.0, DirNW = 0.0; + if (BCv_t[c3-nxvz] == 11) DirSW = 1.0; + if (BCv_t[c3] == 11) DirNW = 1.0; + double DirSE = 0.0, DirNE = 0.0; + if (BCv_t[c3-nxvz+1] == 11) DirSE = 1.0; + if (BCv_t[c3+1] == 11) DirNE = 1.0; + + // Linear equation coefficient (AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC.ipynb) + uW = (1.0/3.0)*D11W*inW*(1 - SxxBCW)*(OOP*comp*inW - 3)/pow(dx, 2); + uC = -1.0/6.0*(SxxBCE*(2*D11W*pow(dz, 2)*inW*(OOP*comp*inW - 3) - 3*pow(dx, 2)*(D33N*inN*(2*DirN + RegN) + D33S*inS*(2*DirS + RegS))) + SxxBCW*(2*D11E*pow(dz, 2)*inE*(OOP*comp*inE - 3) - 3*pow(dx, 2)*(D33N*inN*(2*DirN + RegN) + D33S*inS*(2*DirS + RegS))) + 2*(3*pow(dx, 2)*(D33N*inN*(2*DirN + RegN) + D33S*inS*(2*DirS + RegS)) - pow(dz, 2)*(D11E*inE*(OOP*comp*inE - 3) + D11W*inW*(OOP*comp*inW - 3)))*(SxxBCE + SxxBCW - 1))/(pow(dx, 2)*pow(dz, 2)); + uE = (1.0/3.0)*D11E*inE*(1 - SxxBCE)*(OOP*comp*inE - 3)/pow(dx, 2); + uS = (1.0/2.0)*D33S*RegS*inS*(SxxBCE + SxxBCW - 2)/pow(dz, 2); + uN = (1.0/2.0)*D33N*RegN*inN*(SxxBCE + SxxBCW - 2)/pow(dz, 2); + vSW = (1.0/6.0)*(-3*D33S*SxxBCW*inS*(2*DirSE*SxxBCE - SxxBCE - SxxBCW + 1) + SxxBCE*(2*D11W*OOP*comp*pow(inW, 2) - 3*D33S*inS*(2*DirSE*SxxBCE - SxxBCE - SxxBCW + 1)) - 2*(D11W*OOP*comp*pow(inW, 2) - 3*D33S*inS*(2*DirSE*SxxBCE - SxxBCE - SxxBCW + 1))*(SxxBCE + SxxBCW - 1))/(dx*dz); + vSE = (1.0/6.0)*(3*D33S*SxxBCE*inS*(2*DirSW*SxxBCW - SxxBCE - SxxBCW + 1) - SxxBCW*(2*D11E*OOP*comp*pow(inE, 2) - 3*D33S*inS*(2*DirSW*SxxBCW - SxxBCE - SxxBCW + 1)) + 2*(D11E*OOP*comp*pow(inE, 2) - 3*D33S*inS*(2*DirSW*SxxBCW - SxxBCE - SxxBCW + 1))*(SxxBCE + SxxBCW - 1))/(dx*dz); + vNW = (1.0/6.0)*(3*D33N*SxxBCW*inN*(2*DirNE*SxxBCE - SxxBCE - SxxBCW + 1) - SxxBCE*(2*D11W*OOP*comp*pow(inW, 2) - 3*D33N*inN*(2*DirNE*SxxBCE - SxxBCE - SxxBCW + 1)) + 2*(D11W*OOP*comp*pow(inW, 2) - 3*D33N*inN*(2*DirNE*SxxBCE - SxxBCE - SxxBCW + 1))*(SxxBCE + SxxBCW - 1))/(dx*dz); + vNE = (1.0/6.0)*(-3*D33N*SxxBCE*inN*(2*DirNW*SxxBCW - SxxBCE - SxxBCW + 1) + SxxBCW*(2*D11E*OOP*comp*pow(inE, 2) - 3*D33N*inN*(2*DirNW*SxxBCW - SxxBCE - SxxBCW + 1)) - 2*(D11E*OOP*comp*pow(inE, 2) - 3*D33N*inN*(2*DirNW*SxxBCW - SxxBCE - SxxBCW + 1))*(SxxBCE + SxxBCW - 1))/(dx*dz); + pW = (SxxBCW - 1)/dx; + pE = (1 - SxxBCE)/dx; // Stabilisation with density gradients if ( stab==1 ) { @@ -956,102 +970,158 @@ void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars uC += uC_corr; } - // Add contribution from non-conforming Dirichlets - if ( mesh->BCu.type[iVxS] == 11 ) uC -= uS; - if ( mesh->BCu.type[iVxN] == 11 ) uC -= uN; - if ( Assemble == 1 ) { + const bool SetVxS = BCu_t[iVxS] != 30; + const bool SetVxW = BCu_t[iVxW] != 30 && (ix>0 || BCu_t[iVxC]==-2); + const bool SetVxC = BCu_t[iVxC] != 30; + const bool SetVxE = BCu_t[iVxE] != 30 && (ix0 || BCu_t[iVxC]==-2); + const bool SetVzSE = BCv_t[iVzSE] != 30 && (ix0 || BCu_t[iVxC]==-2); + const bool SetVzNE = BCv_t[iVzNE] != 30 && (ix0 || BCu_t[iVxC]==-2); + const bool SetPrE = BCp_t[iPrE] != 30 && (ixb[eqn] *= celvol; StokesB->b[eqn] *= celvol; //-------------------- // dsxx/dx - normal stencil - if (mesh->BCu.type[iVxS] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxS], &(nnzc2A[ith]), uS*celvol, mesh->BCu.type[iVxS], mesh->BCu.val[iVxS], StokesA->bbc); - if (mesh->BCu.type[iVxW] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxW], &(nnzc2A[ith]), uW*celvol, mesh->BCu.type[iVxW], mesh->BCu.val[iVxW], StokesA->bbc); - AddCoeff3( JtempA[ith], AtempA[ith], eqn, eqn, &(nnzc2A[ith]), uC*celvol, mesh->BCu.type[iVxC], mesh->BCu.val[iVxC], StokesA->bbc); - if (mesh->BCu.type[iVxE] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxE], &(nnzc2A[ith]), uE*celvol, mesh->BCu.type[iVxE], mesh->BCu.val[iVxE], StokesA->bbc); - if (mesh->BCu.type[iVxN] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxN], &(nnzc2A[ith]), uN*celvol, mesh->BCu.type[iVxN], mesh->BCu.val[iVxN], StokesA->bbc); + if (SetVxS) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxS], &(nnzc2A[ith]), uS*celvol, BCu_t[iVxS], mesh->BCu.val[iVxS], StokesA->bbc); + if (SetVxW) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxW], &(nnzc2A[ith]), uW*celvol, BCu_t[iVxW], mesh->BCu.val[iVxW], StokesA->bbc); + if (SetVxC) AddCoeff3( JtempA[ith], AtempA[ith], eqn, eqn, &(nnzc2A[ith]), uC*celvol, BCu_t[iVxC], mesh->BCu.val[iVxC], StokesA->bbc); + if (SetVxE) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxE], &(nnzc2A[ith]), uE*celvol, BCu_t[iVxE], mesh->BCu.val[iVxE], StokesA->bbc); + if (SetVxN) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxN], &(nnzc2A[ith]), uN*celvol, BCu_t[iVxN], mesh->BCu.val[iVxN], StokesA->bbc); //-------------------- - if ( mesh->BCv.type[iVzSW] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzSW], &(nnzc2A[ith]), vSW*celvol, mesh->BCv.type[iVzSW], mesh->BCv.val[iVzSW], StokesA->bbc); - if ( mesh->BCv.type[iVzSE] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzSE], &(nnzc2A[ith]), vSE*celvol, mesh->BCv.type[iVzSE], mesh->BCv.val[iVzSE], StokesA->bbc); - if ( mesh->BCv.type[iVzNW] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzNW], &(nnzc2A[ith]), vNW*celvol, mesh->BCv.type[iVzNW], mesh->BCv.val[iVzNW], StokesA->bbc); - if ( mesh->BCv.type[iVzNE] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzNE], &(nnzc2A[ith]), vNE*celvol, mesh->BCv.type[iVzNE], mesh->BCv.val[iVzNE], StokesA->bbc); + if (SetVzSW) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzSW], &(nnzc2A[ith]), vSW*celvol, BCv_t[iVzSW], mesh->BCv.val[iVzSW], StokesA->bbc); + if (SetVzSE) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzSE], &(nnzc2A[ith]), vSE*celvol, BCv_t[iVzSE], mesh->BCv.val[iVzSE], StokesA->bbc); + if (SetVzNW) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzNW], &(nnzc2A[ith]), vNW*celvol, BCv_t[iVzNW], mesh->BCv.val[iVzNW], StokesA->bbc); + if (SetVzNE) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzNE], &(nnzc2A[ith]), vNE*celvol, BCv_t[iVzNE], mesh->BCv.val[iVzNE], StokesA->bbc); //-------------------- - if ( mesh->BCp.type[iPrE] != 30 && mesh->BCp.type[iPrW] != 30 ) { - AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrW] - Stokes->neq_mom, &(nnzc2B[ith]), pW*celvol, mesh->BCp.type[iPrW], mesh->BCp.val[iPrW], StokesB->bbc); - AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrE] - Stokes->neq_mom, &(nnzc2B[ith]), pE*celvol, mesh->BCp.type[iPrE], mesh->BCp.val[iPrE], StokesB->bbc); - } + // if ( BCp_t[iPrE] != 30 && BCp_t[iPrW] != 30 ) { + if (SetPrW) AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrW] - Stokes->neq_mom, &(nnzc2B[ith]), pW*celvol, BCp_t[iPrW], mesh->BCp.val[iPrW], StokesB->bbc); + if (SetPrE) AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrE] - Stokes->neq_mom, &(nnzc2B[ith]), pE*celvol, BCp_t[iPrE], mesh->BCp.val[iPrE], StokesB->bbc); + // } } else { // Residual function - StokesA->F[eqn] = 0.0; - // StokesA->F[eqn] += (inE*mesh->sxxd[iPrE] - inW*mesh->sxxd[iPrW] )/dx; - // StokesA->F[eqn] -= (inE*mesh->p_corr[iPrE] - inW*mesh->p_corr[iPrW])/dx; - // StokesA->F[eqn] += (inN*mesh->sxz[ixyN] - inS*mesh->sxz[ixyS] )/dz; - StokesA->F[eqn] += (mesh->sxxd[iPrE] - mesh->sxxd[iPrW] )/dx; - StokesA->F[eqn] -= (mesh->p_corr[iPrE] - mesh->p_corr[iPrW])/dx; - StokesA->F[eqn] += (mesh->sxz[ixyN] - mesh->sxz[ixyS] )/dz; - StokesA->F[eqn] *= -1.0; - StokesA->F[eqn] -= StokesA->b[eqn] - uC_corr*u[iVxC]; // no body force - StokesA->F[eqn] *= celvol; + double F_Reg = 0., F_W =0., F_E =0., SxxE, SxxW; + + // Residual: regular nodes + F_Reg += (mesh->sxxd[iPrE] - mesh->sxxd[iPrW] )/dx; + F_Reg -= (mesh->p_corr[iPrE] - mesh->p_corr[iPrW])/dx; + F_Reg += (mesh->sxz[ixyN] - mesh->sxz[ixyS] )/dz; + + // Residual: stress BC west + SxxE = mesh->sxxd[iPrE] - mesh->p_corr[iPrE]; + SxxW = 2.0*mesh->BCu.val[iVxC] - SxxE; + F_W += (SxxE - SxxW)/dx; + F_W += (mesh->sxz[ixyN] - mesh->sxz[ixyS])/dz; + F_W *= 0.5; + + // Residual: stress BC east + SxxW = mesh->sxxd[iPrW] - mesh->p_corr[iPrW]; + SxxE = 2.0*mesh->BCu.val[iVxC] - SxxW; + F_E += (SxxE - SxxW)/dx; + F_E += (mesh->sxz[ixyN] - mesh->sxz[ixyS])/dz; + F_E *= 0.5; + + // Residual + StokesA->F[eqn] = (1.0 - SxxBCE - SxxBCW)*F_Reg + SxxBCW*F_W + SxxBCE*F_E; + StokesA->F[eqn] += StokesA->b[eqn] - uC_corr*u[iVxC]; // no body force + StokesA->F[eqn] *= -celvol; } -} +} /*--------------------------------------------------------------------------------------------------------------------*/ /*------------------------------------------------------ M-Doodz -----------------------------------------------------*/ /*--------------------------------------------------------------------------------------------------------------------*/ -void Zmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, SparseMat *StokesB, int Assemble, int lev, int stab, int comp, double om, int sign, params model, double one_dx, double one_dz, double one_dx_dx, double one_dz_dz, double one_dx_dz, double celvol, grid* mesh, int ith, int c1, int c2, int c3, int nx, int ncx, int nxvz, int eqn, double* u, double* v, double* p, int **JtempA, double **AtempA, int *nnzc2A, int **JtempB, double **AtempB, int *nnzc2B ) { +void Zmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, SparseMat *StokesB, int Assemble, int ix, int iz, int stab, int comp, double om, int sign, params model, double one_dx, double one_dz, double one_dx_dx, double one_dz_dz, double one_dx_dz, double celvol, grid* mesh, int ith, int c1, int c2, int c3, int nx, int nz, int eqn, double* u, double* v, double* p, int **JtempA, double **AtempA, int *nnzc2A, int **JtempB, double **AtempB, int *nnzc2B ) { - double dx = mesh->dx; - double dz = mesh->dz; - double out_of_plane = 1.0; - if ( model.out_of_plane==1 ) out_of_plane = 3.0/2.0; - double vC_corr = 0.0; + const int nxvz = nx+1, ncx = nx-1; + const double dx = mesh->dx, dz = mesh->dz; + double OOP = 1.0, vC_corr = 0.0; + if ( model.out_of_plane==1 ) OOP = 3.0/2.0; + char *BCv_t = mesh->BCv.type, *BCu_t = mesh->BCu.type, *BCg_t = mesh->BCg.type, *BCp_t = mesh->BCp.type; + + // Regular or stress BC nodes? + double SzzBCS = 0., SzzBCN = 0.; + if (iz==0 && BCv_t[c3]==2) SzzBCS = 1.0; + if (iz==nz-1 && BCv_t[c3]==2) SzzBCN = 1.0; + // Stencil int iVzC = c3; int iVzW = iVzC-1; int iVzE = iVzC+1; int iVzS = iVzC-nxvz; int iVzN = iVzC+nxvz; + // ------------------ int iVxNW = c1+nx-1; int iVxNE = c1+nx; int iVxSW = c1-1; int iVxSE = c1; + // ------------------ int iPrS = c2; int iPrN = c2+ncx; int ixyE = c1; int ixyW = c1-1; - - double D22S = mesh->D22_n[iPrS]; - double D22N = mesh->D22_n[iPrN]; - double D33E = mesh->D33_s[ixyE]; - double D33W = mesh->D33_s[ixyW]; - comp = 1.0; - + + if ( BCv_t[iVzW] == -12 && ix==1 ) iVzW = c3+nxvz-3; + if ( BCv_t[iVzE] == -12 && ix==nxvz-2 ) iVzE = c3-nxvz+3; + + // Linear equation coefficient double vS=0.0, vN=0.0, vW=0.0, vE=0.0, vC=0.0, uSW=0.0, uSE=0.0, uNW=0.0, uNE=0.0, pN=0.0, pS=0.0; + // Variable PDE coefficient + const double D22S = mesh->D22_n[iPrS]; + const double D22N = mesh->D22_n[iPrN]; + const double D33E = mesh->D33_s[ixyE]; + const double D33W = mesh->D33_s[ixyW]; + // Flags double inN=0.0, inS=0.0; - if (mesh->BCp.type[iPrS] == -1) inS = 1.0; - if (mesh->BCp.type[iPrN] == -1) inN = 1.0; + if (BCp_t[iPrS] == -1) inS = 1.0; + if (BCp_t[iPrN] == -1) inN = 1.0; - double inW=0.0, inE = 0.0; - if (mesh->BCg.type[ixyW] != 30 && mesh->BCv.type[iVzW] != 13) inW = 1.0; - if (mesh->BCg.type[ixyE] != 30 && mesh->BCv.type[iVzE] != 13) inE = 1.0; + double inW=1.0, inE = 1.0, RegW = 1., RegE = 1., NeuW = 0., NeuE = 0, DirW = 0., DirE = 0.; + if (BCg_t[ixyW] == 30) inW = 0.0; + if (BCg_t[ixyE] == 30) inE = 0.0; + if (BCv_t[iVzW] == 13) {NeuW = 1.0; RegW = 0.0;} + if (BCv_t[iVzE] == 13) {NeuE = 1.0; RegE = 0.0;} + if (BCv_t[iVzW] == 11) {DirW = 1.0; RegW = 0.0;} + if (BCv_t[iVzE] == 11) {DirE = 1.0; RegE = 0.0;} + + double DirSW =0., DirNW = 0. ; + if (mesh->BCu.type[c1-1] == 11) DirSW = 1.0; + if (mesh->BCu.type[c1+nx-1] == 11) DirNW = 1.0; + + double DirSE = 0., DirNE = 0.; + if (mesh->BCu.type[c1] == 11) DirSE = 1.0; + if (mesh->BCu.type[c1+nx] == 11) DirNE = 1.0; // Simpler - vW = -D33W*inW/pow(dx, 2); - vC = (1.0/3.0)*(-pow(dx, 2)*(D22N*inN + D22S*inS)*(comp*out_of_plane - 3) + 3*pow(dz, 2)*(D33E*inE + D33W*inW))/(pow(dx, 2)*pow(dz, 2)); - vE = -D33E*inE/pow(dx, 2); - vS = (1.0/3.0)*D22S*inS*(comp*out_of_plane - 3)/pow(dz, 2); - vN = (1.0/3.0)*D22N*inN*(comp*out_of_plane - 3)/pow(dz, 2); - uSW = ((1.0/3.0)*D22S*comp*inS*out_of_plane - D33W*inW)/(dx*dz); - uSE = (-1.0/3.0*D22S*comp*inS*out_of_plane + D33E*inE)/(dx*dz); - uNW = (-1.0/3.0*D22N*comp*inN*out_of_plane + D33W*inW)/(dx*dz); - uNE = ((1.0/3.0)*D22N*comp*inN*out_of_plane - D33E*inE)/(dx*dz); - pS = -inS*one_dz; - pN = inN*one_dz; - + vW = (1.0/2.0)*D33W*RegW*inW*(SzzBCN + SzzBCS - 2)/pow(dx, 2); + vC = (1.0/6.0)*(-SzzBCN*(2*D22S*pow(dx, 2)*(OOP*comp - 3) - 3*pow(dz, 2)*(D33E*inE*(2*DirE + RegE) + D33W*inW*(2*DirW + RegW))) - SzzBCS*(2*D22N*pow(dx, 2)*(OOP*comp - 3) - 3*pow(dz, 2)*(D33E*inE*(2*DirE + RegE) + D33W*inW*(2*DirW + RegW))) + 2*(pow(dx, 2)*(D22N + D22S)*(OOP*comp - 3) - 3*pow(dz, 2)*(D33E*inE*(2*DirE + RegE) + D33W*inW*(2*DirW + RegW)))*(SzzBCN + SzzBCS - 1))/(pow(dx, 2)*pow(dz, 2)); + vE = (1.0/2.0)*D33E*RegE*inE*(SzzBCN + SzzBCS - 2)/pow(dx, 2); + vS = (1.0/3.0)*D22S*(1 - SzzBCS)*(OOP*comp - 3)/pow(dz, 2); + vN = (1.0/3.0)*D22N*(1 - SzzBCN)*(OOP*comp - 3)/pow(dz, 2); + uSW = (1.0/6.0)*(-3*D33W*SzzBCS*inW*(2*DirNW*SzzBCN - SzzBCN - SzzBCS + 1) + SzzBCN*(2*D22S*OOP*comp - 3*D33W*inW*(2*DirNW*SzzBCN - SzzBCN - SzzBCS + 1)) - 2*(D22S*OOP*comp - 3*D33W*inW*(2*DirNW*SzzBCN - SzzBCN - SzzBCS + 1))*(SzzBCN + SzzBCS - 1))/(dx*dz); + uSE = (1.0/6.0)*(3*D33E*SzzBCS*inE*(2*DirNE*SzzBCN - SzzBCN - SzzBCS + 1) - SzzBCN*(2*D22S*OOP*comp - 3*D33E*inE*(2*DirNE*SzzBCN - SzzBCN - SzzBCS + 1)) + 2*(D22S*OOP*comp - 3*D33E*inE*(2*DirNE*SzzBCN - SzzBCN - SzzBCS + 1))*(SzzBCN + SzzBCS - 1))/(dx*dz); + uNW = (1.0/6.0)*(3*D33W*SzzBCN*inW*(2*DirSW*SzzBCS - SzzBCN - SzzBCS + 1) - SzzBCS*(2*D22N*OOP*comp - 3*D33W*inW*(2*DirSW*SzzBCS - SzzBCN - SzzBCS + 1)) + 2*(D22N*OOP*comp - 3*D33W*inW*(2*DirSW*SzzBCS - SzzBCN - SzzBCS + 1))*(SzzBCN + SzzBCS - 1))/(dx*dz); + uNE = (1.0/6.0)*(-3*D33E*SzzBCN*inE*(2*DirSE*SzzBCS - SzzBCN - SzzBCS + 1) + SzzBCS*(2*D22N*OOP*comp - 3*D33E*inE*(2*DirSE*SzzBCS - SzzBCN - SzzBCS + 1)) - 2*(D22N*OOP*comp - 3*D33E*inE*(2*DirSE*SzzBCS - SzzBCN - SzzBCS + 1))*(SzzBCN + SzzBCS - 1))/(dx*dz); + pS = (SzzBCS - 1)/dz; + pN = (1 - SzzBCN)/dz; + + // if (eqn==100) printf("%d %d %d\n", eqn, mesh->BCu.type[iVxSW], mesh->BCu.type[iVxNW]); + + // if (mesh->BCu.type[iVxSW] == 2 && BCv_t[iVzW] == 13) uSW = 0.; + // if (mesh->BCu.type[iVxNW] == 2 && BCv_t[iVzW] == 13) uNW = 0.; + // if (mesh->BCu.type[iVxSW] == 2) uSE = 0.; + // if (mesh->BCu.type[iVxNW] == 2) uNE = 0.; + // Stabilisation with density gradients if ( stab==1 ) { double drhodz = (mesh->rho_n[c2+ncx] - mesh->rho_n[c2])*one_dz; @@ -1060,88 +1130,77 @@ void Zmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars if ((vC+vC_corr)<0.0 || (vC+vC_corr)BCv.type[iVzW] == 11 ) vC -= vW; - if ( mesh->BCv.type[iVzE] == 11 ) vC -= vE; if ( Assemble == 1 ) { - + const bool SetVzS = BCv_t[iVzS] != 30 && iz>0; + const bool SetVzW = BCv_t[iVzW] != 30 ; + const bool SetVzC = BCv_t[iVzC] != 30; + const bool SetVzE = BCv_t[iVzE] != 30 ; + const bool SetVzN = BCv_t[iVzN] != 30 && iz0 ; + const bool SetVxSE = BCu_t[iVxSE] != 30 && iz>0 ; + const bool SetVxNW = BCu_t[iVxNW] != 30 && iz0; + const bool SetPrN = BCp_t[iPrN] != 30 && izb[eqn] *= celvol; StokesB->b[eqn] *= celvol; //-------------------- - if ( mesh->BCu.type[iVxSW] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxSW], &(nnzc2A[ith]), uSW*celvol, mesh->BCu.type[iVxSW], mesh->BCu.val[iVxSW], StokesA->bbc ); - if ( mesh->BCu.type[iVxSE] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxSE], &(nnzc2A[ith]), uSE*celvol, mesh->BCu.type[iVxSE], mesh->BCu.val[iVxSE], StokesA->bbc ); - if ( mesh->BCu.type[iVxNW] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxNW], &(nnzc2A[ith]), uNW*celvol, mesh->BCu.type[iVxNW], mesh->BCu.val[iVxNW], StokesA->bbc ); - if ( mesh->BCu.type[iVxNE] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxNE], &(nnzc2A[ith]), uNE*celvol, mesh->BCu.type[iVxNE], mesh->BCu.val[iVxNE], StokesA->bbc ); + if (SetVxSW) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxSW], &(nnzc2A[ith]), uSW*celvol, BCu_t[iVxSW], mesh->BCu.val[iVxSW], StokesA->bbc ); + if (SetVxSE) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxSE], &(nnzc2A[ith]), uSE*celvol, BCu_t[iVxSE], mesh->BCu.val[iVxSE], StokesA->bbc ); + if (SetVxNW) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxNW], &(nnzc2A[ith]), uNW*celvol, BCu_t[iVxNW], mesh->BCu.val[iVxNW], StokesA->bbc ); + if (SetVxNE) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxNE], &(nnzc2A[ith]), uNE*celvol, BCu_t[iVxNE], mesh->BCu.val[iVxNE], StokesA->bbc ); //-------------------- - if ( mesh->BCv.type[iVzS] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzS], &(nnzc2A[ith]), vS*celvol, mesh->BCv.type[iVzS], mesh->BCv.val[iVzS], StokesA->bbc ); - if ( mesh->BCv.type[iVzW] != 30 && mesh->BCv.type[iVzW] != -12 ) { - AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzW], &(nnzc2A[ith]), vW*celvol, mesh->BCv.type[iVzW], mesh->BCv.val[iVzW], StokesA->bbc ); - } - // Periodic - if ( mesh->BCv.type[iVzW] == -12 ) { - AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3+nxvz-3], &(nnzc2A[ith]), vW*celvol, mesh->BCv.type[c3+nxvz-3], mesh->BCv.val[c3+nxvz-3], StokesA->bbc ); - } - AddCoeff3( JtempA[ith], AtempA[ith], eqn, eqn, &(nnzc2A[ith]), vC*celvol, mesh->BCv.type[iVzC], mesh->BCv.val[iVzC], StokesA->bbc ); - if ( mesh->BCv.type[iVzE] != 30 && mesh->BCv.type[iVzE] != -12 ) { - AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzE], &(nnzc2A[ith]), vE*celvol, mesh->BCv.type[iVzE], mesh->BCv.val[iVzE], StokesA->bbc ); - } - // Periodic - if ( mesh->BCv.type[iVzE] == -12 ) { - AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3-nxvz+3], &(nnzc2A[ith]), vE*celvol, mesh->BCv.type[c3-nxvz+3], mesh->BCv.val[c3-nxvz+3], StokesA->bbc ); - } - if ( mesh->BCv.type[iVzN] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzN], &(nnzc2A[ith]), vN*celvol, mesh->BCv.type[iVzN], mesh->BCv.val[iVzN], StokesA->bbc ); + if (SetVzS) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzS], &(nnzc2A[ith]), vS*celvol, BCv_t[iVzS], mesh->BCv.val[iVzS], StokesA->bbc ); + if (SetVzW) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzW], &(nnzc2A[ith]), vW*celvol, BCv_t[iVzW], mesh->BCv.val[iVzW], StokesA->bbc ); + if (SetVzC) AddCoeff3( JtempA[ith], AtempA[ith], eqn, eqn, &(nnzc2A[ith]), vC*celvol, BCv_t[iVzC], mesh->BCv.val[iVzC], StokesA->bbc ); + if (SetVzE) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzE], &(nnzc2A[ith]), vE*celvol, BCv_t[iVzE], mesh->BCv.val[iVzE], StokesA->bbc ); + if (SetVzN) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzN], &(nnzc2A[ith]), vN*celvol, BCv_t[iVzN], mesh->BCv.val[iVzN], StokesA->bbc ); //-------------------- - if ( mesh->BCp.type[iPrS] != 30 && mesh->BCp.type[iPrN] != 30) { - AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrS] - Stokes->neq_mom, &(nnzc2B[ith]), pS*celvol, mesh->BCp.type[iPrS], mesh->BCp.val[iPrS], StokesB->bbc ); - AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrN] - Stokes->neq_mom, &(nnzc2B[ith]), pN*celvol, mesh->BCp.type[iPrN], mesh->BCp.val[iPrN], StokesB->bbc ); - } - - // StokesA->b[eqn] *= celvol; - // StokesB->b[eqn] *= celvol; - //-------------------- - - // if ( mesh->BCu.type[c1+nx] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+nx], &(nnzc2A[ith]), uNE*celvol, mesh->BCu.type[c1+nx], mesh->BCu.val[c1+nx], StokesA->bbc ); - // if ( mesh->BCu.type[c1] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1], &(nnzc2A[ith]), uSE*celvol, mesh->BCu.type[c1], mesh->BCu.val[c1], StokesA->bbc ); - // if ( mesh->BCu.type[c1+nx-1] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+nx-1], &(nnzc2A[ith]), uNW*celvol, mesh->BCu.type[c1+nx-1], mesh->BCu.val[c1+nx-1], StokesA->bbc ); - // if ( mesh->BCu.type[c1-1] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1-1], &(nnzc2A[ith]), uSW*celvol, mesh->BCu.type[c1-1], mesh->BCu.val[c1-1], StokesA->bbc ); - //-------------------- - // if ( mesh->BCv.type[c3-nxvz] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3-nxvz], &(nnzc2A[ith]), vS*celvol, mesh->BCv.type[c3-nxvz], mesh->BCv.val[c3-nxvz], StokesA->bbc ); - // if ( mesh->BCv.type[c3-1] != 30 && mesh->BCv.type[c3-1] != -12 ) { - // AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3-1], &(nnzc2A[ith]), vW*celvol, mesh->BCv.type[c3-1], mesh->BCv.val[c3-1], StokesA->bbc ); - // } - // // Periodic - // if ( mesh->BCv.type[c3-1] == -12 ) { - // AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3+nxvz-3], &(nnzc2A[ith]), vW*celvol, mesh->BCv.type[c3+nxvz-3], mesh->BCv.val[c3+nxvz-3], StokesA->bbc ); - // } - // AddCoeff3( JtempA[ith], AtempA[ith], eqn, eqn, &(nnzc2A[ith]), vC*celvol, mesh->BCv.type[c3], mesh->BCv.val[c3], StokesA->bbc ); - // if ( mesh->BCv.type[c3+1] != 30 && mesh->BCv.type[c3+1] != -12 ) { - // AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3+1], &(nnzc2A[ith]), vE*celvol, mesh->BCv.type[c3+1], mesh->BCv.val[c3+1], StokesA->bbc ); - // } - // // Periodic - // if ( mesh->BCv.type[c3+1] == -12 ) { - // AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3-nxvz+3], &(nnzc2A[ith]), vE*celvol, mesh->BCv.type[c3-nxvz+3], mesh->BCv.val[c3-nxvz+3], StokesA->bbc ); - // } - // if ( mesh->BCv.type[c3+nxvz] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3+nxvz], &(nnzc2A[ith]), vN*celvol, mesh->BCv.type[c3+nxvz], mesh->BCv.val[c3+nxvz], StokesA->bbc ); - // //-------------------- - // if ( mesh->BCp.type[c2] != 30 && mesh->BCp.type[c2+ncx] != 30) { - // AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[c2] - Stokes->neq_mom, &(nnzc2B[ith]), pS*celvol, mesh->BCp.type[c2], mesh->BCp.val[c2], StokesB->bbc ); - // AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[c2+ncx] - Stokes->neq_mom, &(nnzc2B[ith]), pN*celvol, mesh->BCp.type[c2+ncx], mesh->BCp.val[c2+ncx], StokesB->bbc ); + // if ( BCp_t[iPrS] != 30 && BCp_t[iPrN] != 30) { + if (SetPrS) AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrS] - Stokes->neq_mom, &(nnzc2B[ith]), pS*celvol, BCp_t[iPrS], mesh->BCp.val[iPrS], StokesB->bbc ); + if (SetPrN) AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrN] - Stokes->neq_mom, &(nnzc2B[ith]), pN*celvol, BCp_t[iPrN], mesh->BCp.val[iPrN], StokesB->bbc ); // } } else { + // Residual function - StokesA->F[eqn] = 0.0; - // StokesA->F[eqn] += (inN*mesh->szzd[iPrN] - inS*mesh->szzd[iPrS])/dz; - // StokesA->F[eqn] -= (inN*mesh->p_corr[iPrN] - inS*mesh->p_corr[iPrS])/dz; - // StokesA->F[eqn] += (inE*mesh->sxz[ixyE] - inW*mesh->sxz[ixyW]) /dx; - StokesA->F[eqn] += (mesh->szzd[iPrN] - mesh->szzd[iPrS])/dz; - StokesA->F[eqn] -= (mesh->p_corr[iPrN] - mesh->p_corr[iPrS])/dz; - StokesA->F[eqn] += (mesh->sxz[ixyE] - mesh->sxz[ixyW]) /dx; - StokesA->F[eqn] *= -1.0; - StokesA->F[eqn] -= StokesA->b[eqn] - vC_corr*v[iVzC]; - StokesA->F[eqn] *= celvol; + double F_Reg = 0., F_S =0., F_N =0., SzzS, SzzN; + + // Residual: regular nodes + F_Reg += (mesh->szzd[iPrN] - mesh->szzd[iPrS])/dz; + F_Reg -= (mesh->p_corr[iPrN] - mesh->p_corr[iPrS])/dz; + F_Reg += (mesh->sxz[ixyE] - mesh->sxz[ixyW]) /dx; + + // Residual: stress BC south + SzzN = mesh->szzd[iPrN] - mesh->p_corr[iPrN]; + SzzS = 2.0*mesh->BCv.val[iVzC] - SzzN; + F_S += (SzzN - SzzS)/dz; + F_S += (mesh->sxz[ixyE] - mesh->sxz[ixyW]) /dx; + F_S *= 0.5; + + // Residual: stress BC north + SzzS = mesh->szzd[iPrS] - mesh->p_corr[iPrS]; + SzzN = 2.0*mesh->BCv.val[iVzC] - SzzS; + F_N += (SzzN - SzzS)/dz; + F_N += (mesh->sxz[ixyE] - mesh->sxz[ixyW]) /dx; + F_N *= 0.5; + + // Residual + StokesA->F[eqn] = (1.0 - SzzBCS - SzzBCN)*F_Reg + SzzBCS*F_S + SzzBCN*F_N; + StokesA->F[eqn] += StokesA->b[eqn] - vC_corr*v[iVzC]; + StokesA->F[eqn] *= -celvol; + + // StokesA->F[eqn] = 0.0; + // StokesA->F[eqn] += (mesh->szzd[iPrN] - mesh->szzd[iPrS])/dz; + // StokesA->F[eqn] -= (mesh->p_corr[iPrN] - mesh->p_corr[iPrS])/dz; + // StokesA->F[eqn] += (mesh->sxz[ixyE] - mesh->sxz[ixyW]) /dx; + // StokesA->F[eqn] *= -1.0; + // StokesA->F[eqn] -= StokesA->b[eqn] - vC_corr*v[iVzC]; + // StokesA->F[eqn] *= celvol; } } @@ -1314,9 +1373,7 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ if (model.free_surface_stab>0) stab = 1; if (model.compressible==1 ) comp = 1; double theta = model.free_surface_stab; - - // if ( stab==1 ) printf("It tastes like Bo'\n"); - + // Decompose domain int n_th, N, ith; int *DD, *estart, *eend, *last_eqn; @@ -1339,10 +1396,11 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ n_th = omp_get_num_threads(); } #pragma omp barrier + + if (Assemble==1) printf("Assemble Picard operator on %d threads\n", n_th); // Matrix initialisation if ( Assemble == 1 ) { - nnzA = 9*((mesh->Nx-1) * mesh->Nz + (mesh->Nz-1) * mesh->Nx); nnzB = 5*((mesh->Nx-1) * (mesh->Nz-1)); nnzC = nnzB; @@ -1356,7 +1414,6 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ AllocMat( StokesB, nnzB ); AllocMat( StokesC, nnzC ); AllocMat( StokesD, nnzD ); - } // Build velocity block RHS @@ -1431,18 +1488,14 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ } //--------------------- INNER NODES ---------------------// - if ( l>0 && l0 && kBCu.type[c1]==-2 ) { - Xmomentum_InnerNodesDecoupled( Stokes, StokesA, StokesB, Assemble, lev, stab, comp, theta, sign, model, one_dx, one_dz, one_dx_dx, one_dz_dz, one_dx_dz, celvol, mesh, ith, c1, c2, c3, nx, ncx, nxvz, eqn, u, v, p_corr, JtempA, AtempA, nnzc2A, JtempB, AtempB, nnzc2B ); + if ( k>=0 && k<=nx-1 && l>0 && l0 && kBCu.type[c1] != 2 ?? + Xmomentum_InnerNodesDecoupled( Stokes, StokesA, StokesB, Assemble, k, l, stab, comp, theta, sign, model, one_dx, one_dz, one_dx_dx, one_dz_dz, one_dx_dz, celvol, mesh, ith, c1, c2, c3, nx, nz, eqn, u, v, p_corr, JtempA, AtempA, nnzc2A, JtempB, AtempB, nnzc2B ); } - } + } } } - + //-----------------------------------------------------------------// if ( Assemble == 1 ) { @@ -1493,10 +1546,10 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ } //--------------------- INNER NODES ---------------------// - if ( k>0 && k0 && l0 && k=0 && l<=nz-1 ) { + // Maybe one should avoid sides (l>0 && lBCv.type[c3] != 2 ?? + Zmomentum_InnerNodesDecoupled( Stokes, StokesA, StokesB, Assemble, k, l, stab, comp, theta, sign, model, one_dx, one_dz, one_dx_dx, one_dz_dz, one_dx_dz, celvol, mesh, ith, c1, c2, c3, nx, nz, eqn, u, v, p_corr, JtempA, AtempA, nnzc2A, JtempB, AtempB, nnzc2B ); + } } } } @@ -1528,10 +1581,9 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ AllocateTempMatArraysDecoupled( &AtempC, &ItempC, &JtempC, n_th, nnzC, Stokes->neq_cont, DD, &nnzc2C ); AllocateTempMatArraysDecoupled( &AtempD, &ItempD, &JtempD, n_th, nnzD, Stokes->neq_cont, DD, &nnzc2D ); } - + #pragma omp parallel shared( eend, estart, mesh, Stokes, StokesC, StokesD, u, v, p, nnzc2C, AtempC, JtempC, ItempC, nnzc2D, AtempD, JtempD, ItempD, last_eqn ) private( ith, l, k, c1, c2, c3, eqn, comp ) firstprivate( model, Assemble, lev, one_dx_dx, one_dz_dz, one_dx_dz, one_dx, one_dz, sign, theta, stab, celvol, nx, ncx, nxvz, nzvx ) - { - + { ith = omp_get_thread_num(); for( c2=estart[ith]; c2b[k] = StokesC->b[k]; } - - // Final index StokesA->Ic[StokesA->neq] = nnzcA; StokesB->Ic[StokesB->neq] = nnzcB; @@ -1709,7 +1757,6 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ //--------------------------------------// -#ifndef _VG_ if ( model.writer_debug == 1 ) { char *filename; @@ -1718,33 +1765,33 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ // Fill in DD data structure OutputSparseMatrix OutputDDA, OutputDDB, OutputDDC, OutputDDD; - - OutputDDA.V = StokesA->A; - OutputDDA.Ic = StokesA->Ic; - OutputDDA.J = StokesA->J; - OutputDDA.b = StokesA->b; - OutputDDA.F = StokesA->F; - OutputDDB.V = StokesB->A; - OutputDDB.Ic = StokesB->Ic; - OutputDDB.J = StokesB->J; - OutputDDB.b = StokesB->b; - OutputDDC.V = StokesC->A; - OutputDDC.Ic = StokesC->Ic; - OutputDDC.J = StokesC->J; - OutputDDC.b = StokesC->b; - OutputDDC.F = StokesC->F; - OutputDDD.V = StokesD->A; - OutputDDD.Ic = StokesD->Ic; - OutputDDD.J = StokesD->J; - OutputDDD.b = StokesD->b; - OutputDDA.eta_cell = mesh->eta_n; + + OutputDDA.V = StokesA->A; + OutputDDA.Ic = StokesA->Ic; + OutputDDA.J = StokesA->J; + OutputDDA.b = StokesA->b; + OutputDDA.F = StokesA->F; + OutputDDB.V = StokesB->A; + OutputDDB.Ic = StokesB->Ic; + OutputDDB.J = StokesB->J; + OutputDDB.b = StokesB->b; + OutputDDC.V = StokesC->A; + OutputDDC.Ic = StokesC->Ic; + OutputDDC.J = StokesC->J; + OutputDDC.b = StokesC->b; + OutputDDC.F = StokesC->F; + OutputDDD.V = StokesD->A; + OutputDDD.Ic = StokesD->Ic; + OutputDDD.J = StokesD->J; + OutputDDD.b = StokesD->b; + OutputDDA.eta_cell = mesh->eta_n; OutputDDA.params[0] = nx; OutputDDA.params[1] = nz; OutputDDA.params[2] = mesh->dx; OutputDDA.params[3] = mesh->dz; - OutputDDA.eqn_u = DoodzMalloc(nx*nzvx*sizeof(int)); - OutputDDA.eqn_v = DoodzMalloc(nxvz*nz*sizeof(int)); - OutputDDA.eqn_p = DoodzMalloc(ncx*ncz*sizeof(int)); + OutputDDA.eqn_u = DoodzMalloc(nx*nzvx*sizeof(int)); + OutputDDA.eqn_v = DoodzMalloc(nxvz*nz*sizeof(int)); + OutputDDA.eqn_p = DoodzMalloc(ncx*ncz*sizeof(int)); for( l=0; lneq = Stokes->neq_mom; StokesB->neq = Stokes->neq_mom; @@ -1879,7 +1927,6 @@ void BuildJacobianOperatorDecoupled( grid *mesh, params model, int lev, double * AllocMat( StokesB, nnzB ); AllocMat( StokesC, nnzC ); AllocMat( StokesD, nnzD ); - } // Build velocity block RHS @@ -1960,6 +2007,15 @@ void BuildJacobianOperatorDecoupled( grid *mesh, params model, int lev, double * if ( k==0 && mesh->BCu.type[c1]==-2 ) { Xjacobian_InnerNodesDecoupled3( Stokes, StokesA, StokesB, Assemble, lev, stab, comp, theta, sign, model, one_dx, one_dz, one_dx_dx, one_dz_dz, one_dx_dz, celvol, mesh, ith, c1, c2, c3, nx, ncx, nxvz, eqn, u, v, p_corr, JtempA, AtempA, nnzc2A, JtempB, AtempB, nnzc2B, k, l ); } + + // // NOT CODED YET + // if ( k==0 && mesh->BCu.type[c1]==2 ) { + // Xmomentum_WestNeumannDecoupled( Stokes, StokesA, StokesB, Assemble, lev, stab, comp, theta, sign, model, one_dx, one_dz, one_dx_dx, one_dz_dz, one_dx_dz, celvol, mesh, ith, c1, c2, c3, nx, ncx, nxvz, eqn, u, v, p_corr, JtempA, AtempA, nnzc2A, JtempB, AtempB, nnzc2B ); + // } + + // if ( k==nx-1 && mesh->BCu.type[c1]==2 ) { + // Xmomentum_EastNeumannDecoupled( Stokes, StokesA, StokesB, Assemble, lev, stab, comp, theta, sign, model, one_dx, one_dz, one_dx_dx, one_dz_dz, one_dx_dz, celvol, mesh, ith, c1, c2, c3, nx, ncx, nxvz, eqn, u, v, p_corr, JtempA, AtempA, nnzc2A, JtempB, AtempB, nnzc2B ); + // } } } } @@ -2017,6 +2073,24 @@ void BuildJacobianOperatorDecoupled( grid *mesh, params model, int lev, double * if ( k>0 && k0 && l0 && k0 && lBCv.type[c3] == 2 ) { + // Zmomentum_SouthNeumannDecoupled( Stokes, StokesA, StokesB, Assemble, lev, stab, comp, theta, sign, model, one_dx, one_dz, one_dx_dx, one_dz_dz, one_dx_dz, celvol, mesh, ith, c1, c2, c3, nx, ncx, nxvz, eqn, u, v, p_corr, JtempA, AtempA, nnzc2A, JtempB, AtempB, nnzc2B ); + // } + + // //--------------------- INNER NODES ---------------------// + + // if ( l==nz-1 && mesh->BCv.type[c3] == 2 ) { + // Zmomentum_NorthNeumannDecoupled( Stokes, StokesA, StokesB, Assemble, lev, stab, comp, theta, sign, model, one_dx, one_dz, one_dx_dx, one_dz_dz, one_dx_dz, celvol, mesh, ith, c1, c2, c3, nx, ncx, nxvz, eqn, u, v, p_corr, JtempA, AtempA, nnzc2A, JtempB, AtempB, nnzc2B ); + // } + } } } diff --git a/MDLIB/include/mdoodz.h b/MDLIB/include/mdoodz.h index 089965c5..1833d654 100644 --- a/MDLIB/include/mdoodz.h +++ b/MDLIB/include/mdoodz.h @@ -321,7 +321,7 @@ struct MdoodzInput { params model; mat_prop materials; scale scaling; - LateralFlux *flux; + LateralFlux *flux, *stress; CrazyConductivity *crazyConductivity; Geometry *geometry; }; @@ -355,5 +355,7 @@ typedef struct { bool IsEllipseCoordinates(Coordinates coordinates, Ellipse ellipse, double scalingL); bool IsRectangleCoordinates(Coordinates coordinates, Rectangle rectangle, double scalingL); +char *DefaultTextFilename(char*); + #endif diff --git a/MDLIB/mdoodz-private.h b/MDLIB/mdoodz-private.h index 1b5f2e17..833a053c 100755 --- a/MDLIB/mdoodz-private.h +++ b/MDLIB/mdoodz-private.h @@ -129,6 +129,7 @@ typedef struct { double *kc_x, *kc_z; double *FreeSurfW_s, *FreeSurfW_n; double *noise_n, *noise_s; + double *sxx_W, *sxx_E, *szz_S, *szz_N; } grid; diff --git a/SETS/ShearBoxLaeti.c b/SETS/ShearBoxLaeti.c index 4316203f..38802796 100755 --- a/SETS/ShearBoxLaeti.c +++ b/SETS/ShearBoxLaeti.c @@ -36,6 +36,168 @@ double SetTxz(MdoodzInput *input, Coordinates coordinates, int phase) { return input->model.preload_sxz; } + +SetBC SetBCVx(MdoodzInput *instance, POSITION position, Coordinates coord) { + SetBC bc; + const double radius = instance->model.user1 / instance->scaling.L; + double x, z; + // ----------- Pure shear ----------- // + if (instance->model.shear_style==0) { + if (position == W) { + x = coord.x; + z = coord.z; + bc.type = 2; + bc.value = -1.;//instance->stress->west*0.0;//1*z*4; + } else if (position == E) { + x = coord.x; + z = coord.z; + bc.type = 2; + bc.value = -1.;//instance->stress->east*0.0;//1*z*4; + } else if (position == S || position == SE || position == SW) { + x = coord.x; + z = coord.z + instance->model.dx / 2.0;// make sure it lives on the boundary + bc.type = 13; + bc.value = 0*3e-3; + // Pick a point in the middle + if ( (fabs(x)-0.0) < instance->model.dx/2) { + printf("COUCOU x S\n"); + bc.type = 11; + bc.value = 0*3e-3; + } + } else if (position == N || position == NE || position == NW) { + x = coord.x; + z = coord.z - instance->model.dx / 2.0;// make sure it lives on the boundary + bc.type = 13; + bc.value = 0.0; + // Pick a point in the middle + if ( (fabs(x)-0.0) < instance->model.dx/2) { + printf("COUCOU x N\n"); + bc.type = 11; + bc.value = 0*3e-3; + } + } else { + bc.type = -1; + bc.value = 0.0; + } + } + // ----------- Simple shear ----------- // + else { + const double Lz = (double) (instance->model.zmax - instance->model.zmin); + if (position == S || position == SE || position == SW) { + bc.type = 11; + bc.value = -instance->model.bkg_strain_rate * Lz; + } else if (position == N || position == NE || position == NW) { + bc.type = 11; + bc.value = instance->model.bkg_strain_rate * Lz; + } else if (position == E) { + // bc.type = -12; + // bc.value = 0.0; + bc.type = 0; + bc.value = 0.0; + } else if (position == W) { + bc.type = 0; + bc.value = 0.0; + // bc.type = -2; + // bc.value = 0.0; + } else { + bc.type = -1; + bc.value = 0.0; + } + } + return bc; +} + +SetBC SetBCVz(MdoodzInput *instance, POSITION position, Coordinates coord) { + SetBC bc; + double x, z; + double stress = instance->model.user3 / instance->scaling.S; + // ----------- Pure shear ----------- // + if (instance->model.shear_style==0) { + if (position == N) { + x = coord.x; + z = coord.z; + bc.type = 2; + bc.value = 1.; + } + else if (position == S) { + x = coord.x; + z = coord.z; + bc.type = 2; + bc.value = 1.; + } + else if (position == W || position == SE || position == SW) { + x = coord.x + instance->model.dx / 2.0; + z = coord.z; + bc.type = 13; + bc.value = 0.0; + // Pick a point in the middle + if ( (fabs(z)-0.0) < instance->model.dz/2) { + printf("COUCOU z W\n"); + bc.type = 11; + bc.value = 0; + } + } + else if (position == E || position == NE || position == NW) { + x = coord.x - instance->model.dx / 2.0; + z = coord.z; + bc.type = 13; + bc.value = 0.0; + // Pick a point in the middle + if ( (fabs(z)-0.0) < instance->model.dz/2) { + printf("COUCOU z E\n"); + bc.type = 11; + bc.value = 0; + } + } else { + bc.type = -1; + bc.value = 0.0; + } + } + // ----------- Simple shear ----------- // + else { + if ( position == W || position == NW || position == SW) { + // bc.type = -12; + // bc.value = 0.0; + bc.type = 13; + bc.value = 0.0; + // // Pick a point in the middle + // if ( (fabs(coord.z)-0.0) < instance->model.dz/2) { + // printf("COUCOU z W\n"); + // bc.type = 11; + // bc.value = 0; + // } + } + else if (position == E || position == NE || position == SE) { + // bc.type = -12; + // bc.value = 0.0; + bc.type = 13; + bc.value = 0.0; + } else if (position == S) { + bc.type = 2; + bc.value = stress; + // if ( (fabs(coord.x)-0.0) < instance->model.dx/2) { + // bc.type = 2; + // bc.value = 1.0; + // } + } + else if (position == N) { + bc.type = 0; + bc.value = 0.0; + // if ( (fabs(coord.x)-0.0) < instance->model.dx/2) { + // bc.type = 2; + // bc.value = 1.0; + // } + + } else { + bc.type = -1; + bc.value = 0.0; + } + } + return bc; +} + + + int main(int nargs, char *args[]) { // Input file name char *input_file; @@ -56,8 +218,8 @@ int main(int nargs, char *args[]) { .SetTxz = SetTxz, }, .SetBCs = &(SetBCs_ff){ - .SetBCVx = SetPureOrSimpleShearBCVx, - .SetBCVz = SetPureOrSimpleShearBCVz, + .SetBCVx = SetBCVx, + .SetBCVz = SetBCVz, }, }; RunMDOODZ(input_file, &setup); diff --git a/SETS/ShearTemplate.c b/SETS/ShearTemplate.c index 64c8e469..049476f4 100755 --- a/SETS/ShearTemplate.c +++ b/SETS/ShearTemplate.c @@ -24,7 +24,7 @@ int main(int nargs, char *args[]) { // Input file name char *input_file; if ( nargs < 2 ) { - asprintf(&input_file, "ShearTemplate.txt"); // Default + asprintf(&input_file, DefaultTextFilename(__FILE__)); // Default } else { asprintf(&input_file, "%s", args[1]); // Custom diff --git a/SETS/ShearTemplateCI/LinearPureshearAnisotropic.txt b/SETS/ShearTemplateCI/LinearPureshearAnisotropic.txt new file mode 100755 index 00000000..5ee62985 --- /dev/null +++ b/SETS/ShearTemplateCI/LinearPureshearAnisotropic.txt @@ -0,0 +1,136 @@ +/**** RESTART ****/ +istep = 000001 +irestart = 0 + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 1 +writer_markers = 0 +writer_debug = 0 +writer_energies = 0 +gnuplot_log_res = 0 +writer_subfolder = LinearPureshearAnisotropic + +/**** SCALES ****/ +eta = 1 +L = 1 +V = 1 +T = 1 + +/**** SPACE ****/ +Nx = 21 +Nz = 21 +xmin = -5.000000e-1 +zmin = -5.000000e-1 +xmax = 5.000000e-1 +zmax = 5.000000e-1 + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 + +/**** TIME ****/ +Nt = 1 +constant_dt = 1 +Courant = 0.3 +dt = 0.5 +RK = 4 + +/**** LINEAR SOLVER ****/ +noisy = 1 +penalty = 1e5 +lin_solver = -1 +diag_scaling = 0 +lin_abs_div = 1e-9 +lin_rel_div = 1e-9 +lin_abs_mom = 1e-11 +lin_rel_mom = 1e-11 + +/**** NON-LINEAR SOLVER ****/ +Newton = 1 +line_search = 1 +nit_max = 10 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-9 +nonlin_rel_mom = 1e-9 +nonlin_abs_div = 1e-9 +nonlin_rel_div = 1e-9 + +/**** VISCOSITY CUT-OFF ****/ +min_eta = 1e-3 +max_eta = 1e6 + +/**** SWITCHES ****/ +anisotropy = 1 +elastic = 0 +free_surface = 0 +free_surface_stab = 0 +finite_strain = 1 +advection = 0 + +/**** SETUP-DEPENDANT ****/ +shear_style = 0 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1 +user0 = 1 / temperature [°C] +user1 = 0.05 / inclusion radius [m] +user2 = 0 +user3 = 0 + +/**** GRAVITY ****/ +gx = 0.0000 +gz = -0e0 + +/**** MAT PROPERTIES ****/ +Nb_phases = 2 + +/**** PHASE 1 ****/ +ID = 0 +rho = 2700.00 / matrix +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1 +npwl = 3.0 +Qpwl = 0 +aniso_factor = 5.0 / THIS IS DEPRECATED USE LINE BELOW: +ani_fac_v = 3.0 +ani_fac_e = 1.0 +aniso_angle = 30.0 + +/**** PHASE 2 ****/ +ID = 1 +rho = 2750.00 / inclusion +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e1 +npwl = 0 +Qpwl = 0 diff --git a/SETS/ShearTemplateCI/LinearPureshearIsotropic.txt b/SETS/ShearTemplateCI/LinearPureshearIsotropic.txt new file mode 100755 index 00000000..304c9a96 --- /dev/null +++ b/SETS/ShearTemplateCI/LinearPureshearIsotropic.txt @@ -0,0 +1,136 @@ +/**** RESTART ****/ +istep = 000001 +irestart = 0 + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 1 +writer_markers = 0 +writer_debug = 0 +writer_energies = 0 +gnuplot_log_res = 0 +writer_subfolder = LinearPureshearIsotropic + +/**** SCALES ****/ +eta = 1 +L = 1 +V = 1 +T = 1 + +/**** SPACE ****/ +Nx = 21 +Nz = 21 +xmin = -5.000000e-1 +zmin = -5.000000e-1 +xmax = 5.000000e-1 +zmax = 5.000000e-1 + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 + +/**** TIME ****/ +Nt = 1 +constant_dt = 1 +Courant = 0.3 +dt = 0.5 +RK = 4 + +/**** LINEAR SOLVER ****/ +noisy = 1 +penalty = 1e5 +lin_solver = -1 +diag_scaling = 0 +lin_abs_div = 1e-9 +lin_rel_div = 1e-9 +lin_abs_mom = 1e-11 +lin_rel_mom = 1e-11 + +/**** NON-LINEAR SOLVER ****/ +Newton = 1 +line_search = 1 +nit_max = 10 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-9 +nonlin_rel_mom = 1e-9 +nonlin_abs_div = 1e-9 +nonlin_rel_div = 1e-9 + +/**** VISCOSITY CUT-OFF ****/ +min_eta = 1e-3 +max_eta = 1e6 + +/**** SWITCHES ****/ +anisotropy = 0 +elastic = 0 +free_surface = 0 +free_surface_stab = 0 +finite_strain = 1 +advection = 0 + +/**** SETUP-DEPENDANT ****/ +shear_style = 0 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1 +user0 = 1 / temperature [°C] +user1 = 0.05 / inclusion radius [m] +user2 = 0 +user3 = 0 + +/**** GRAVITY ****/ +gx = 0.0000 +gz = -0e0 + +/**** MAT PROPERTIES ****/ +Nb_phases = 2 + +/**** PHASE 1 ****/ +ID = 0 +rho = 2700.00 / matrix +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1 +npwl = 3.0 +Qpwl = 0 +aniso_factor = 5.0 / THIS IS DEPRECATED USE LINE BELOW: +ani_fac_v = 3.0 +ani_fac_e = 1.0 +aniso_angle = 30.0 + +/**** PHASE 2 ****/ +ID = 1 +rho = 2750.00 / inclusion +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e1 +npwl = 0 +Qpwl = 0 diff --git a/SETS/ShearTemplateCI/LinearSimpleshearAnisotropic.txt b/SETS/ShearTemplateCI/LinearSimpleshearAnisotropic.txt new file mode 100755 index 00000000..c1cd84a8 --- /dev/null +++ b/SETS/ShearTemplateCI/LinearSimpleshearAnisotropic.txt @@ -0,0 +1,136 @@ +/**** RESTART ****/ +istep = 000001 +irestart = 0 + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 1 +writer_markers = 0 +writer_debug = 0 +writer_energies = 0 +gnuplot_log_res = 0 +writer_subfolder = LinearSimpleshearAnisotropic + +/**** SCALES ****/ +eta = 1 +L = 1 +V = 1 +T = 1 + +/**** SPACE ****/ +Nx = 21 +Nz = 21 +xmin = -5.000000e-1 +zmin = -5.000000e-1 +xmax = 5.000000e-1 +zmax = 5.000000e-1 + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 + +/**** TIME ****/ +Nt = 1 +constant_dt = 1 +Courant = 0.3 +dt = 0.5 +RK = 4 + +/**** LINEAR SOLVER ****/ +noisy = 1 +penalty = 1e5 +lin_solver = -1 +diag_scaling = 0 +lin_abs_div = 1e-9 +lin_rel_div = 1e-9 +lin_abs_mom = 1e-11 +lin_rel_mom = 1e-11 + +/**** NON-LINEAR SOLVER ****/ +Newton = 1 +line_search = 1 +nit_max = 10 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-9 +nonlin_rel_mom = 1e-9 +nonlin_abs_div = 1e-9 +nonlin_rel_div = 1e-9 + +/**** VISCOSITY CUT-OFF ****/ +min_eta = 1e-3 +max_eta = 1e6 + +/**** SWITCHES ****/ +anisotropy = 1 +elastic = 0 +free_surface = 0 +free_surface_stab = 0 +finite_strain = 1 +advection = 0 + +/**** SETUP-DEPENDANT ****/ +shear_style = 1 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1 +user0 = 1 / temperature [°C] +user1 = 0.05 / inclusion radius [m] +user2 = 0 +user3 = 0 + +/**** GRAVITY ****/ +gx = 0.0000 +gz = -0e0 + +/**** MAT PROPERTIES ****/ +Nb_phases = 2 + +/**** PHASE 1 ****/ +ID = 0 +rho = 2700.00 / matrix +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1 +npwl = 3.0 +Qpwl = 0 +aniso_factor = 5.0 / THIS IS DEPRECATED USE LINE BELOW: +ani_fac_v = 3.0 +ani_fac_e = 1.0 +aniso_angle = 30.0 + +/**** PHASE 2 ****/ +ID = 1 +rho = 2750.00 / inclusion +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e1 +npwl = 0 +Qpwl = 0 diff --git a/SETS/ShearTemplateCI/LinearSimpleshearIsotropic.txt b/SETS/ShearTemplateCI/LinearSimpleshearIsotropic.txt new file mode 100755 index 00000000..ad234b75 --- /dev/null +++ b/SETS/ShearTemplateCI/LinearSimpleshearIsotropic.txt @@ -0,0 +1,136 @@ +/**** RESTART ****/ +istep = 000001 +irestart = 0 + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 1 +writer_markers = 0 +writer_debug = 0 +writer_energies = 0 +gnuplot_log_res = 0 +writer_subfolder = LinearSimpleshearIsotropic + +/**** SCALES ****/ +eta = 1 +L = 1 +V = 1 +T = 1 + +/**** SPACE ****/ +Nx = 21 +Nz = 21 +xmin = -5.000000e-1 +zmin = -5.000000e-1 +xmax = 5.000000e-1 +zmax = 5.000000e-1 + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 + +/**** TIME ****/ +Nt = 1 +constant_dt = 1 +Courant = 0.3 +dt = 0.5 +RK = 4 + +/**** LINEAR SOLVER ****/ +noisy = 1 +penalty = 1e5 +lin_solver = -1 +diag_scaling = 0 +lin_abs_div = 1e-9 +lin_rel_div = 1e-9 +lin_abs_mom = 1e-11 +lin_rel_mom = 1e-11 + +/**** NON-LINEAR SOLVER ****/ +Newton = 1 +line_search = 1 +nit_max = 10 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-9 +nonlin_rel_mom = 1e-9 +nonlin_abs_div = 1e-9 +nonlin_rel_div = 1e-9 + +/**** VISCOSITY CUT-OFF ****/ +min_eta = 1e-3 +max_eta = 1e6 + +/**** SWITCHES ****/ +anisotropy = 0 +elastic = 0 +free_surface = 0 +free_surface_stab = 0 +finite_strain = 1 +advection = 0 + +/**** SETUP-DEPENDANT ****/ +shear_style = 1 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1 +user0 = 1 / temperature [°C] +user1 = 0.05 / inclusion radius [m] +user2 = 0 +user3 = 0 + +/**** GRAVITY ****/ +gx = 0.0000 +gz = -0e0 + +/**** MAT PROPERTIES ****/ +Nb_phases = 2 + +/**** PHASE 1 ****/ +ID = 0 +rho = 2700.00 / matrix +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1 +npwl = 3.0 +Qpwl = 0 +aniso_factor = 5.0 / THIS IS DEPRECATED USE LINE BELOW: +ani_fac_v = 3.0 +ani_fac_e = 1.0 +aniso_angle = 30.0 + +/**** PHASE 2 ****/ +ID = 1 +rho = 2750.00 / inclusion +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e1 +npwl = 0 +Qpwl = 0 diff --git a/SETS/ShearTemplateCI/NonLinearPureshearAnisotropic.txt b/SETS/ShearTemplateCI/NonLinearPureshearAnisotropic.txt new file mode 100755 index 00000000..730c9d5c --- /dev/null +++ b/SETS/ShearTemplateCI/NonLinearPureshearAnisotropic.txt @@ -0,0 +1,136 @@ +/**** RESTART ****/ +istep = 000001 +irestart = 0 + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 1 +writer_markers = 0 +writer_debug = 0 +writer_energies = 0 +gnuplot_log_res = 0 +writer_subfolder = NonLinearPureshearAnisotropic + +/**** SCALES ****/ +eta = 1 +L = 1 +V = 1 +T = 1 + +/**** SPACE ****/ +Nx = 21 +Nz = 21 +xmin = -5.000000e-1 +zmin = -5.000000e-1 +xmax = 5.000000e-1 +zmax = 5.000000e-1 + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 + +/**** TIME ****/ +Nt = 1 +constant_dt = 1 +Courant = 0.3 +dt = 0.5 +RK = 4 + +/**** LINEAR SOLVER ****/ +noisy = 1 +penalty = 1e5 +lin_solver = -1 +diag_scaling = 0 +lin_abs_div = 1e-9 +lin_rel_div = 1e-9 +lin_abs_mom = 1e-11 +lin_rel_mom = 1e-11 + +/**** NON-LINEAR SOLVER ****/ +Newton = 1 +line_search = 1 +nit_max = 10 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-9 +nonlin_rel_mom = 1e-9 +nonlin_abs_div = 1e-9 +nonlin_rel_div = 1e-9 + +/**** VISCOSITY CUT-OFF ****/ +min_eta = 1e-3 +max_eta = 1e6 + +/**** SWITCHES ****/ +anisotropy = 1 +elastic = 0 +free_surface = 0 +free_surface_stab = 0 +finite_strain = 1 +advection = 0 + +/**** SETUP-DEPENDANT ****/ +shear_style = 0 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1 +user0 = 1 / temperature [°C] +user1 = 0.05 / inclusion radius [m] +user2 = 0 +user3 = 0 + +/**** GRAVITY ****/ +gx = 0.0000 +gz = -0e0 + +/**** MAT PROPERTIES ****/ +Nb_phases = 2 + +/**** PHASE 1 ****/ +ID = 0 +rho = 2700.00 / matrix +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 0 / constant visc law +pwlv = 1 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1 +npwl = 3.0 +Qpwl = 0 +aniso_factor = 5.0 / THIS IS DEPRECATED USE LINE BELOW: +ani_fac_v = 3.0 +ani_fac_e = 1.0 +aniso_angle = 30.0 + +/**** PHASE 2 ****/ +ID = 1 +rho = 2750.00 / inclusion +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e1 +npwl = 0 +Qpwl = 0 diff --git a/SETS/ShearTemplateCI/NonLinearPureshearIsotropic.txt b/SETS/ShearTemplateCI/NonLinearPureshearIsotropic.txt new file mode 100755 index 00000000..fd4c5f7b --- /dev/null +++ b/SETS/ShearTemplateCI/NonLinearPureshearIsotropic.txt @@ -0,0 +1,136 @@ +/**** RESTART ****/ +istep = 000001 +irestart = 0 + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 1 +writer_markers = 0 +writer_debug = 0 +writer_energies = 0 +gnuplot_log_res = 0 +writer_subfolder = NonLinearPureshearIsotropic + +/**** SCALES ****/ +eta = 1 +L = 1 +V = 1 +T = 1 + +/**** SPACE ****/ +Nx = 21 +Nz = 21 +xmin = -5.000000e-1 +zmin = -5.000000e-1 +xmax = 5.000000e-1 +zmax = 5.000000e-1 + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 + +/**** TIME ****/ +Nt = 1 +constant_dt = 1 +Courant = 0.3 +dt = 0.5 +RK = 4 + +/**** LINEAR SOLVER ****/ +noisy = 1 +penalty = 1e5 +lin_solver = -1 +diag_scaling = 0 +lin_abs_div = 1e-9 +lin_rel_div = 1e-9 +lin_abs_mom = 1e-11 +lin_rel_mom = 1e-11 + +/**** NON-LINEAR SOLVER ****/ +Newton = 1 +line_search = 1 +nit_max = 10 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-9 +nonlin_rel_mom = 1e-9 +nonlin_abs_div = 1e-9 +nonlin_rel_div = 1e-9 + +/**** VISCOSITY CUT-OFF ****/ +min_eta = 1e-3 +max_eta = 1e6 + +/**** SWITCHES ****/ +anisotropy = 0 +elastic = 0 +free_surface = 0 +free_surface_stab = 0 +finite_strain = 1 +advection = 0 + +/**** SETUP-DEPENDANT ****/ +shear_style = 0 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1 +user0 = 1 / temperature [°C] +user1 = 0.05 / inclusion radius [m] +user2 = 0 +user3 = 0 + +/**** GRAVITY ****/ +gx = 0.0000 +gz = -0e0 + +/**** MAT PROPERTIES ****/ +Nb_phases = 2 + +/**** PHASE 1 ****/ +ID = 0 +rho = 2700.00 / matrix +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 0 / constant visc law +pwlv = 1 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1 +npwl = 3.0 +Qpwl = 0 +aniso_factor = 5.0 / THIS IS DEPRECATED USE LINE BELOW: +ani_fac_v = 3.0 +ani_fac_e = 1.0 +aniso_angle = 30.0 + +/**** PHASE 2 ****/ +ID = 1 +rho = 2750.00 / inclusion +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e1 +npwl = 0 +Qpwl = 0 diff --git a/SETS/ShearTemplateCI/NonLinearSimpleshearAnisotropi.txt b/SETS/ShearTemplateCI/NonLinearSimpleshearAnisotropi.txt new file mode 100755 index 00000000..20dc9b28 --- /dev/null +++ b/SETS/ShearTemplateCI/NonLinearSimpleshearAnisotropi.txt @@ -0,0 +1,136 @@ +/**** RESTART ****/ +istep = 000001 +irestart = 0 + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 1 +writer_markers = 0 +writer_debug = 0 +writer_energies = 0 +gnuplot_log_res = 0 +writer_subfolder = NonLinearSimpleshearAnisotropi + +/**** SCALES ****/ +eta = 1 +L = 1 +V = 1 +T = 1 + +/**** SPACE ****/ +Nx = 21 +Nz = 21 +xmin = -5.000000e-1 +zmin = -5.000000e-1 +xmax = 5.000000e-1 +zmax = 5.000000e-1 + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 + +/**** TIME ****/ +Nt = 1 +constant_dt = 1 +Courant = 0.3 +dt = 0.5 +RK = 4 + +/**** LINEAR SOLVER ****/ +noisy = 1 +penalty = 1e5 +lin_solver = -1 +diag_scaling = 0 +lin_abs_div = 1e-9 +lin_rel_div = 1e-9 +lin_abs_mom = 1e-11 +lin_rel_mom = 1e-11 + +/**** NON-LINEAR SOLVER ****/ +Newton = 1 +line_search = 1 +nit_max = 10 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-9 +nonlin_rel_mom = 1e-9 +nonlin_abs_div = 1e-9 +nonlin_rel_div = 1e-9 + +/**** VISCOSITY CUT-OFF ****/ +min_eta = 1e-3 +max_eta = 1e6 + +/**** SWITCHES ****/ +anisotropy = 1 +elastic = 0 +free_surface = 0 +free_surface_stab = 0 +finite_strain = 1 +advection = 0 + +/**** SETUP-DEPENDANT ****/ +shear_style = 1 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1 +user0 = 1 / temperature [°C] +user1 = 0.05 / inclusion radius [m] +user2 = 0 +user3 = 0 + +/**** GRAVITY ****/ +gx = 0.0000 +gz = -0e0 + +/**** MAT PROPERTIES ****/ +Nb_phases = 2 + +/**** PHASE 1 ****/ +ID = 0 +rho = 2700.00 / matrix +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 0 / constant visc law +pwlv = 1 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1 +npwl = 3.0 +Qpwl = 0 +aniso_factor = 5.0 / THIS IS DEPRECATED USE LINE BELOW: +ani_fac_v = 3.0 +ani_fac_e = 1.0 +aniso_angle = 30.0 + +/**** PHASE 2 ****/ +ID = 1 +rho = 2750.00 / inclusion +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e1 +npwl = 0 +Qpwl = 0 diff --git a/SETS/ShearTemplateCI/NonLinearSimpleshearIsotropic.txt b/SETS/ShearTemplateCI/NonLinearSimpleshearIsotropic.txt new file mode 100755 index 00000000..48e9b318 --- /dev/null +++ b/SETS/ShearTemplateCI/NonLinearSimpleshearIsotropic.txt @@ -0,0 +1,136 @@ +/**** RESTART ****/ +istep = 000001 +irestart = 0 + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 1 +writer_markers = 0 +writer_debug = 0 +writer_energies = 0 +gnuplot_log_res = 0 +writer_subfolder = NonLinearSimpleshearIsotropic + +/**** SCALES ****/ +eta = 1 +L = 1 +V = 1 +T = 1 + +/**** SPACE ****/ +Nx = 21 +Nz = 21 +xmin = -5.000000e-1 +zmin = -5.000000e-1 +xmax = 5.000000e-1 +zmax = 5.000000e-1 + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 + +/**** TIME ****/ +Nt = 1 +constant_dt = 1 +Courant = 0.3 +dt = 0.5 +RK = 4 + +/**** LINEAR SOLVER ****/ +noisy = 1 +penalty = 1e5 +lin_solver = -1 +diag_scaling = 0 +lin_abs_div = 1e-9 +lin_rel_div = 1e-9 +lin_abs_mom = 1e-11 +lin_rel_mom = 1e-11 + +/**** NON-LINEAR SOLVER ****/ +Newton = 1 +line_search = 1 +nit_max = 10 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-9 +nonlin_rel_mom = 1e-9 +nonlin_abs_div = 1e-9 +nonlin_rel_div = 1e-9 + +/**** VISCOSITY CUT-OFF ****/ +min_eta = 1e-3 +max_eta = 1e6 + +/**** SWITCHES ****/ +anisotropy = 0 +elastic = 0 +free_surface = 0 +free_surface_stab = 0 +finite_strain = 1 +advection = 0 + +/**** SETUP-DEPENDANT ****/ +shear_style = 1 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1 +user0 = 1 / temperature [°C] +user1 = 0.05 / inclusion radius [m] +user2 = 0 +user3 = 0 + +/**** GRAVITY ****/ +gx = 0.0000 +gz = -0e0 + +/**** MAT PROPERTIES ****/ +Nb_phases = 2 + +/**** PHASE 1 ****/ +ID = 0 +rho = 2700.00 / matrix +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 0 / constant visc law +pwlv = 1 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1 +npwl = 3.0 +Qpwl = 0 +aniso_factor = 5.0 / THIS IS DEPRECATED USE LINE BELOW: +ani_fac_v = 3.0 +ani_fac_e = 1.0 +aniso_angle = 30.0 + +/**** PHASE 2 ****/ +ID = 1 +rho = 2750.00 / inclusion +G = 1 +Cp = 1050 +k = 2.5 +Qr = 0 +C = 1e90 +phi = 30 +Slim = 1e100 +alp = 10.0e-6 +bet = 1e-11 +drho = 0 +cstv = 1 / constant visc law +pwlv = 0 / disloc. creep +linv = 0 / diff. creep +gbsv = 0 / grain boundary sliding +expv = 0 / peierls creep +gsel = 0 / grain size evo. +eta0 = 1e1 +npwl = 0 +Qpwl = 0 diff --git a/SETS/StressBC_FreeSurf.c b/SETS/StressBC_FreeSurf.c new file mode 100755 index 00000000..863ece16 --- /dev/null +++ b/SETS/StressBC_FreeSurf.c @@ -0,0 +1,138 @@ +#include "mdoodz.h" +#include "stdio.h" +#include "stdlib.h" +#include "math.h" + +double SetSurfaceZCoord(MdoodzInput *instance, double x_coord) { + const double TopoLevel = .4 / instance->scaling.L; + const double h_pert = instance->model.user3 / instance->scaling.L; + return TopoLevel + 0.0*x_coord;// + h_pert * (3330.0 - 2800.0) / 2800.0 * cos(2 * M_PI * x_coord / (instance->model.xmax - instance->model.xmin)); +} + +int SetPhase(MdoodzInput *input, Coordinates coordinates) { + const double radius = input->model.user1 / input->scaling.L; + if (coordinates.x * coordinates.x + coordinates.z * coordinates.z < radius * radius) { + return 1; + } else { + return 0; + } +} + +double SetDensity(MdoodzInput *input, Coordinates coordinates, int phase) { + const double T_init = (input->model.user0 + zeroC) / input->scaling.T; + if (1 == 0) { + return input->materials.rho[phase] * (1 - input->materials.alp[phase] * (T_init - input->materials.T0[phase])); + } else { + return input->materials.rho[phase]; + } +} + + +SetBC SetBCVx(MdoodzInput *instance, POSITION position, Coordinates coord) { + SetBC bc; + const double radius = instance->model.user1 / instance->scaling.L; + double x, z; + if (position == W) { + x = coord.x; + z = coord.z; + bc.type = 2; + // bc.value = instance->stress->west*1.1; + bc.value = 1*z*4; + // bc.type = 0; + // bc.value = 0.0; + } else if (position == E) { + x = coord.x; + z = coord.z; + bc.type = 2; + // bc.value = instance->stress->east*1.1;//1*z*4; + bc.value = 1*z*4; + } else if (position == SE || position == SW) { + x = coord.x; + z = coord.z + instance->model.dx / 2.0;// make sure it lives on the boundary + bc.type = 13; + bc.value = 0.0; + } else if (position == S) { + x = coord.x; + z = coord.z + instance->model.dx / 2.0;// make sure it lives on the boundary + bc.type = 13; + bc.value = 0*3e-3; + // Pick a point in the middle + if ( (fabs(x)-0.0) < instance->model.dx/2) { + bc.type = 11; + bc.value = 0*3e-3; + } + } else if (position == N || position == NE || position == NW) { + x = coord.x; + z = coord.z - instance->model.dx / 2.0;// make sure it lives on the boundary + bc.type = 11; + bc.value = 0.0; + } else { + bc.type = -1; + bc.value = 0.0; + } + return bc; +} + +SetBC SetBCVz(MdoodzInput *instance, POSITION position, Coordinates coord) { + SetBC bc; + double x, z; + if (position == N || position == NE || position == NW || position == SE || position == SW) { + x = coord.x; + z = coord.z; + bc.type = 0; + bc.value = 0.0; + } + else if (position == S) { + x = coord.x; + z = coord.z; + bc.type = 0; + bc.value = 0.0; + if (fabs(x)<0.1) { + bc.type = 2; + bc.value = -1.5; + } + } + else if (position == W) { + x = coord.x + instance->model.dx / 2.0; + z = coord.z; + bc.type = 13; + bc.value = 0.0; + } + else if (position == E) { + x = coord.x - instance->model.dx / 2.0; + z = coord.z; + bc.type = 13; + bc.value = 0.0; + } else { + bc.type = -1; + bc.value = 0.0; + } + return bc; +} + +int main(int nargs, char *args[]) { + // Input file name + char *input_file; + if ( nargs < 2 ) { + asprintf(&input_file, DefaultTextFilename(__FILE__)); // Default + } + else { + asprintf(&input_file, "%s", args[1]); // Custom + } + printf("Running MDoodz7.0 using %s\n", input_file); + MdoodzSetup setup = { + .BuildInitialTopography = &(BuildInitialTopography_ff){ + .SetSurfaceZCoord = SetSurfaceZCoord, + }, + .SetParticles = &(SetParticles_ff){ + .SetPhase = SetPhase, + .SetDensity = SetDensity, + }, + .SetBCs = &(SetBCs_ff){ + .SetBCVx = SetBCVx, + .SetBCVz = SetBCVz, + }, + }; + RunMDOODZ(input_file, &setup); + free(input_file); +} diff --git a/SETS/StressBC_FreeSurf.txt b/SETS/StressBC_FreeSurf.txt new file mode 100755 index 00000000..8c013085 --- /dev/null +++ b/SETS/StressBC_FreeSurf.txt @@ -0,0 +1,102 @@ +/**** RESTART ****/ +istep = 00001 +irestart = 0 + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 1 +writer_markers = 1 +writer_debug = 1 +writer_energies = 0 + +/**** SCALES ****/ +eta = 1 +L = 1 +V = 1 +T = 1 + +/**** SPACE ****/ +Nx = 101 +Nz = 101 +Nt = 1 +xmin = -5.000000e-1 +zmin = -5.000000e-1 +xmax = 5.000000e-1 +zmax = 5.000000e-1 + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 + +/**** TIME ****/ +constant_dt = 1 +Courant = 0.3 +dt = 0.5 +RK = 4 + +/**** LINEAR SOLVER ****/ +noisy = 0 +penalty = 1e5 +diag_scaling = 0 +lin_abs_div = 1e-9 +lin_rel_div = 1e-9 +preconditioner = -1 + +/**** NON-LINEAR SOLVER ****/ +Newton = 0 +line_search = 0 +nit_max = 1 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-8 +nonlin_abs_mom = 1e-8 +nonlin_rel_div = 1e-8 +nonlin_rel_div = 1e-8 + +/**** VISCOSITY CUT-OFF ****/ +min_eta = 1e-3 +max_eta = 1e6 + +/**** SWITCHES ****/ +compressible = 1 +elastic = 0 +free_surface = 1 +free_surface_stab = 1 +finite_strain = 1 +advection = 0 +gnuplot_log_res = 0 +eta_average = 0 + +/**** SETUP-DEPENDANT ****/ +shear_style = 0 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1e-10 +user0 = 0 +user1 = 0.15 / inclusion radius [m] +user2 = 0 +user3 = 0 + +/**** GRAVITY ****/ +gx = 0.0000 +gz = -1.000 + +/**** MAT PROPERTIES ****/ +Nb_phases = 2 + +/**** PHASE 1 ****/ +ID = 0 +rho = 1e-0 +plast = 0 +cstv = 1 / constant visc law +eta0 = 1 +npwl = 1.0 +Qpwl = 0 + +/**** PHASE 2 ****/ +ID = 1 +rho = 1e-0 +plast = 0 +cstv = 1 / constant visc law +eta0 = 1 +npwl = 1.0 +Qpwl = 0 diff --git a/SETS/StressBC_ShearTemplate.c b/SETS/StressBC_ShearTemplate.c new file mode 100755 index 00000000..b04be35b --- /dev/null +++ b/SETS/StressBC_ShearTemplate.c @@ -0,0 +1,216 @@ +#include "mdoodz.h" +#include "stdio.h" +#include "stdlib.h" +#include "math.h" +#include "string.h" + +double SetSurfaceZCoord(MdoodzInput *instance, double x_coord) { + const double TopoLevel = 10 / instance->scaling.L; + const double h_pert = instance->model.user3 / instance->scaling.L; + return TopoLevel + 0.0*x_coord;// + h_pert * (3330.0 - 2800.0) / 2800.0 * cos(2 * M_PI * x_coord / (instance->model.xmax - instance->model.xmin)); +} + +int SetPhase(MdoodzInput *input, Coordinates coordinates) { + const double radius = input->model.user1 / input->scaling.L; + if (coordinates.x * coordinates.x + coordinates.z * coordinates.z < radius * radius) { + return 1; + } else { + return 0; + } +} + +double SetDensity(MdoodzInput *input, Coordinates coordinates, int phase) { + const double T_init = (input->model.user0 + zeroC) / input->scaling.T; + if (1 == 0) { + return input->materials.rho[phase] * (1 - input->materials.alp[phase] * (T_init - input->materials.T0[phase])); + } else { + return input->materials.rho[phase]; + } +} + + +SetBC SetBCVx(MdoodzInput *instance, POSITION position, Coordinates coord) { + SetBC bc; + const double radius = instance->model.user1 / instance->scaling.L; + double x, z; + // ----------- Pure shear ----------- // + if (instance->model.shear_style==0) { + if (position == W) { + x = coord.x; + z = coord.z; + bc.type = 2; + bc.value = -1.;//instance->stress->west*0.0;//1*z*4; + } else if (position == E) { + x = coord.x; + z = coord.z; + bc.type = 2; + bc.value = -1.;//instance->stress->east*0.0;//1*z*4; + } else if (position == S || position == SE || position == SW) { + x = coord.x; + z = coord.z + instance->model.dx / 2.0;// make sure it lives on the boundary + bc.type = 13; + bc.value = 0*3e-3; + // Pick a point in the middle + if ( (fabs(x)-0.0) < instance->model.dx/2) { + printf("COUCOU x S\n"); + bc.type = 11; + bc.value = 0*3e-3; + } + } else if (position == N || position == NE || position == NW) { + x = coord.x; + z = coord.z - instance->model.dx / 2.0;// make sure it lives on the boundary + bc.type = 13; + bc.value = 0.0; + // Pick a point in the middle + if ( (fabs(x)-0.0) < instance->model.dx/2) { + printf("COUCOU x N\n"); + bc.type = 11; + bc.value = 0*3e-3; + } + } else { + bc.type = -1; + bc.value = 0.0; + } + } + // ----------- Simple shear ----------- // + else { + const double Lz = (double) (instance->model.zmax - instance->model.zmin); + if (position == S || position == SE || position == SW) { + bc.type = 11; + bc.value = -instance->model.bkg_strain_rate * Lz; + } else if (position == N || position == NE || position == NW) { + bc.type = 11; + bc.value = instance->model.bkg_strain_rate * Lz; + } else if (position == E) { + // bc.type = 0; + // bc.value = 0.0; + bc.type = -12; + bc.value = 0.0; + } else if (position == W) { + // bc.type = 0; + // bc.value = 0.0; + bc.type = -2; + bc.value = 0.0; + } else { + bc.type = -1; + bc.value = 0.0; + } + } + return bc; +} + +SetBC SetBCVz(MdoodzInput *instance, POSITION position, Coordinates coord) { + SetBC bc; + double x, z; + // ----------- Pure shear ----------- // + if (instance->model.shear_style==0) { + if (position == N) { + x = coord.x; + z = coord.z; + bc.type = 2; + bc.value = 1.; + } + else if (position == S) { + x = coord.x; + z = coord.z; + bc.type = 2; + bc.value = 1.; + } + else if (position == W || position == SE || position == SW) { + x = coord.x + instance->model.dx / 2.0; + z = coord.z; + bc.type = 13; + bc.value = 0.0; + // Pick a point in the middle + if ( (fabs(z)-0.0) < instance->model.dz/2) { + printf("COUCOU z W\n"); + bc.type = 11; + bc.value = 0; + } + } + else if (position == E || position == NE || position == NW) { + x = coord.x - instance->model.dx / 2.0; + z = coord.z; + bc.type = 13; + bc.value = 0.0; + // Pick a point in the middle + if ( (fabs(z)-0.0) < instance->model.dz/2) { + printf("COUCOU z E\n"); + bc.type = 11; + bc.value = 0; + } + } else { + bc.type = -1; + bc.value = 0.0; + } + } + // ----------- Simple shear ----------- // + else { + if ( position == W || position == NW || position == SW) { + bc.type = -12; + bc.value = 0.0; + // bc.type = 13; + // bc.value = 0.0; + // // Pick a point in the middle + // if ( (fabs(coord.z)-0.0) < instance->model.dz/2) { + // printf("COUCOU z W\n"); + // bc.type = 11; + // bc.value = 0; + // } + } + else if (position == E || position == NE || position == SE) { + bc.type = -12; + bc.value = 0.0; + // bc.type = 13; + // bc.value = 0.0; + } else if (position == S) { + bc.type = 2; + bc.value = 1.0; + // if ( (fabs(coord.x)-0.0) < instance->model.dx/2) { + // bc.type = 2; + // bc.value = 1.0; + // } + } + else if (position == N) { + bc.type = 0; + bc.value = 0.0; + // if ( (fabs(coord.x)-0.0) < instance->model.dx/2) { + // bc.type = 2; + // bc.value = 1.0; + // } + + } else { + bc.type = -1; + bc.value = 0.0; + } + } + return bc; +} + +int main(int nargs, char *args[]) { + // Input file name + char *input_file; + + if ( nargs < 2 ) { + asprintf(&input_file, DefaultTextFilename(__FILE__)); // Default + } + else { + asprintf(&input_file, "%s", args[1]); // Custom + } + printf("Running MDoodz7.0 using %s\n", input_file); + MdoodzSetup setup = { + .BuildInitialTopography = &(BuildInitialTopography_ff){ + .SetSurfaceZCoord = SetSurfaceZCoord, + }, + .SetParticles = &(SetParticles_ff){ + .SetPhase = SetPhase, + .SetDensity = SetDensity, + }, + .SetBCs = &(SetBCs_ff){ + .SetBCVx = SetBCVx, + .SetBCVz = SetBCVz, + }, + }; + RunMDOODZ(input_file, &setup); + free(input_file); +} diff --git a/SETS/StressBC_ShearTemplate.txt b/SETS/StressBC_ShearTemplate.txt new file mode 100755 index 00000000..0d67623a --- /dev/null +++ b/SETS/StressBC_ShearTemplate.txt @@ -0,0 +1,102 @@ +/**** RESTART ****/ +istep = 00001 +irestart = 0 + +/**** OUTPUT FILES ****/ +writer = 1 +writer_step = 1 +writer_markers = 1 +writer_debug = 1 +writer_energies = 0 + +/**** SCALES ****/ +eta = 1 +L = 1 +V = 1 +T = 1 + +/**** SPACE ****/ +Nx = 31 +Nz = 31 +Nt = 1 +xmin = -5.000000e-1 +zmin = -5.000000e-1 +xmax = 5.000000e-1 +zmax = 5.000000e-1 + +/**** PARTICLES ****/ +Nx_part = 4 +Nz_part = 4 + +/**** TIME ****/ +constant_dt = 1 +Courant = 0.3 +dt = 0.5 +RK = 4 + +/**** LINEAR SOLVER ****/ +noisy = 1 +penalty = 1e5 +diag_scaling = 0 +lin_abs_div = 1e-9 +lin_rel_div = 1e-9 +preconditioner = -1 + +/**** NON-LINEAR SOLVER ****/ +Newton = 0 +line_search = 0 +nit_max = 1 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-10 +nonlin_abs_mom = 1e-10 +nonlin_rel_div = 1e-10 +nonlin_rel_div = 1e-10 + +/**** VISCOSITY CUT-OFF ****/ +min_eta = 1e-3 +max_eta = 1e6 + +/**** SWITCHES ****/ +compressible = 1 +elastic = 0 +free_surface = 0 +free_surface_stab = 0 +finite_strain = 1 +advection = 0 +gnuplot_log_res = 0 +eta_average = 0 + +/**** SETUP-DEPENDANT ****/ +shear_style = 1 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1.0 +user0 = 0 +user1 = 0.15 / inclusion radius [m] +user2 = 0 +user3 = 0 + +/**** GRAVITY ****/ +gx = 0.0000 +gz = 0.000 + +/**** MAT PROPERTIES ****/ +Nb_phases = 2 + +/**** PHASE 1 ****/ +ID = 0 +rho = 1e-0 +plast = 0 +cstv = 1 / constant visc law +eta0 = 1 +npwl = 1.0 +Qpwl = 0 + +/**** PHASE 2 ****/ +ID = 1 +rho = 1e-0 +plast = 0 +cstv = 1 / constant visc law +eta0 = 1e3 +npwl = 1.0 +Qpwl = 0 diff --git a/SETS/ViscousInclusion.c b/SETS/ViscousInclusion.c index 8e7d350b..a5d4a800 100755 --- a/SETS/ViscousInclusion.c +++ b/SETS/ViscousInclusion.c @@ -89,13 +89,15 @@ SetBC SetBCVx(MdoodzInput *instance, POSITION position, Coordinates coord) { const double mc = 1e3; double Vx, Vz, P, eta, sxx, szz, x, z; if (instance->model.shear_style == 0) { - if (position == W || position == E) { - x = coord.x; - z = coord.z; - eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, x, z, 1, radius, mm, mc); + if (position == W) { + eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, coord.x, coord.z, 1, radius, mm, mc); + bc.type = 2; + bc.value = sxx; + } else if (position == E) { + eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, coord.x, coord.z, 1, radius, mm, mc); bc.type = 0; bc.value = Vx; - } else if (position == S || position == SE || position == SW) { + } else if (position == S || position == SE || position == SW) { x = coord.x; z = coord.z + instance->model.dx / 2.0;// make sure it lives on the boundary eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, x, z, 1, radius, mm, mc); @@ -134,29 +136,6 @@ SetBC SetBCVx(MdoodzInput *instance, POSITION position, Coordinates coord) { bc.type = -1; bc.value = 0.0; } - // const double Lz = (double) (instance->model.zmax - instance->model.zmin); - // if (position == S || position == SE || position == SW) { - // x = coord.x; - // z = coord.z + instance->model.dz / 2.0;// make sure it lives on the boundary - // eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, x, z, 0, radius, mm, mc); - // bc.type = 11; - // bc.value = Vx; - // } else if (position == N || position == NE || position == NW) { - // x = coord.x; - // z = coord.z - instance->model.dz / 2.0;// make sure it lives on the boundary - // eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, x, z, 0, radius, mm, mc); - // bc.type = 11; - // bc.value = Vx; - // } else if (position == E) { - // bc.value = 0.0; - // bc.type = -12; - // } else if (position == W) { - // bc.value = 0.0; - // bc.type = -2; - // } else { - // bc.value = 0.0; - // bc.type = -1; - // } } return bc; } @@ -168,21 +147,24 @@ SetBC SetBCVz(MdoodzInput *instance, POSITION position, Coordinates coord) { const double mc = 1e3; double Vx, Vz, P, eta, sxx, szz, x, z; if (instance->model.shear_style == 0) { - if (position == N || position == NE || position == NW || position == S || position == SE || position == SW) { - x = coord.x; - z = coord.z; - eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, x, z, 1, radius, mm, mc); + if ( position == S) { + eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, coord.x, coord.z, 1, radius, mm, mc); bc.type = 0; bc.value = Vz; } - else if (position == W) { + else if (position == N ) { + eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, coord.x, coord.z, 1, radius, mm, mc); + bc.type = 2; + bc.value = szz; + } + else if (position == W || position == SW || position == NW) { x = coord.x + instance->model.dx / 2.0; z = coord.z; eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, x, z, 1, radius, mm, mc); bc.type = 11; bc.value = Vz; } - else if (position == E) { + else if (position == E || position == SE || position == NE) { x = coord.x - instance->model.dx / 2.0; z = coord.z; eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, x, z, 1, radius, mm, mc); diff --git a/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC.ipynb b/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC.ipynb index aa5271b8..b6b3b26a 100644 --- a/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC.ipynb +++ b/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC.ipynb @@ -320,7 +320,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3.9.5 ('base')", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -335,11 +335,6 @@ "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.6" - }, - "vscode": { - "interpreter": { - "hash": "309baa410e37084482bcba9b39a5b9e635e78b91cd7553c4b48a2d89780d0f88" - } } }, "nbformat": 4, diff --git a/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC_v2.ipynb b/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC_v2.ipynb new file mode 100644 index 00000000..ace47115 --- /dev/null +++ b/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC_v2.ipynb @@ -0,0 +1,897 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Txx = D11*Exx;\n", + "Tzz = D22*Ezz;\n", + "Txz = D33*Gxz;\n", + "P = 0;\n" + ] + } + ], + "source": [ + "# %reset\n", + "from sympy import *\n", + "import numpy as np \n", + "init_printing()\n", + "\n", + "comp, out_of_plane = symbols( 'comp, OOP' )\n", + "iso = 1\n", + "\n", + "# Definitions\n", + "eta, dx, dz = symbols( 'eta, dx, dz' )\n", + "D11, D12, D13, D14, D21, D22, D23, D24 ,D31, D32, D33, D34, D41, D42, D43, D44 = symbols('D11, D12, D13, D14, D21, D22, D23, D24 ,D31, D32, D33, D34, D41, D42, D43, D44')\n", + "Exx, Ezz, Gxz, divV = symbols('Exx, Ezz, Gxz, divV')\n", + "stress_labels = ['Txx', 'Tzz', 'Txz', 'P']\n", + "\n", + "if iso==1:\n", + " # Isotropic\n", + " Dv = Matrix( [ [D11, 0, 0, 0], [0, D22, 0, 0], [0, 0, D33, 0], [0, 0, 0, 0] ] )\n", + "else:\n", + " # Anisotropic\n", + " Dv = Matrix( [ [D11, D12, D13, D14], [D21, D22, D23, D24], [D31, D32, D33, D34], [D41, D42, D43, D44] ] )\n", + "\n", + "# Check constitutive\n", + "Ed = Matrix( [ [Exx], [Ezz], [Gxz], [divV] ] )\n", + "T = Dv*Ed\n", + "for i in range(4):\n", + " print(stress_labels[i] + ' = ' + ccode( T[i]) + ';' )" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "uW = (1.0/3.0)*D11W*inW*(OOP*comp*inW - 3)/pow(dx, 2);\n", + "uC = (1.0/3.0)*(3*pow(dx, 2)*(D33N*inN*(2*DirS + RegN) + D33S*inS*(2*DirN + RegS)) - pow(dz, 2)*(D11E*inE*(OOP*comp*inE - 3) + D11W*inW*(OOP*comp*inW - 3)))/(pow(dx, 2)*pow(dz, 2));\n", + "uE = (1.0/3.0)*D11E*inE*(OOP*comp*inE - 3)/pow(dx, 2);\n", + "uS = -D33S*RegS*inS/pow(dz, 2);\n", + "uN = -D33N*RegN*inN/pow(dz, 2);\n", + "vSW = ((1.0/3.0)*D11W*OOP*comp*pow(inW, 2) - D33S*inS)/(dx*dz);\n", + "vSE = (-1.0/3.0*D11E*OOP*comp*pow(inE, 2) + D33S*inS)/(dx*dz);\n", + "vNW = (-1.0/3.0*D11W*OOP*comp*pow(inW, 2) + D33N*inN)/(dx*dz);\n", + "vNE = ((1.0/3.0)*D11E*OOP*comp*pow(inE, 2) - D33N*inN)/(dx*dz);\n", + "uSW = 0;\n", + "uSE = 0;\n", + "uNW = 0;\n", + "uNE = 0;\n", + "vSWW = 0;\n", + "vSEE = 0;\n", + "vNWW = 0;\n", + "vNEE = 0;\n", + "vSSW = 0;\n", + "vSSE = 0;\n", + "vNNW = 0;\n", + "vNNE = 0;\n", + "pW = -1/dx;\n", + "pE = 1.0/dx;\n", + "pSW = 0;\n", + "pSE = 0;\n", + "pNW = 0;\n", + "pNE = 0;\n" + ] + } + ], + "source": [ + "# Stencil x\n", + "uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE = symbols('uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE')\n", + "dofs = Matrix([uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE])\n", + "\n", + "# Stencil extension 1\n", + "uSW,uSE,uNW,uNE,vSWW,vSEE,vNWW,vNEE,vSSW,vSSE,vNNW,vNNE = symbols('uSW,uSE,uNW,uNE,vSWW,vSEE,vNWW,vNEE,vSSW,vSSE,vNNW,vNNE')\n", + "dofs_ext = Matrix([uSW,uSE,uNW,uNE,vSWW,vSEE,vNWW,vNEE,vSSW,vSSE,vNNW,vNNE])\n", + "dofs = np.append(dofs,dofs_ext)\n", + "\n", + "# Stencil extension 2\n", + "pW, pE, pSW, pSE, pNW, pNE = symbols('pW, pE, pSW, pSE, pNW, pNE')\n", + "dofs_ext = Matrix([pW, pE, pSW, pSE, pNW, pNE])\n", + "dofs = np.append(dofs,dofs_ext)\n", + "\n", + "# Flags\n", + "inS,inN,inW,inE = symbols('inS,inN,inW,inE')\n", + "inSEc,inSWc,inNEc,inNWc = symbols('inSEc,inSWc,inNEc,inNWc')\n", + "inSEv,inSWv,inNEv,inNWv = symbols('inSEv,inSWv,inNEv,inNWv')\n", + "NeuN, NeuS, DirN, DirS, RegN, RegS = symbols('NeuN, NeuS, DirN, DirS, RegN, RegS')\n", + "DirSW, DirSE, NeuSW, NeuSE = symbols('DirSW, DirSE, NeuSW, NeuSE')\n", + "DirNW, DirNE, NeuNW, NeuNE = symbols('DirNW, DirNE, NeuNW, NeuNE')\n", + "\n", + "# North/South dVxdz\n", + "dVxdzN_Reg = (uN - uC )/dz\n", + "dVxdzS_Reg = (uC - uS )/dz\n", + "dVxdzN_Dir = ((-uC) - uC )/dz\n", + "dVxdzS_Dir = (uC - (-uC) )/dz\n", + "dVxdzN_Neu = ((uC) - uC )/dz\n", + "dVxdzS_Neu = (uC - (uC) )/dz\n", + "\n", + "dVzdxN_Reg = (vNE -vNW )/dx\n", + "dVzdxS_Reg = (vSE -vSW )/dx\n", + "\n", + "# Divergences\n", + "divW = inW *out_of_plane*( (uC - uW )/dx + (vNW - vSW )/dz ) \n", + "divE = inE *out_of_plane*( (uE - uC )/dx + (vNE - vSE )/dz )\n", + "divSW = inSWc*out_of_plane*( (uS - uSW)/dx + (vSW - vSSW)/dz )\n", + "divSE = inSEc*out_of_plane*( (uSE - uS )/dx + (vSE - vSSE)/dz )\n", + "divNW = inNWc*out_of_plane*( (uN - uNW)/dx + (vNNW- vNW )/dz )\n", + "divNE = inNEc*out_of_plane*( (uNE - uN )/dx + (vNNE- vNE )/dz )\n", + "\n", + "# Deviatoric normal strain rate\n", + "ExxW = inW *((uC -uW )/dx - comp*Rational(1,3)*divW) \n", + "EzzW = inW *((vNW-vSW)/dz - comp*Rational(1,3)*divW)\n", + "ExxE = inE *((uE -uC )/dx - comp*Rational(1,3)*divE)\n", + "EzzE = inE *((vNE-vSE)/dz - comp*Rational(1,3)*divE)\n", + "\n", + "# Additional missing stress tensor components needed for interpolation\n", + "GxzNE = inNEv*((uNE - uE )/dz + ( vNEE - vNE )/dx)\n", + "GxzSE = inSEv*((uE - uSE)/dz + ( vSEE - vSE )/dx)\n", + "GxzNW = inNWv*((uNW - uW )/dz + ( vNW - vNWW )/dx)\n", + "GxzSW = inSWv*((uW - uSW)/dz + ( vSW - vSWW )/dx)\n", + "\n", + "\n", + "# Shear strain rate\n", + "dVxdzN = RegN*dVxdzN_Reg + DirS*dVxdzN_Dir + NeuS*dVxdzN_Neu\n", + "dVxdzS = RegS*dVxdzS_Reg + DirN*dVxdzS_Dir + NeuN*dVxdzS_Neu\n", + "dVzdxN = dVzdxN_Reg\n", + "dVzdxS = dVzdxS_Reg\n", + "\n", + "GxzN = inN*(dVxdzN + dVzdxN) \n", + "GxzS = inS*(dVxdzS + dVzdxS) \n", + "ExxNE = inNEc*((uNE -uN )/dx - comp*Rational(1,3)*divNE)\n", + "ExxNW = inNWc*((uN -uNW)/dx - comp*Rational(1,3)*divNW)\n", + "ExxSE = inSEc*((uSE -uS )/dx - comp*Rational(1,3)*divSE)\n", + "ExxSW = inSWc*((uS -uSW)/dx - comp*Rational(1,3)*divSW)\n", + "EzzNE = inNEc*((vNNE-vNE)/dz - comp*Rational(1,3)*divNE)\n", + "EzzNW = inNWc*((vNNW-vNW)/dz - comp*Rational(1,3)*divNW)\n", + "EzzSE = inSEc*((vSE-vSSE)/dz - comp*Rational(1,3)*divSE)\n", + "EzzSW = inSWc*((vSW-vSSW)/dz - comp*Rational(1,3)*divSW)\n", + "# INV 0\n", + "GxzE = 0.25*( inN*GxzN + inS*GxzS + inNEv*GxzNE + inSEv*GxzSE) # clear\n", + "GxzW = 0.25*( inN*GxzN + inS*GxzS + inNWv*GxzNW + inSWv*GxzSW)\n", + "ExxN = 0.25*( inE*ExxE + inW*ExxW + inNWc*ExxNW + inNEc*ExxNE)\n", + "ExxS = 0.25*( inE*ExxE + inW*ExxW + inSWc*ExxSW + inSEc*ExxSE)\n", + "EzzN = 0.25*( inE*EzzE + inW*EzzW + inNWc*EzzNW + inNEc*EzzNE) \n", + "EzzS = 0.25*( inE*EzzE + inW*EzzW + inSWc*EzzSW + inSEc*EzzSE)\n", + "pS = 0.25*( inE*pE + inW*pW + inSWc*pSW + inSEc*pSE)\n", + "pN = 0.25*( inE*pE + inW*pW + inSWc*pNW + inSEc*pNE) \n", + "\n", + "#----------------------------------------#\n", + "EdW = Matrix([[ExxW], [EzzW], [GxzW], [pW]])\n", + "TW = Dv*EdW\n", + "TxxW = TW[0].subs({D11:'D11W', D12:'D12W', D13:'D13W', D14:'D14W'})\n", + "#----------------------------------------#\n", + "EdE = Matrix([[ExxE], [EzzE], [GxzE], [pE]])\n", + "TE = Dv*EdE\n", + "TxxE = TE[0].subs({D11:'D11E', D12:'D12E', D13:'D13E', D14:'D14E'})\n", + "#----------------------------------------#\n", + "EdS = Matrix([[ExxS], [EzzS], [GxzS], [pS]])\n", + "TS = Dv*EdS\n", + "TxzS = TS[2].subs({D31:'D31S', D32:'D32S', D33:'D33S', D34:'D34S'})\n", + "#----------------------------------------#\n", + "EdN = Matrix([[ExxN], [EzzN], [GxzN], [pN]])\n", + "TN = Dv*EdN\n", + "TxzN = TN[2].subs({D31:'D31N', D32:'D32N', D33:'D33N', D34:'D34N'})\n", + "#----------------------------------------#\n", + "Fx = 1/dx*(TxxE - TxxW) + 1/dz*(TxzN - TxzS) - 1/dx*(pE - pW)\n", + "Fx = -Fx\n", + "\n", + "for i in range(len(dofs)):\n", + " cUc = Fx.diff(dofs[i])\n", + " print(str(dofs[i]) + ' = ' + ccode(simplify(cUc)) + ';')" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "uW = 0;\n", + "uC = (1.0/6.0)*(-2*D11E*pow(dz, 2)*inE*(OOP*comp*inE - 3) + 3*pow(dx, 2)*(D33N*inN*(2*DirS + RegN) + D33S*inS*(2*DirN + RegS)))/(pow(dx, 2)*pow(dz, 2));\n", + "uE = (1.0/3.0)*D11E*inE*(OOP*comp*inE - 3)/pow(dx, 2);\n", + "uS = -1.0/2.0*D33S*RegS*inS/pow(dz, 2);\n", + "uN = -1.0/2.0*D33N*RegN*inN/pow(dz, 2);\n", + "vSW = 0;\n", + "vSE = (-1.0/3.0*D11E*OOP*comp*pow(inE, 2) + D33S*DirSW*inS)/(dx*dz);\n", + "vNW = 0;\n", + "vNE = ((1.0/3.0)*D11E*OOP*comp*pow(inE, 2) - D33N*DirNW*inN)/(dx*dz);\n", + "uSW = 0;\n", + "uSE = 0;\n", + "uNW = 0;\n", + "uNE = 0;\n", + "vSWW = 0;\n", + "vSEE = 0;\n", + "vNWW = 0;\n", + "vNEE = 0;\n", + "vSSW = 0;\n", + "vSSE = 0;\n", + "vNNW = 0;\n", + "vNNE = 0;\n", + "pW = 0;\n", + "pE = 1.0/dx;\n", + "pSW = 0;\n", + "pSE = 0;\n", + "pNW = 0;\n", + "pNE = 0;\n" + ] + } + ], + "source": [ + "# Compressibility switch\n", + "SxxBC = symbols('SxxBC') # boundary stress\n", + "\n", + "# Stencil x\n", + "uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE = symbols('uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE')\n", + "dofs = Matrix([uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE])\n", + "\n", + "# Stencil extension 1\n", + "uSW,uSE,uNW,uNE,vSWW,vSEE,vNWW,vNEE,vSSW,vSSE,vNNW,vNNE = symbols('uSW,uSE,uNW,uNE,vSWW,vSEE,vNWW,vNEE,vSSW,vSSE,vNNW,vNNE')\n", + "dofs_ext = Matrix([uSW,uSE,uNW,uNE,vSWW,vSEE,vNWW,vNEE,vSSW,vSSE,vNNW,vNNE])\n", + "dofs = np.append(dofs,dofs_ext)\n", + "\n", + "# Stencil extension 2\n", + "pW, pE, pSW, pSE, pNW, pNE = symbols('pW, pE, pSW, pSE, pNW, pNE')\n", + "dofs_ext = Matrix([pW, pE, pSW, pSE, pNW, pNE])\n", + "dofs = np.append(dofs,dofs_ext)\n", + "\n", + "# Flags\n", + "inS,inN,inW,inE = symbols('inS,inN,inW,inE')\n", + "inSEc,inSWc,inNEc,inNWc = symbols('inSEc,inSWc,inNEc,inNWc')\n", + "inSEv,inSWv,inNEv,inNWv = symbols('inSEv,inSWv,inNEv,inNWv')\n", + "NeuN, NeuS, DirN, DirS, RegN, RegS = symbols('NeuN, NeuS, DirN, DirS, RegN, RegS')\n", + "DirSW, DirSE, NeuSW, NeuSE = symbols('DirSW, DirSE, NeuSW, NeuSE')\n", + "DirNW, DirNE, NeuNW, NeuNE = symbols('DirNW, DirNE, NeuNW, NeuNE')\n", + "\n", + "# North/South dVxdz\n", + "dVxdzN_Reg = (uN - uC )/dz\n", + "dVxdzS_Reg = (uC - uS )/dz\n", + "dVxdzN_Dir = ((-uC) - uC )/dz\n", + "dVxdzS_Dir = (uC - (-uC) )/dz\n", + "dVxdzN_Neu = ((uC) - uC )/dz\n", + "dVxdzS_Neu = (uC - (uC) )/dz\n", + "\n", + "dVzdxN_Reg = (vNE -vNW )/dx\n", + "dVzdxN_DirW = (vNE -(-vNE) )/dx\n", + "dVzdxN_NeuW = (vNE -vNE )/dx\n", + "dVzdxN_DirE = (-vNW -vNW )/dx\n", + "dVzdxN_NeuE = (vNW -vNW )/dx\n", + "dVzdxS_Reg = (vSE -vSW )/dx\n", + "dVzdxS_DirW = (vSE -(-vSE) )/dx\n", + "dVzdxS_NeuW = (vSE -vSE )/dx\n", + "dVzdxS_DirE = ((-vSW) -vSW )/dx\n", + "dVzdxS_NeuE = (vSW -vSW )/dx\n", + "\n", + "# Divergences\n", + "divW = inW *out_of_plane*( (uC - uW )/dx + (vNW - vSW )/dz ) \n", + "divE = inE *out_of_plane*( (uE - uC )/dx + (vNE - vSE )/dz )\n", + "divSW = inSWc*out_of_plane*( (uS - uSW)/dx + (vSW - vSSW)/dz )\n", + "divSE = inSEc*out_of_plane*( (uSE - uS )/dx + (vSE - vSSE)/dz )\n", + "divNW = inNWc*out_of_plane*( (uN - uNW)/dx + (vNNW- vNW )/dz )\n", + "divNE = inNEc*out_of_plane*( (uNE - uN )/dx + (vNNE- vNE )/dz )\n", + "\n", + "# Deviatoric normal strain rate\n", + "ExxW = inW *((uC -uW )/dx - comp*Rational(1,3)*divW) \n", + "EzzW = inW *((vNW-vSW)/dz - comp*Rational(1,3)*divW)\n", + "ExxE = inE *((uE -uC )/dx - comp*Rational(1,3)*divE)\n", + "EzzE = inE *((vNE-vSE)/dz - comp*Rational(1,3)*divE)\n", + "\n", + "# Additional missing stress tensor components needed for interpolation\n", + "GxzNE = inNEv*((uNE - uE )/dz + ( vNEE - vNE )/dx)\n", + "GxzSE = inSEv*((uE - uSE)/dz + ( vSEE - vSE )/dx)\n", + "GxzNW = inNWv*((uNW - uW )/dz + ( vNW - vNWW )/dx)\n", + "GxzSW = inSWv*((uW - uSW)/dz + ( vSW - vSWW )/dx)\n", + "\n", + "\n", + "# Shear strain rate\n", + "dVxdzN = RegN*dVxdzN_Reg + DirS*dVxdzN_Dir + NeuS*dVxdzN_Neu\n", + "dVxdzS = RegS*dVxdzS_Reg + DirN*dVxdzS_Dir + NeuN*dVxdzS_Neu\n", + "\n", + "dVzdxN = DirNW*dVzdxN_DirW + NeuNW*dVzdxN_NeuW\n", + "dVzdxS = DirSW*dVzdxS_DirW + NeuSW*dVzdxS_NeuW\n", + "# dVzdxN = dVzdxN_Reg\n", + "# dVzdxS = dVzdxS_Reg\n", + "\n", + "GxzN = inN*(dVxdzN + dVzdxN) \n", + "GxzS = inS*(dVxdzS + dVzdxS) \n", + "ExxNE = inNEc*((uNE -uN )/dx - comp*Rational(1,3)*divNE)\n", + "ExxNW = inNWc*((uN -uNW)/dx - comp*Rational(1,3)*divNW)\n", + "ExxSE = inSEc*((uSE -uS )/dx - comp*Rational(1,3)*divSE)\n", + "ExxSW = inSWc*((uS -uSW)/dx - comp*Rational(1,3)*divSW)\n", + "EzzNE = inNEc*((vNNE-vNE)/dz - comp*Rational(1,3)*divNE)\n", + "EzzNW = inNWc*((vNNW-vNW)/dz - comp*Rational(1,3)*divNW)\n", + "EzzSE = inSEc*((vSE-vSSE)/dz - comp*Rational(1,3)*divSE)\n", + "EzzSW = inSWc*((vSW-vSSW)/dz - comp*Rational(1,3)*divSW)\n", + "# INV 0\n", + "GxzE = 0.25*( inN*GxzN + inS*GxzS + inNEv*GxzNE + inSEv*GxzSE) # clear\n", + "GxzW = 0.25*( inN*GxzN + inS*GxzS + inNWv*GxzNW + inSWv*GxzSW)\n", + "ExxN = 0.25*( inE*ExxE + inW*ExxW + inNWc*ExxNW + inNEc*ExxNE)\n", + "ExxS = 0.25*( inE*ExxE + inW*ExxW + inSWc*ExxSW + inSEc*ExxSE)\n", + "EzzN = 0.25*( inE*EzzE + inW*EzzW + inNWc*EzzNW + inNEc*EzzNE) \n", + "EzzS = 0.25*( inE*EzzE + inW*EzzW + inSWc*EzzSW + inSEc*EzzSE)\n", + "pS = 0.25*( inE*pE + inW*pW + inSWc*pSW + inSEc*pSE)\n", + "pN = 0.25*( inE*pE + inW*pW + inSWc*pNW + inSEc*pNE) \n", + "\n", + "#----------------------------------------#\n", + "EdW = Matrix([[ExxW], [EzzW], [GxzW], [pW]])\n", + "TW = Dv*EdW\n", + "TxxW = TW[0].subs({D11:'D11W', D12:'D12W', D13:'D13W', D14:'D14W'})\n", + "#----------------------------------------#\n", + "EdE = Matrix([[ExxE], [EzzE], [GxzE], [pE]])\n", + "TE = Dv*EdE\n", + "TxxE = TE[0].subs({D11:'D11E', D12:'D12E', D13:'D13E', D14:'D14E'})\n", + "#----------------------------------------#\n", + "EdS = Matrix([[ExxS], [EzzS], [GxzS], [pS]])\n", + "TS = Dv*EdS\n", + "TxzS = TS[2].subs({D31:'D31S', D32:'D32S', D33:'D33S', D34:'D34S'})\n", + "#----------------------------------------#\n", + "EdN = Matrix([[ExxN], [EzzN], [GxzN], [pN]])\n", + "TN = Dv*EdN\n", + "TxzN = TN[2].subs({D31:'D31N', D32:'D32N', D33:'D33N', D34:'D34N'})\n", + "#----------------------------------------#\n", + "Fx = 1/dx*(TxxE - TxxW) + 1/dz*(TxzN - TxzS) - 1/dx*(pE - pW)\n", + "Fx = -Fx\n", + "\n", + "# Stress boundary\n", + "SxxE = TxxE - pE\n", + "SxxW = 2*SxxBC - SxxE\n", + "Fx = 1/dx*(SxxE - SxxW) + 1/dz*(TxzN - TxzS)\n", + "Fx = -Fx/2\n", + "\n", + "for i in range(len(dofs)):\n", + " cUc = Fx.diff(dofs[i])\n", + " print(str(dofs[i]) + ' = ' + ccode(simplify(cUc)) + ';')\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "vW = -D33W*RegW*inW/pow(dx, 2);\n", + "vC = (1.0/3.0)*(-pow(dx, 2)*(D22N + D22S)*(OOP*comp - 3) + 3*pow(dz, 2)*(D33E*inE*(2*DirE + RegE) + D33W*inW*(2*DirW + RegW)))/(pow(dx, 2)*pow(dz, 2));\n", + "vE = -D33E*RegE*inE/pow(dx, 2);\n", + "vS = (1.0/3.0)*D22S*(OOP*comp - 3)/pow(dz, 2);\n", + "vN = (1.0/3.0)*D22N*(OOP*comp - 3)/pow(dz, 2);\n", + "uSW = ((1.0/3.0)*D22S*OOP*comp - D33W*inW)/(dx*dz);\n", + "uSE = (-1.0/3.0*D22S*OOP*comp + D33E*inE)/(dx*dz);\n", + "uNW = (-1.0/3.0*D22N*OOP*comp + D33W*inW)/(dx*dz);\n", + "uNE = ((1.0/3.0)*D22N*OOP*comp - D33E*inE)/(dx*dz);\n", + "vSW = 0;\n", + "vSE = 0;\n", + "vNW = 0;\n", + "vNE = 0;\n", + "uSWW = 0;\n", + "uSEE = 0;\n", + "uNWW = 0;\n", + "uNEE = 0;\n", + "uSSW = 0;\n", + "uSSE = 0;\n", + "uNNW = 0;\n", + "uNNE = 0;\n", + "pS = -1/dz;\n", + "pN = 1.0/dz;\n", + "pSW = 0;\n", + "pSE = 0;\n", + "pNW = 0;\n", + "pNE = 0;\n" + ] + } + ], + "source": [ + "# Compressibility switch\n", + "comp = symbols('comp')\n", + "\n", + "# Stencil z\n", + "vW,vC,vE,vS,vN,uSW,uSE,uNW,uNE = symbols('vW,vC,vE,vS,vN,uSW,uSE,uNW,uNE')\n", + "dofs = Matrix([vW,vC,vE,vS,vN,uSW,uSE,uNW,uNE])\n", + "\n", + "# Stencil extension 1\n", + "vSW,vSE,vNW,vNE,uSWW,uSEE,uNWW,uNEE,uSSW,uSSE,uNNW,uNNE = symbols('vSW,vSE,vNW,vNE,uSWW,uSEE,uNWW,uNEE,uSSW,uSSE,uNNW,uNNE')\n", + "dofs_ext = Matrix([vSW,vSE,vNW,vNE,uSWW,uSEE,uNWW,uNEE,uSSW,uSSE,uNNW,uNNE])\n", + "dofs = np.append(dofs,dofs_ext)\n", + "\n", + "# Stencil extension 2\n", + "pS, pN, pSW, pSE, pNW, pNE = symbols('pS, pN, pSW, pSE, pNW, pNE')\n", + "dofs_ext = Matrix([pS, pN, pSW, pSE, pNW, pNE])\n", + "dofs = np.append(dofs,dofs_ext)\n", + "\n", + "# Flags\n", + "inS,inN,inW,inE = symbols('inS,inN,inW,inE')\n", + "inSEc,inSWc,inNEc,inNWc = symbols('inSEc,inSWc,inNEc,inNWc')\n", + "inSEv,inSWv,inNEv,inNWv = symbols('inSEv,inSWv,inNEv,inNWv')\n", + "DirE, DirW, NeuE, NeuW, RegE, RegW = symbols('DirE, DirW, NeuE, NeuW, RegE, RegW')\n", + "\n", + "# East/West velocity gradients\n", + "dVzdxE_Reg = (vE - vC )/dx\n", + "dVzdxW_Reg = (vC - vW )/dx\n", + "dVzdxE_Dir = (-vC - vC )/dx\n", + "dVzdxW_Dir = (vC - -vC )/dx\n", + "dVzdxE_Neu = (vC - vC )/dx\n", + "dVzdxW_Neu = (vC - vC )/dx\n", + "\n", + "# Divergences\n", + "divS = out_of_plane*( (uSE-uSW)/dx + (vC -vS )/dz ) \n", + "divN = out_of_plane*( (uNE-uNW)/dx + (vN -vC )/dz ) \n", + "divSW = out_of_plane*( (uSW-uSWW)/dx + (vW -vSW )/dz )\n", + "divSE = out_of_plane*( (uSEE-uSE)/dx + (vE -vSE )/dz )\n", + "divNW = out_of_plane*( (uNW-uNWW)/dx + (vNW -vW )/dz )\n", + "divNE = out_of_plane*( (uNEE-uNE)/dx + (vNE -vE )/dz ) \n", + "\n", + "# Deviatoric normal strain\n", + "ExxS = ((uSE-uSW)/dx - comp*Rational(1,3)*divS)\n", + "EzzS = ((vC -vS )/dz - comp*Rational(1,3)*divS)\n", + "ExxN = ((uNE-uNW)/dx - comp*Rational(1,3)*divN)\n", + "EzzN = ((vN -vC )/dz - comp*Rational(1,3)*divN)\n", + "\n", + "# Additional missing stress tensor components via interpolation\n", + "GxzSW = ((uSW - uSSW)/dz + (vS - vSW)/dx)\n", + "GxzSE = ((uSE - uSSE)/dz + (vSE - vS )/dx)\n", + "GxzNW = ((uNNW - uNW )/dz + (vN - vNW)/dx)\n", + "GxzNE = ((uNNE - uNE )/dz + (vNE - vN )/dx)\n", + "\n", + "# Shear strain rate\n", + "GxzE = inE*((uNE-uSE)/dz + RegE*dVzdxE_Reg + DirE*dVzdxE_Dir + NeuE*dVzdxE_Neu)\n", + "GxzW = inW*((uNW-uSW)/dz + RegW*dVzdxW_Reg + DirW*dVzdxW_Dir + NeuW*dVzdxW_Neu)\n", + "\n", + "# Additional missing stress tensor components for interpolation\n", + "ExxNE = ((uNEE-uNE)/dx - comp*1/3*(divNE))\n", + "ExxNW = ((uNW-uNWW)/dx - comp*1/3*(divNW))\n", + "ExxSE = ((uSEE-uSE)/dx - comp*1/3*(divSE))\n", + "ExxSW = ((uSW-uSWW)/dx - comp*1/3*(divSW))\n", + "EzzNE = ((vNE -vE )/dz - comp*1/3*(divNE))\n", + "EzzNW = ((vNW -vW )/dz - comp*1/3*(divNW))\n", + "EzzSE = ((vE -vSE )/dz - comp*1/3*(divSE))\n", + "EzzSW = ((vW -vSW )/dz - comp*1/3*(divSW))\n", + "\n", + "# INV 0\n", + "GxzN = Rational(1,4)*( inW*GxzW + inE*GxzE + inNWv*GxzNW + inNEv*GxzNE )\n", + "GxzS = Rational(1,4)*( inW*GxzW + inE*GxzE + inSWv*GxzSW + inSEv*GxzSE ) \n", + "ExxE = Rational(1,4)*( inS*ExxS + inN*ExxN + inNEc*ExxNE + inSEc*ExxSE )\n", + "ExxW = Rational(1,4)*( inS*ExxS + inN*ExxN + inNWc*ExxNW + inSWc*ExxSW )\n", + "EzzE = Rational(1,4)*( inS*EzzS + inN*EzzN + inNEc*EzzNE + inSEc*EzzSE )\n", + "EzzW = Rational(1,4)*( inS*EzzS + inN*EzzN + inNWc*EzzNW + inSWc*EzzSW )\n", + "pE = Rational(1,4)*( inS*pS + inN*pN + inNEc*pNE + inSEc*pSE )\n", + "pW = Rational(1,4)*( inS*pS + inN*pN + inNWc*pNW + inSWc*pSW )\n", + "\n", + "#----------------------------------------#\n", + "EdS = Matrix([[ExxS], [EzzS], [GxzS], [pS]])\n", + "TS = Dv*EdS\n", + "TzzS = TS[1].subs({D21:'D21S', D22:'D22S', D23:'D23S', D24:'D24S'})\n", + "#----------------------------------------#\n", + "EdN = Matrix([[ExxN], [EzzN], [GxzN], [pN]])\n", + "TN = Dv*EdN\n", + "TzzN = TN[1].subs({D21:'D21N', D22:'D22N', D23:'D23N', D24:'D24N'})\n", + "#----------------------------------------#\n", + "EdW = Matrix([[ExxW], [EzzW], [GxzW], [pW]])\n", + "TW = Dv*EdW\n", + "TxzW = TW[2].subs({D31:'D31W', D32:'D32W', D33:'D33W', D34:'D34W'})\n", + "#----------------------------------------#\n", + "EdE = Matrix([[ExxE], [EzzE], [GxzE], [pE]])\n", + "TE = Dv*EdE\n", + "TxzE = TE[2].subs({D31:'D31E', D32:'D32E', D33:'D33E', D34:'D34E'})\n", + "#----------------------------------------#\n", + "Fz = 1/dz*(TzzN - TzzS) + 1/dx*(TxzE - TxzW) - (pN - pS)/dz \n", + "Fz = -Fz\n", + "for i in range(len(dofs)):\n", + " cVc = Fz.diff(dofs[i])\n", + " print(str(dofs[i]) + ' = ' + ccode(simplify(cVc)) + ';') " + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "uW = (1.0/3.0)*D11W*inW*(1 - SxxBCW)*(OOP*comp*inW - 3)/pow(dx, 2);\n", + "uC = -1.0/6.0*(SxxBCE*(2*D11W*pow(dz, 2)*inW*(OOP*comp*inW - 3) - 3*pow(dx, 2)*(D33N*inN*(2*DirN + RegN) + D33S*inS*(2*DirS + RegS))) + SxxBCW*(2*D11E*pow(dz, 2)*inE*(OOP*comp*inE - 3) - 3*pow(dx, 2)*(D33N*inN*(2*DirN + RegN) + D33S*inS*(2*DirS + RegS))) + 2*(3*pow(dx, 2)*(D33N*inN*(2*DirN + RegN) + D33S*inS*(2*DirS + RegS)) - pow(dz, 2)*(D11E*inE*(OOP*comp*inE - 3) + D11W*inW*(OOP*comp*inW - 3)))*(SxxBCE + SxxBCW - 1))/(pow(dx, 2)*pow(dz, 2));\n", + "uE = (1.0/3.0)*D11E*inE*(1 - SxxBCE)*(OOP*comp*inE - 3)/pow(dx, 2);\n", + "uS = (1.0/2.0)*D33S*RegS*inS*(SxxBCE + SxxBCW - 2)/pow(dz, 2);\n", + "uN = (1.0/2.0)*D33N*RegN*inN*(SxxBCE + SxxBCW - 2)/pow(dz, 2);\n", + "vSW = (1.0/6.0)*(-3*D33S*SxxBCW*inS*(2*DirSE*SxxBCE - SxxBCE - SxxBCW + 1) + SxxBCE*(2*D11W*OOP*comp*pow(inW, 2) - 3*D33S*inS*(2*DirSE*SxxBCE - SxxBCE - SxxBCW + 1)) - 2*(D11W*OOP*comp*pow(inW, 2) - 3*D33S*inS*(2*DirSE*SxxBCE - SxxBCE - SxxBCW + 1))*(SxxBCE + SxxBCW - 1))/(dx*dz);\n", + "vSE = (1.0/6.0)*(3*D33S*SxxBCE*inS*(2*DirSW*SxxBCW - SxxBCE - SxxBCW + 1) - SxxBCW*(2*D11E*OOP*comp*pow(inE, 2) - 3*D33S*inS*(2*DirSW*SxxBCW - SxxBCE - SxxBCW + 1)) + 2*(D11E*OOP*comp*pow(inE, 2) - 3*D33S*inS*(2*DirSW*SxxBCW - SxxBCE - SxxBCW + 1))*(SxxBCE + SxxBCW - 1))/(dx*dz);\n", + "vNW = (1.0/6.0)*(3*D33N*SxxBCW*inN*(2*DirNE*SxxBCE - SxxBCE - SxxBCW + 1) - SxxBCE*(2*D11W*OOP*comp*pow(inW, 2) - 3*D33N*inN*(2*DirNE*SxxBCE - SxxBCE - SxxBCW + 1)) + 2*(D11W*OOP*comp*pow(inW, 2) - 3*D33N*inN*(2*DirNE*SxxBCE - SxxBCE - SxxBCW + 1))*(SxxBCE + SxxBCW - 1))/(dx*dz);\n", + "vNE = (1.0/6.0)*(-3*D33N*SxxBCE*inN*(2*DirNW*SxxBCW - SxxBCE - SxxBCW + 1) + SxxBCW*(2*D11E*OOP*comp*pow(inE, 2) - 3*D33N*inN*(2*DirNW*SxxBCW - SxxBCE - SxxBCW + 1)) - 2*(D11E*OOP*comp*pow(inE, 2) - 3*D33N*inN*(2*DirNW*SxxBCW - SxxBCE - SxxBCW + 1))*(SxxBCE + SxxBCW - 1))/(dx*dz);\n", + "uSW = 0;\n", + "uSE = 0;\n", + "uNW = 0;\n", + "uNE = 0;\n", + "vSWW = 0;\n", + "vSEE = 0;\n", + "vNWW = 0;\n", + "vNEE = 0;\n", + "vSSW = 0;\n", + "vSSE = 0;\n", + "vNNW = 0;\n", + "vNNE = 0;\n", + "pW = (SxxBCW - 1)/dx;\n", + "pE = (1 - SxxBCE)/dx;\n", + "pSW = 0;\n", + "pSE = 0;\n", + "pNW = 0;\n", + "pNE = 0;\n" + ] + } + ], + "source": [ + "SxxBC = symbols('SxxBC') # boundary stress\n", + "SxxBCW = symbols('SxxBCW') # Flag to switch between BC and internal nodes\n", + "SxxBCE = symbols('SxxBCE') # Flag to switch between BC and internal nodes\n", + "\n", + "# Stencil x\n", + "uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE = symbols('uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE')\n", + "dofs = Matrix([uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE])\n", + "\n", + "# Stencil extension 1\n", + "uSW,uSE,uNW,uNE,vSWW,vSEE,vNWW,vNEE,vSSW,vSSE,vNNW,vNNE = symbols('uSW,uSE,uNW,uNE,vSWW,vSEE,vNWW,vNEE,vSSW,vSSE,vNNW,vNNE')\n", + "dofs_ext = Matrix([uSW,uSE,uNW,uNE,vSWW,vSEE,vNWW,vNEE,vSSW,vSSE,vNNW,vNNE])\n", + "dofs = np.append(dofs,dofs_ext)\n", + "\n", + "# Stencil extension 2\n", + "pW, pE, pSW, pSE, pNW, pNE = symbols('pW, pE, pSW, pSE, pNW, pNE')\n", + "dofs_ext = Matrix([pW, pE, pSW, pSE, pNW, pNE])\n", + "dofs = np.append(dofs,dofs_ext)\n", + "\n", + "# Flags\n", + "inS,inN,inW,inE = symbols('inS,inN,inW,inE')\n", + "inSEc,inSWc,inNEc,inNWc = symbols('inSEc,inSWc,inNEc,inNWc')\n", + "inSEv,inSWv,inNEv,inNWv = symbols('inSEv,inSWv,inNEv,inNWv')\n", + "NeuN, NeuS, DirN, DirS, RegN, RegS = symbols('NeuN, NeuS, DirN, DirS, RegN, RegS')\n", + "DirSW, DirSE, NeuSW, NeuSE = symbols('DirSW, DirSE, NeuSW, NeuSE')\n", + "DirNW, DirNE, NeuNW, NeuNE = symbols('DirNW, DirNE, NeuNW, NeuNE')\n", + "\n", + "# North/South dVxdz\n", + "dVxdzN_Reg = (uN - uC )/dz\n", + "dVxdzS_Reg = (uC - uS )/dz\n", + "dVxdzN_Dir = ((-uC) - uC )/dz\n", + "dVxdzS_Dir = (uC - (-uC) )/dz\n", + "dVxdzN_Neu = ((uC) - uC )/dz\n", + "dVxdzS_Neu = (uC - (uC) )/dz\n", + "\n", + "# North/South dVzdx\n", + "dVzdxN_Reg = (vNE -vNW )/dx\n", + "dVzdxN_DirW = (vNE -(-vNE) )/dx\n", + "dVzdxN_NeuW = (vNE -vNE )/dx\n", + "dVzdxN_DirE = (-vNW -vNW )/dx\n", + "dVzdxN_NeuE = (vNW -vNW )/dx\n", + "dVzdxS_Reg = (vSE -vSW )/dx\n", + "dVzdxS_DirW = (vSE -(-vSE) )/dx\n", + "dVzdxS_NeuW = (vSE -vSE )/dx\n", + "dVzdxS_DirE = ((-vSW) -vSW )/dx\n", + "dVzdxS_NeuE = (vSW -vSW )/dx\n", + "\n", + "# Divergences\n", + "divW = inW *out_of_plane*( (uC - uW )/dx + (vNW - vSW )/dz ) \n", + "divE = inE *out_of_plane*( (uE - uC )/dx + (vNE - vSE )/dz )\n", + "divSW = inSWc*out_of_plane*( (uS - uSW)/dx + (vSW - vSSW)/dz )\n", + "divSE = inSEc*out_of_plane*( (uSE - uS )/dx + (vSE - vSSE)/dz )\n", + "divNW = inNWc*out_of_plane*( (uN - uNW)/dx + (vNNW- vNW )/dz )\n", + "divNE = inNEc*out_of_plane*( (uNE - uN )/dx + (vNNE- vNE )/dz )\n", + "\n", + "# Deviatoric normal strain rate\n", + "ExxW = inW *((uC -uW )/dx - comp*Rational(1,3)*divW) \n", + "EzzW = inW *((vNW-vSW)/dz - comp*Rational(1,3)*divW)\n", + "ExxE = inE *((uE -uC )/dx - comp*Rational(1,3)*divE)\n", + "EzzE = inE *((vNE-vSE)/dz - comp*Rational(1,3)*divE)\n", + "\n", + "# Additional missing stress tensor components needed for interpolation\n", + "GxzNE = inNEv*((uNE - uE )/dz + ( vNEE - vNE )/dx)\n", + "GxzSE = inSEv*((uE - uSE)/dz + ( vSEE - vSE )/dx)\n", + "GxzNW = inNWv*((uNW - uW )/dz + ( vNW - vNWW )/dx)\n", + "GxzSW = inSWv*((uW - uSW)/dz + ( vSW - vSWW )/dx)\n", + "\n", + "# Shear strain rate\n", + "dVxdzN = RegN*dVxdzN_Reg + DirN*dVxdzN_Dir + NeuN*dVxdzN_Neu\n", + "dVxdzS = RegS*dVxdzS_Reg + DirS*dVxdzS_Dir + NeuS*dVxdzS_Neu\n", + "\n", + "dVzdxN = SxxBCW*(DirNW*dVzdxN_DirW + NeuNW*dVzdxN_NeuW) + (1-SxxBCW-SxxBCE)*dVzdxN_Reg + SxxBCE*(DirNE*dVzdxN_DirE + NeuNE*dVzdxN_NeuE) \n", + "dVzdxS = SxxBCW*(DirSW*dVzdxS_DirW + NeuSW*dVzdxS_NeuW) + (1-SxxBCW-SxxBCE)*dVzdxS_Reg + SxxBCE*(DirSE*dVzdxS_DirE + NeuSE*dVzdxS_NeuE)\n", + "\n", + "GxzN = inN*(dVxdzN + dVzdxN) \n", + "GxzS = inS*(dVxdzS + dVzdxS) \n", + "ExxNE = inNEc*((uNE -uN )/dx - comp*Rational(1,3)*divNE)\n", + "ExxNW = inNWc*((uN -uNW)/dx - comp*Rational(1,3)*divNW)\n", + "ExxSE = inSEc*((uSE -uS )/dx - comp*Rational(1,3)*divSE)\n", + "ExxSW = inSWc*((uS -uSW)/dx - comp*Rational(1,3)*divSW)\n", + "EzzNE = inNEc*((vNNE-vNE)/dz - comp*Rational(1,3)*divNE)\n", + "EzzNW = inNWc*((vNNW-vNW)/dz - comp*Rational(1,3)*divNW)\n", + "EzzSE = inSEc*((vSE-vSSE)/dz - comp*Rational(1,3)*divSE)\n", + "EzzSW = inSWc*((vSW-vSSW)/dz - comp*Rational(1,3)*divSW)\n", + "# INV 0\n", + "GxzE = 0.25*( inN*GxzN + inS*GxzS + inNEv*GxzNE + inSEv*GxzSE) # clear\n", + "GxzW = 0.25*( inN*GxzN + inS*GxzS + inNWv*GxzNW + inSWv*GxzSW)\n", + "ExxN = 0.25*( inE*ExxE + inW*ExxW + inNWc*ExxNW + inNEc*ExxNE)\n", + "ExxS = 0.25*( inE*ExxE + inW*ExxW + inSWc*ExxSW + inSEc*ExxSE)\n", + "EzzN = 0.25*( inE*EzzE + inW*EzzW + inNWc*EzzNW + inNEc*EzzNE) \n", + "EzzS = 0.25*( inE*EzzE + inW*EzzW + inSWc*EzzSW + inSEc*EzzSE)\n", + "pS = 0.25*( inE*pE + inW*pW + inSWc*pSW + inSEc*pSE)\n", + "pN = 0.25*( inE*pE + inW*pW + inSWc*pNW + inSEc*pNE) \n", + "\n", + "#----------------------------------------#\n", + "EdW = Matrix([[ExxW], [EzzW], [GxzW], [pW]])\n", + "TW = Dv*EdW\n", + "TxxW = TW[0].subs({D11:'D11W', D12:'D12W', D13:'D13W', D14:'D14W'})\n", + "#----------------------------------------#\n", + "EdE = Matrix([[ExxE], [EzzE], [GxzE], [pE]])\n", + "TE = Dv*EdE\n", + "TxxE = TE[0].subs({D11:'D11E', D12:'D12E', D13:'D13E', D14:'D14E'})\n", + "#----------------------------------------#\n", + "EdS = Matrix([[ExxS], [EzzS], [GxzS], [pS]])\n", + "TS = Dv*EdS\n", + "TxzS = TS[2].subs({D31:'D31S', D32:'D32S', D33:'D33S', D34:'D34S'})\n", + "#----------------------------------------#\n", + "EdN = Matrix([[ExxN], [EzzN], [GxzN], [pN]])\n", + "TN = Dv*EdN\n", + "TxzN = TN[2].subs({D31:'D31N', D32:'D32N', D33:'D33N', D34:'D34N'})\n", + "#----------------------------------------#\n", + "Fx_Reg = 1/dx*(TxxE - TxxW) + 1/dz*(TxzN - TxzS) - 1/dx*(pE - pW)\n", + "Fx_Reg = -Fx_Reg\n", + "\n", + "# Stress boundary\n", + "SxxE = TxxE - pE\n", + "SxxW = 2*SxxBC - SxxE\n", + "Fx_SxxBCW = 1/dx*(SxxE - SxxW) + 1/dz*(TxzN - TxzS)\n", + "Fx_SxxBCW = -Fx_SxxBCW/2\n", + "\n", + "SxxW = TxxW - pW\n", + "SxxE = 2*SxxBC - SxxW\n", + "Fx_SxxBCE = 1/dx*(SxxE - SxxW) + 1/dz*(TxzN - TxzS)\n", + "Fx_SxxBCE = -Fx_SxxBCE/2\n", + "\n", + "Fx = (1-SxxBCW-SxxBCE)*Fx_Reg + SxxBCW*Fx_SxxBCW + SxxBCE*Fx_SxxBCE\n", + "\n", + "for i in range(len(dofs)):\n", + " cUc = Fx.diff(dofs[i])\n", + " print(str(dofs[i]) + ' = ' + ccode(simplify(cUc)) + ';')" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "vW = (1.0/2.0)*D33W*RegW*inW*(SzzBCN + SzzBCS - 2)/pow(dx, 2);\n", + "vC = (1.0/6.0)*(-SzzBCN*(2*D22S*pow(dx, 2)*(OOP*comp - 3) - 3*pow(dz, 2)*(D33E*inE*(2*DirE + RegE) + D33W*inW*(2*DirW + RegW))) - SzzBCS*(2*D22N*pow(dx, 2)*(OOP*comp - 3) - 3*pow(dz, 2)*(D33E*inE*(2*DirE + RegE) + D33W*inW*(2*DirW + RegW))) + 2*(pow(dx, 2)*(D22N + D22S)*(OOP*comp - 3) - 3*pow(dz, 2)*(D33E*inE*(2*DirE + RegE) + D33W*inW*(2*DirW + RegW)))*(SzzBCN + SzzBCS - 1))/(pow(dx, 2)*pow(dz, 2));\n", + "vE = (1.0/2.0)*D33E*RegE*inE*(SzzBCN + SzzBCS - 2)/pow(dx, 2);\n", + "vS = (1.0/3.0)*D22S*(1 - SzzBCS)*(OOP*comp - 3)/pow(dz, 2);\n", + "vN = (1.0/3.0)*D22N*(1 - SzzBCN)*(OOP*comp - 3)/pow(dz, 2);\n", + "uSW = (1.0/6.0)*(-3*D33W*SzzBCS*inW*(2*DirNW*SzzBCN - SzzBCN - SzzBCS + 1) + SzzBCN*(2*D22S*OOP*comp - 3*D33W*inW*(2*DirNW*SzzBCN - SzzBCN - SzzBCS + 1)) - 2*(D22S*OOP*comp - 3*D33W*inW*(2*DirNW*SzzBCN - SzzBCN - SzzBCS + 1))*(SzzBCN + SzzBCS - 1))/(dx*dz);\n", + "uSE = (1.0/6.0)*(3*D33E*SzzBCS*inE*(2*DirNE*SzzBCN - SzzBCN - SzzBCS + 1) - SzzBCN*(2*D22S*OOP*comp - 3*D33E*inE*(2*DirNE*SzzBCN - SzzBCN - SzzBCS + 1)) + 2*(D22S*OOP*comp - 3*D33E*inE*(2*DirNE*SzzBCN - SzzBCN - SzzBCS + 1))*(SzzBCN + SzzBCS - 1))/(dx*dz);\n", + "uNW = (1.0/6.0)*(3*D33W*SzzBCN*inW*(2*DirSW*SzzBCS - SzzBCN - SzzBCS + 1) - SzzBCS*(2*D22N*OOP*comp - 3*D33W*inW*(2*DirSW*SzzBCS - SzzBCN - SzzBCS + 1)) + 2*(D22N*OOP*comp - 3*D33W*inW*(2*DirSW*SzzBCS - SzzBCN - SzzBCS + 1))*(SzzBCN + SzzBCS - 1))/(dx*dz);\n", + "uNE = (1.0/6.0)*(-3*D33E*SzzBCN*inE*(2*DirSE*SzzBCS - SzzBCN - SzzBCS + 1) + SzzBCS*(2*D22N*OOP*comp - 3*D33E*inE*(2*DirSE*SzzBCS - SzzBCN - SzzBCS + 1)) - 2*(D22N*OOP*comp - 3*D33E*inE*(2*DirSE*SzzBCS - SzzBCN - SzzBCS + 1))*(SzzBCN + SzzBCS - 1))/(dx*dz);\n", + "vSW = 0;\n", + "vSE = 0;\n", + "vNW = 0;\n", + "vNE = 0;\n", + "uSWW = 0;\n", + "uSEE = 0;\n", + "uNWW = 0;\n", + "uNEE = 0;\n", + "uSSW = 0;\n", + "uSSE = 0;\n", + "uNNW = 0;\n", + "uNNE = 0;\n", + "pS = (SzzBCS - 1)/dz;\n", + "pN = (1 - SzzBCN)/dz;\n", + "pSW = 0;\n", + "pSE = 0;\n", + "pNW = 0;\n", + "pNE = 0;\n" + ] + } + ], + "source": [ + "# Compressibility switch\n", + "comp = symbols('comp')\n", + "\n", + "SzzBC = symbols('SzzBC') # boundary stress\n", + "SxzBC = symbols('SxzBC') # boundary stress\n", + "SzzBCS = symbols('SzzBCS') # Flag to switch between BC and internal nodes\n", + "SzzBCN = symbols('SzzBCN') # Flag to switch between BC and internal nodes\n", + "SxzBCW = symbols('SxzBCW') # Flag to switch between BC and internal nodes\n", + "SxzBCE = symbols('SxzBCE') # Flag to switch between BC and internal nodes\n", + "\n", + "# Stencil z\n", + "vW,vC,vE,vS,vN,uSW,uSE,uNW,uNE = symbols('vW,vC,vE,vS,vN,uSW,uSE,uNW,uNE')\n", + "dofs = Matrix([vW,vC,vE,vS,vN,uSW,uSE,uNW,uNE])\n", + "\n", + "# Stencil extension 1\n", + "vSW,vSE,vNW,vNE,uSWW,uSEE,uNWW,uNEE,uSSW,uSSE,uNNW,uNNE = symbols('vSW,vSE,vNW,vNE,uSWW,uSEE,uNWW,uNEE,uSSW,uSSE,uNNW,uNNE')\n", + "dofs_ext = Matrix([vSW,vSE,vNW,vNE,uSWW,uSEE,uNWW,uNEE,uSSW,uSSE,uNNW,uNNE])\n", + "dofs = np.append(dofs,dofs_ext)\n", + "\n", + "# Stencil extension 2\n", + "pS, pN, pSW, pSE, pNW, pNE = symbols('pS, pN, pSW, pSE, pNW, pNE')\n", + "dofs_ext = Matrix([pS, pN, pSW, pSE, pNW, pNE])\n", + "dofs = np.append(dofs,dofs_ext)\n", + "\n", + "# Flags\n", + "inS,inN,inW,inE = symbols('inS,inN,inW,inE')\n", + "inSEc,inSWc,inNEc,inNWc = symbols('inSEc,inSWc,inNEc,inNWc')\n", + "inSEv,inSWv,inNEv,inNWv = symbols('inSEv,inSWv,inNEv,inNWv')\n", + "DirE, DirW, NeuE, NeuW, RegE, RegW = symbols('DirE, DirW, NeuE, NeuW, RegE, RegW')\n", + "\n", + "# East/West dVzdx\n", + "dVzdxE_Reg = (vE - vC )/dx\n", + "dVzdxW_Reg = (vC - vW )/dx\n", + "dVzdxE_Dir = (-vC - vC )/dx\n", + "dVzdxW_Dir = (vC - -vC )/dx\n", + "dVzdxE_Neu = (vC - vC )/dx\n", + "dVzdxW_Neu = (vC - vC )/dx\n", + "\n", + "# East/West dVxdz\n", + "dVxdzE_Reg = (uNE-uSE)/dz\n", + "dVxdzW_Reg = (uNW-uSW)/dz\n", + "dVxdzE_DirS = (uNE-(-uNE))/dz\n", + "dVxdzE_NeuS = 0.0\n", + "dVxdzE_DirN = ((-uSE)-uSE)/dz\n", + "dVxdzE_NeuN = 0.0\n", + "dVxdzW_DirS = (uNW-(-uNW))/dz\n", + "dVxdzW_NeuS = 0.0\n", + "dVxdzW_DirN = ((-uSW)-uSW)/dz\n", + "dVxdzW_NeuN = 0.0\n", + "\n", + "# Divergences\n", + "divS = out_of_plane*( (uSE-uSW)/dx + (vC -vS )/dz ) \n", + "divN = out_of_plane*( (uNE-uNW)/dx + (vN -vC )/dz ) \n", + "divSW = out_of_plane*( (uSW-uSWW)/dx + (vW -vSW )/dz )\n", + "divSE = out_of_plane*( (uSEE-uSE)/dx + (vE -vSE )/dz )\n", + "divNW = out_of_plane*( (uNW-uNWW)/dx + (vNW -vW )/dz )\n", + "divNE = out_of_plane*( (uNEE-uNE)/dx + (vNE -vE )/dz ) \n", + "\n", + "# Deviatoric normal strain\n", + "ExxS = ((uSE-uSW)/dx - comp*Rational(1,3)*divS)\n", + "EzzS = ((vC -vS )/dz - comp*Rational(1,3)*divS)\n", + "ExxN = ((uNE-uNW)/dx - comp*Rational(1,3)*divN)\n", + "EzzN = ((vN -vC )/dz - comp*Rational(1,3)*divN)\n", + "\n", + "# Additional missing stress tensor components via interpolation\n", + "GxzSW = ((uSW - uSSW)/dz + (vS - vSW)/dx)\n", + "GxzSE = ((uSE - uSSE)/dz + (vSE - vS )/dx)\n", + "GxzNW = ((uNNW - uNW )/dz + (vN - vNW)/dx)\n", + "GxzNE = ((uNNE - uNE )/dz + (vNE - vN )/dx)\n", + "\n", + "# Shear strain rate\n", + "# GxzE = inE*((uNE-uSE)/dz + RegE*dVzdxE_Reg + DirE*dVzdxE_Dir + NeuE*dVzdxE_Neu)\n", + "# GxzW = inW*((uNW-uSW)/dz + RegW*dVzdxW_Reg + DirW*dVzdxW_Dir + NeuW*dVzdxW_Neu)\n", + "\n", + "dVzdxE = RegE*dVzdxE_Reg + DirE*dVzdxE_Dir + NeuE*dVzdxE_Neu\n", + "dVzdxW = RegW*dVzdxW_Reg + DirW*dVzdxW_Dir + NeuW*dVzdxW_Neu\n", + "\n", + "dVxdzE = SzzBCS*(DirSE*dVxdzE_DirS + NeuSE*dVxdzE_NeuS) + (1-SzzBCS-SzzBCN)*dVxdzE_Reg + SzzBCN*(DirNE*dVxdzE_DirN + NeuNE*dVxdzE_NeuN) \n", + "dVxdzW = SzzBCS*(DirSW*dVxdzW_DirS + NeuSW*dVxdzW_NeuS) + (1-SzzBCS-SzzBCN)*dVxdzW_Reg + SzzBCN*(DirNW*dVxdzW_DirN + NeuNW*dVxdzW_NeuN)\n", + "\n", + "GxzE = inE*(dVxdzE + dVzdxE) \n", + "GxzW = inW*(dVxdzW + dVzdxW) \n", + "\n", + "# Additional missing stress tensor components for interpolation\n", + "ExxNE = ((uNEE-uNE)/dx - comp*1/3*(divNE))\n", + "ExxNW = ((uNW-uNWW)/dx - comp*1/3*(divNW))\n", + "ExxSE = ((uSEE-uSE)/dx - comp*1/3*(divSE))\n", + "ExxSW = ((uSW-uSWW)/dx - comp*1/3*(divSW))\n", + "EzzNE = ((vNE -vE )/dz - comp*1/3*(divNE))\n", + "EzzNW = ((vNW -vW )/dz - comp*1/3*(divNW))\n", + "EzzSE = ((vE -vSE )/dz - comp*1/3*(divSE))\n", + "EzzSW = ((vW -vSW )/dz - comp*1/3*(divSW))\n", + "\n", + "# INV 0\n", + "GxzN = Rational(1,4)*( inW*GxzW + inE*GxzE + inNWv*GxzNW + inNEv*GxzNE )\n", + "GxzS = Rational(1,4)*( inW*GxzW + inE*GxzE + inSWv*GxzSW + inSEv*GxzSE ) \n", + "ExxE = Rational(1,4)*( inS*ExxS + inN*ExxN + inNEc*ExxNE + inSEc*ExxSE )\n", + "ExxW = Rational(1,4)*( inS*ExxS + inN*ExxN + inNWc*ExxNW + inSWc*ExxSW )\n", + "EzzE = Rational(1,4)*( inS*EzzS + inN*EzzN + inNEc*EzzNE + inSEc*EzzSE )\n", + "EzzW = Rational(1,4)*( inS*EzzS + inN*EzzN + inNWc*EzzNW + inSWc*EzzSW )\n", + "pE = Rational(1,4)*( inS*pS + inN*pN + inNEc*pNE + inSEc*pSE )\n", + "pW = Rational(1,4)*( inS*pS + inN*pN + inNWc*pNW + inSWc*pSW )\n", + "\n", + "#----------------------------------------#\n", + "EdS = Matrix([[ExxS], [EzzS], [GxzS], [pS]])\n", + "TS = Dv*EdS\n", + "TzzS = TS[1].subs({D21:'D21S', D22:'D22S', D23:'D23S', D24:'D24S'})\n", + "#----------------------------------------#\n", + "EdN = Matrix([[ExxN], [EzzN], [GxzN], [pN]])\n", + "TN = Dv*EdN\n", + "TzzN = TN[1].subs({D21:'D21N', D22:'D22N', D23:'D23N', D24:'D24N'})\n", + "#----------------------------------------#\n", + "EdW = Matrix([[ExxW], [EzzW], [GxzW], [pW]])\n", + "TW = Dv*EdW\n", + "TxzW = TW[2].subs({D31:'D31W', D32:'D32W', D33:'D33W', D34:'D34W'})\n", + "#----------------------------------------#\n", + "EdE = Matrix([[ExxE], [EzzE], [GxzE], [pE]])\n", + "TE = Dv*EdE\n", + "TxzE = TE[2].subs({D31:'D31E', D32:'D32E', D33:'D33E', D34:'D34E'})\n", + "#----------------------------------------#\n", + "# Shear stress gradient\n", + "dTxzdx = 1/dx*(TxzE - TxzW)\n", + "\n", + "Fz_Reg = 1/dz*(TzzN - TzzS) + 1/dx*(TxzE - TxzW) - (pN - pS)/dz \n", + "Fz_Reg = -Fz_Reg\n", + "\n", + "SzzN = TzzN - pN\n", + "SzzS = 2*SzzBC - SzzN\n", + "Fz_SzzBCS = 1/dz*(SzzN - SzzS)\n", + "Fz_SzzBCS += dTxzdx \n", + "Fz_SzzBCS = -Fz_SzzBCS/2\n", + "\n", + "SzzS = TzzS - pS\n", + "SzzN = 2*SzzBC - SzzS\n", + "Fz_SzzBCN = 1/dz*(SzzN - SzzS)\n", + "Fz_SzzBCN += dTxzdx \n", + "Fz_SzzBCN = -Fz_SzzBCN/2\n", + "\n", + "Fz = (1-SzzBCS-SzzBCN)*Fz_Reg + SzzBCS*Fz_SzzBCS + SzzBCN*Fz_SzzBCN\n", + "\n", + "for i in range(len(dofs)):\n", + " cVc = Fz.diff(dofs[i])\n", + " print(str(dofs[i]) + ' = ' + ccode(simplify(cVc)) + ';') " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.9.5 ('base')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.6" + }, + "vscode": { + "interpreter": { + "hash": "309baa410e37084482bcba9b39a5b9e635e78b91cd7553c4b48a2d89780d0f88" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}