Lines Matching refs:param
103 Parameter *param; member
127 Parameter param; in main() local
147 PetscCall(SetParams(¶m, &grid)); in main()
148 PetscCall(ReportParams(¶m, &grid)); in main()
166 user->param = ¶m; in main()
214 Parameter *param = user->param; in UpdateSolution() local
230 if (param->ivisc == VISC_CONST && !param->stop_solve) { in UpdateSolution()
231 param->ivisc = VISC_CONST; in UpdateSolution()
238 if (param->stop_solve) goto done; in UpdateSolution()
242 if (param->ivisc >= VISC_DIFN && !param->stop_solve) { in UpdateSolution()
246 for (param->continuation = 0.0;; param->continuation += cont_incr) { in UpdateSolution()
247 …Call(PetscPrintf(PETSC_COMM_WORLD, " Continuation parameter = %g\n", (double)param->continuation)); in UpdateSolution()
256 if (param->stop_solve) goto done; in UpdateSolution()
266 if (param->continuation >= 1.0) goto done; in UpdateSolution()
271 if (param->continuation + cont_incr > 1.0) cont_incr = 1.0 - param->continuation; in UpdateSolution()
276 …if (param->stop_solve && !q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: stopping solve.… in UpdateSolution()
308 Parameter *param = user->param; in HorizVelocity() local
310 PetscScalar st, ct, th, c = param->c, d = param->d; in HorizVelocity()
325 Parameter *param = user->param; in VertVelocity() local
327 PetscScalar st, ct, th, c = param->c, d = param->d; in VertVelocity()
342 Parameter *param = user->param; in Pressure() local
344 PetscScalar x, z, r, st, ct, c = param->c, d = param->d; in Pressure()
357 Parameter *param = user->param; in CalcSecInv() local
376 uN = param->cb; in CalcSecInv()
377 wW = param->sb; in CalcSecInv()
394 wN = param->sb; in CalcSecInv()
395 uW = param->cb; in CalcSecInv()
420 static inline PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PetscScalar z, Parameter *param) in Viscosity() argument
423 ViscParam difn = param->diffusion, disl = param->dislocation; in Viscosity()
424 PetscInt iVisc = param->ivisc; in Viscosity()
425 PetscScalar eps_scale = param->V / (param->L * 1000.0); in Viscosity()
429 P = rho_g * (z * param->L * 1000.0); /* Pa */ in Viscosity()
436 …scRealPart(difn.A * PetscExpScalar((difn.Estar + P * difn.Vstar) / R / (T + 273.0)) / param->eta0); in Viscosity()
441 …cExpScalar((disl.Estar + P * disl.Vstar) / disl.n / R / (T + 273.0)) * strain_power / param->eta0); in Viscosity()
446 v1 = difn.A * PetscExpScalar((difn.Estar + P * difn.Vstar) / R / (T + 273.0)) / param->eta0; in Viscosity()
447 …scExpScalar((disl.Estar + P * disl.Vstar) / disl.n / R / (T + 273.0)) * strain_power / param->eta0; in Viscosity()
455 result = PetscMax(result, param->visc_cutoff); in Viscosity()
457 result = PetscPowReal(result, param->continuation); in Viscosity()
464 Parameter *param = user->param; in XMomentumResidual() local
472 z_scale = param->z_scale; in XMomentumResidual()
474 if (param->ivisc == VISC_DIFN || param->ivisc >= VISC_DISL) { /* viscosity is T-dependent */ in XMomentumResidual()
475 TS = param->potentialT * TInterp(x, i, j - 1) * PetscExpScalar((j - 1.0) * dz * z_scale); in XMomentumResidual()
477 else TN = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale); in XMomentumResidual()
478 TW = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale); in XMomentumResidual()
479 TE = param->potentialT * x[j][i + 1].T * PetscExpScalar((j - 0.5) * dz * z_scale); in XMomentumResidual()
480 if (param->ivisc >= VISC_DISL) { /* olivine dislocation creep */ in XMomentumResidual()
487 etaN = Viscosity(TN, epsN, dz * (j + 0.5), param); in XMomentumResidual()
488 etaS = Viscosity(TS, epsS, dz * (j - 0.5), param); in XMomentumResidual()
489 etaW = Viscosity(TW, epsW, dz * j, param); in XMomentumResidual()
490 etaE = Viscosity(TE, epsE, dz * j, param); in XMomentumResidual()
502 if (param->ivisc != VISC_CONST) { in XMomentumResidual()
515 Parameter *param = user->param; in ZMomentumResidual() local
524 z_scale = param->z_scale; in ZMomentumResidual()
527 if (param->ivisc == VISC_DIFN || param->ivisc >= VISC_DISL) { /* viscosity is T-dependent */ in ZMomentumResidual()
528 TN = param->potentialT * x[j + 1][i].T * PetscExpScalar((j + 0.5) * dz * z_scale); in ZMomentumResidual()
529 TS = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale); in ZMomentumResidual()
530 TW = param->potentialT * TInterp(x, i - 1, j) * PetscExpScalar(j * dz * z_scale); in ZMomentumResidual()
532 else TE = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale); in ZMomentumResidual()
533 if (param->ivisc >= VISC_DISL) { /* olivine dislocation creep */ in ZMomentumResidual()
540 etaN = Viscosity(TN, epsN, dz * (j + 1.0), param); in ZMomentumResidual()
541 etaS = Viscosity(TS, epsS, dz * (j + 0.0), param); in ZMomentumResidual()
542 etaW = Viscosity(TW, epsW, dz * (j + 0.5), param); in ZMomentumResidual()
543 etaE = Viscosity(TE, epsE, dz * (j + 0.5), param); in ZMomentumResidual()
556 if (param->ivisc != VISC_CONST) { in ZMomentumResidual()
585 Parameter *param = user->param; in EnergyResidual() local
600 dx * dz / param->peclet; in EnergyResidual()
607 uW = uE = param->cb; in EnergyResidual()
608 wS = wN = param->sb; in EnergyResidual()
617 if (param->adv_scheme == ADVECT_FV || i == ilim - 1 || j == jlim - 1 || i == 1 || j == 1) { in EnergyResidual()
644 Parameter *param = user->param; in ShearStress() local
655 wW = param->sb; in ShearStress()
656 uN = param->cb; in ShearStress()
678 Parameter *param = user->param; in XNormalStress() local
685 ivisc = param->ivisc; in XNormalStress()
686 z_scale = param->z_scale; in XNormalStress()
690 TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale); in XNormalStress()
692 etaC = Viscosity(TC, epsC, dz * j, param); in XNormalStress()
701 TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale); in XNormalStress()
703 etaC = Viscosity(TC, epsC, dz * (j + 0.5), param); in XNormalStress()
705 if (i == j) uW = param->sb; in XNormalStress()
717 Parameter *param = user->param; in ZNormalStress() local
725 ivisc = param->ivisc; in ZNormalStress()
726 z_scale = param->z_scale; in ZNormalStress()
730 TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale); in ZNormalStress()
732 etaC = Viscosity(TC, epsC, dz * j, param); in ZNormalStress()
740 TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale); in ZNormalStress()
742 etaC = Viscosity(TC, epsC, dz * (j + 0.5), param); in ZNormalStress()
743 if (i == j) wN = param->sb; in ZNormalStress()
758 PetscErrorCode SetParams(Parameter *param, GridInfo *grid) in SetParams() argument
765 param->slab_dip = 45.0; in SetParams()
766 param->width = 320.0; /* km */ in SetParams()
767 param->depth = 300.0; /* km */ in SetParams()
768 param->lid_depth = 35.0; /* km */ in SetParams()
769 param->fault_depth = 35.0; /* km */ in SetParams()
771 PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_dip", ¶m->slab_dip, NULL)); in SetParams()
772 PetscCall(PetscOptionsGetReal(NULL, NULL, "-width", ¶m->width, NULL)); in SetParams()
773 PetscCall(PetscOptionsGetReal(NULL, NULL, "-depth", ¶m->depth, NULL)); in SetParams()
774 PetscCall(PetscOptionsGetReal(NULL, NULL, "-lid_depth", ¶m->lid_depth, NULL)); in SetParams()
775 PetscCall(PetscOptionsGetReal(NULL, NULL, "-fault_depth", ¶m->fault_depth, NULL)); in SetParams()
777 param->slab_dip = param->slab_dip * PETSC_PI / 180.0; /* radians */ in SetParams()
784 grid->dx = param->width / ((PetscReal)(grid->ni - 2)); /* km */ in SetParams()
785 grid->dz = grid->dx * PetscTanReal(param->slab_dip); /* km */ in SetParams()
786 grid->nj = (PetscInt)(param->depth / grid->dz + 3.0); /* gridpoints*/ in SetParams()
787 param->depth = grid->dz * (grid->nj - 2); /* km */ in SetParams()
798 param->pv_analytic = PETSC_FALSE; in SetParams()
799 param->ibound = BC_NOSTRESS; in SetParams()
800 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ibound", ¶m->ibound, NULL)); in SetParams()
803 param->slab_velocity = 5.0; /* cm/yr */ in SetParams()
804 param->slab_age = 50.0; /* Ma */ in SetParams()
805 param->lid_age = 50.0; /* Ma */ in SetParams()
806 param->kappa = 0.7272e-6; /* m^2/sec */ in SetParams()
807 param->potentialT = 1300.0; /* degrees C */ in SetParams()
809 PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_velocity", ¶m->slab_velocity, NULL)); in SetParams()
810 PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_age", ¶m->slab_age, NULL)); in SetParams()
811 PetscCall(PetscOptionsGetReal(NULL, NULL, "-lid_age", ¶m->lid_age, NULL)); in SetParams()
812 PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", ¶m->kappa, NULL)); in SetParams()
813 PetscCall(PetscOptionsGetReal(NULL, NULL, "-potentialT", ¶m->potentialT, NULL)); in SetParams()
816 param->ivisc = 3; /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */ in SetParams()
817 param->eta0 = 1e24; /* Pa-s */ in SetParams()
818 param->visc_cutoff = 0.0; /* factor of eta_0 */ in SetParams()
819 param->continuation = 1.0; in SetParams()
822 param->diffusion.A = 1.8e7; /* Pa-s */ in SetParams()
823 param->diffusion.n = 1.0; /* dim'less */ in SetParams()
824 param->diffusion.Estar = 375e3; /* J/mol */ in SetParams()
825 param->diffusion.Vstar = 5e-6; /* m^3/mol */ in SetParams()
828 param->dislocation.A = 2.8969e4; /* Pa-s */ in SetParams()
829 param->dislocation.n = 3.5; /* dim'less */ in SetParams()
830 param->dislocation.Estar = 530e3; /* J/mol */ in SetParams()
831 param->dislocation.Vstar = 14e-6; /* m^3/mol */ in SetParams()
833 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ivisc", ¶m->ivisc, NULL)); in SetParams()
834 PetscCall(PetscOptionsGetReal(NULL, NULL, "-visc_cutoff", ¶m->visc_cutoff, NULL)); in SetParams()
836 param->output_ivisc = param->ivisc; in SetParams()
838 PetscCall(PetscOptionsGetInt(NULL, NULL, "-output_ivisc", ¶m->output_ivisc, NULL)); in SetParams()
839 PetscCall(PetscOptionsGetReal(NULL, NULL, "-vstar", ¶m->dislocation.Vstar, NULL)); in SetParams()
842 param->quiet = PETSC_FALSE; in SetParams()
843 param->param_test = PETSC_FALSE; in SetParams()
845 PetscCall(PetscOptionsHasName(NULL, NULL, "-quiet", ¶m->quiet)); in SetParams()
846 PetscCall(PetscOptionsHasName(NULL, NULL, "-test", ¶m->param_test)); in SetParams()
847 …etscCall(PetscOptionsGetString(NULL, NULL, "-file", param->filename, sizeof(param->filename), &par… in SetParams()
850 param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */ in SetParams()
852 PetscCall(PetscOptionsGetInt(NULL, NULL, "-adv_scheme", ¶m->adv_scheme, NULL)); in SetParams()
855 param->stop_solve = PETSC_FALSE; in SetParams()
856 param->interrupted = PETSC_FALSE; in SetParams()
857 param->kspmon = PETSC_FALSE; in SetParams()
858 param->toggle_kspmon = PETSC_FALSE; in SetParams()
861 param->sb = PetscSinReal(param->slab_dip); in SetParams()
862 param->cb = PetscCosReal(param->slab_dip); in SetParams()
863 …param->c = param->slab_dip * param->sb / (param->slab_dip * param->slab_dip - param->sb * param->… in SetParams()
864 …param->d = (param->slab_dip * param->cb - param->sb) / (param->slab_dip * param->slab_dip - param… in SetParams()
867 param->L = PetscMin(param->width, param->depth); /* km */ in SetParams()
868 param->V = param->slab_velocity / 100.0 / SEC_PER_YR; /* m/sec */ in SetParams()
871 param->scaled_width = param->width / param->L; /* dim'less */ in SetParams()
872 param->scaled_depth = param->depth / param->L; /* dim'less */ in SetParams()
873 param->lid_depth = param->lid_depth / param->L; /* dim'less */ in SetParams()
874 param->fault_depth = param->fault_depth / param->L; /* dim'less */ in SetParams()
875 grid->dx = grid->dx / param->L; /* dim'less */ in SetParams()
876 grid->dz = grid->dz / param->L; /* dim'less */ in SetParams()
877 grid->jlid = (PetscInt)(param->lid_depth / grid->dz); /* gridcells */ in SetParams()
878 grid->jfault = (PetscInt)(param->fault_depth / grid->dz); /* gridcells */ in SetParams()
879 param->lid_depth = grid->jlid * grid->dz; /* dim'less */ in SetParams()
880 param->fault_depth = grid->jfault * grid->dz; /* dim'less */ in SetParams()
882 param->peclet = param->V /* m/sec */ in SetParams()
883 * param->L * 1000.0 /* m */ in SetParams()
884 / param->kappa; /* m^2/sec */ in SetParams()
885 param->z_scale = param->L * alpha_g_on_cp_units_inverse_km; in SetParams()
886 param->skt = PetscSqrtReal(param->kappa * param->slab_age * SEC_PER_YR); in SetParams()
887 PetscCall(PetscOptionsGetReal(NULL, NULL, "-peclet", ¶m->peclet, NULL)); in SetParams()
892 PetscErrorCode ReportParams(Parameter *param, GridInfo *grid) in ReportParams() argument
899 if (!param->quiet) { in ReportParams()
902 …OMM_WORLD, " Width = %g km, Depth = %g km\n", (double)param->width, (double)param->depth)… in ReportParams()
903 …%g degrees, Slab velocity = %g cm/yr\n", (double)(param->slab_dip * 180.0 / PETSC_PI), (double)pa… in ReportParams()
904 …f km, Fault depth = %5.2f km\n", (double)(param->lid_depth * param->L), (double)(param->fault_de… in ReportParams()
907 …] = %g, %g km\n", grid->ni, grid->nj, (double)(grid->dx * param->L), (double)(grid->dz * param->L)… in ReportParams()
909 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Pe = %g\n", (double)param->peclet)); in ReportParams()
912 if (param->ivisc == VISC_CONST) { in ReportParams()
914 …if (param->pv_analytic) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Pressur… in ReportParams()
915 } else if (param->ivisc == VISC_DIFN) { in ReportParams()
917 … Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * par… in ReportParams()
918 } else if (param->ivisc == VISC_DISL) { in ReportParams()
920 … Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * par… in ReportParams()
921 } else if (param->ivisc == VISC_FULL) { in ReportParams()
923 … Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * par… in ReportParams()
930 if (param->ibound == BC_ANALYTIC) { in ReportParams()
932 } else if (param->ibound == BC_NOSTRESS) { in ReportParams()
934 } else if (param->ibound == BC_EXPERMNT) { in ReportParams()
941 if (param->output_to_file) { in ReportParams()
943 …PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Output Destination: Mat file \"%s\"\n", param->file… in ReportParams()
945 …Printf(PETSC_COMM_WORLD, "Output Destination: PETSc binary file \"%s\"\n", param->filename)); in ReportParams()
948 …param->output_ivisc != param->ivisc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " … in ReportParams()
952 if (param->param_test) PetscCall(PetscEnd()); in ReportParams()
963 Parameter *param; in Initialize() local
973 param = user->param; in Initialize()
981 if (i < j) x[j][i].u = param->cb; in Initialize()
985 if (i <= j) x[j][i].w = param->sb; in Initialize()
1005 Parameter *param; in DoOutput() local
1017 param = user->param; in DoOutput()
1019 ivt = param->ivisc; in DoOutput()
1021 param->ivisc = param->output_ivisc; in DoOutput()
1031 if (param->output_to_file) { /* send output to binary file */ in DoOutput()
1040 PetscCall(VecSetValue(pars, 4, (PetscScalar)param->L, INSERT_VALUES)); in DoOutput()
1041 PetscCall(VecSetValue(pars, 5, (PetscScalar)param->V, INSERT_VALUES)); in DoOutput()
1043 PetscCall(VecSetValue(pars, 7, (PetscScalar)param->slab_dip, INSERT_VALUES)); in DoOutput()
1045 PetscCall(VecSetValue(pars, 9, (PetscScalar)param->lid_depth, INSERT_VALUES)); in DoOutput()
1047 PetscCall(VecSetValue(pars, 11, (PetscScalar)param->fault_depth, INSERT_VALUES)); in DoOutput()
1048 PetscCall(VecSetValue(pars, 12, (PetscScalar)param->potentialT, INSERT_VALUES)); in DoOutput()
1049 PetscCall(VecSetValue(pars, 13, (PetscScalar)param->ivisc, INSERT_VALUES)); in DoOutput()
1050 PetscCall(VecSetValue(pars, 14, (PetscScalar)param->visc_cutoff, INSERT_VALUES)); in DoOutput()
1051 PetscCall(VecSetValue(pars, 15, (PetscScalar)param->ibound, INSERT_VALUES)); in DoOutput()
1062 PetscCall(PetscViewerMatlabOpen(PETSC_COMM_WORLD, param->filename, FILE_MODE_WRITE, &viewer)); in DoOutput()
1064 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, param->filename, FILE_MODE_WRITE, &viewer)); in DoOutput()
1085 param->ivisc = ivt; in DoOutput()
1095 Parameter *param; in ViscosityField() local
1104 param = user->param; in ViscosityField()
1106 ivt = param->ivisc; in ViscosityField()
1107 param->ivisc = param->output_ivisc; in ViscosityField()
1125 …T = PetscRealPart(param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * param->z_scale)); in ViscosityField()
1127 …TC = PetscRealPart(param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * param->z_scale)); in ViscosityField()
1136 v[j][i].p = Viscosity(T, eps, dz * (j - 0.5), param); in ViscosityField()
1137 v[j][i].T = Viscosity(TC, epsC, dz * j, param); in ViscosityField()
1144 param->ivisc = ivt; in ViscosityField()
1194 Parameter *param = user->param; in SlabVel() local
1198 if (i < j - 1) return param->cb; in SlabVel()
1200 else return param->cb; in SlabVel()
1203 if (i < j) return param->sb; in SlabVel()
1205 else return param->sb; in SlabVel()
1212 Parameter *param = user->param; in PlateModel() local
1215 else z = (j - 0.5) * user->grid->dz * param->cb; /* PLATE_SLAB */ in PlateModel()
1217 return (PetscReal)(erf((double)PetscRealPart(z * param->L / 2.0 / param->skt))); in PlateModel()
1233 Parameter *param = user->param; in SNESConverged_Interactive() local
1237 if (param->interrupted) { in SNESConverged_Interactive()
1238 param->interrupted = PETSC_FALSE; in SNESConverged_Interactive()
1242 } else if (param->toggle_kspmon) { in SNESConverged_Interactive()
1243 param->toggle_kspmon = PETSC_FALSE; in SNESConverged_Interactive()
1247 if (param->kspmon) { in SNESConverged_Interactive()
1250 param->kspmon = PETSC_FALSE; in SNESConverged_Interactive()
1257 param->kspmon = PETSC_TRUE; in SNESConverged_Interactive()
1271 Parameter *param = user->param; in InteractiveHandler() local
1274 param->toggle_kspmon = PETSC_TRUE; in InteractiveHandler()
1277 param->interrupted = PETSC_TRUE; in InteractiveHandler()
1281 param->stop_solve = PETSC_TRUE; in InteractiveHandler()
1291 Parameter *param = user->param; in FormFunctionLocal() local
1309 /* ivisc = param->ivisc; */ ibound = param->ibound; in FormFunctionLocal()
1397 … - PetscPowRealInt(1.0 - PetscMax(PetscMin(PetscRealPart(x[j][i - 1].u) / param->cb, 1.0), 0.0), 5… in FormFunctionLocal()
1402 … - PetscPowRealInt(1.0 - PetscMax(PetscMin(PetscRealPart(x[j - 1][i].w) / param->sb, 1.0), 0.0), 5… in FormFunctionLocal()