This answer restricts to groups under the multiplication operation modulo a large prime $p$, because the question does (others groups are increasingly common in cryptography, including Elliptic Curve Groups).

So we'll be working within the group $\mathbb Z_p^*$, a notation for the subset of the ring $\mathbb Z/p\mathbb Z$ formed by elements that have a multiplicative inverse, which can be shown to form a group. Since $p$ is prime, $\mathbb Z/p\mathbb Z$ is a field, and $\mathbb Z_p^*$ is that field less the neutral for addition. $\mathbb Z_p^*$ can be assimilated to integers $\{1,2,\ldots,p-2,p-1\}$. It has $p-1$ elements. This group's internal law is multiplication modulo $p$.

We want a group $\mathbb Z_p^*$ (or a subgroup thereof) in which computing the discrete logarithm is hard. Conditions for this (both provably necessary, and conjectured sufficient when both are met):
- The order of the group (that is, it's number of elements) must have a large-enough prime factor $q$ to block [Pohlig-Hellman][1] combined with [Pollard's rho][1], which can compute the logarithm with effort $\mathcal O\,\left(\sqrt q(\ln n)^3\right)$ and moderate memory, once $p$ and $q$ are known. Thus we want $q$ at least 256-bit of 128-bit security.
- $p$ must be large enough to block algorithms of the [index calculus familly][2]. That would require $p$ of 2048 to 4096 bit (depending on conservatism, hypothesis and source) for 128-bit security, further assuming $p$ is not too close to the power of a small integer. 

> How does practical Discrete Logarithm systems get set up?

First thing is to decide what kind of discrete logarithm framework we want to setup, and that depends on what we plan to use it for. There are three kinds commonly used:
1. The full group $\mathbb Z_p^*$, with order $2q=p-1$ and $q$ prime. That one is used when something mandates a generator of the full group (as the question's "generator of multiplicative group $\mathbb Z/p\mathbb Z$" taken to the letter does).
2. The subgroup of quadratic residues, with prime order $q=(p-1)/2$. It's elements are the $y$ in $\mathbb Z_p^*$ that can be written as $y=u^2$ for some element $u$ of $\mathbb Z_p^*$ (and then also for $u'=p-u$, of course). There are standard such groups with generator $g=2$, see [RFC 3526][3]. They can be used for most uses that do not require a generator of the full group.
3. A smaller subset of the subset of quadratic residues, with prime order $q=(p-1)/r$ for some (even) $r$. These are called Schnorr group. They are used in the original Schnorr signature, and DSA.

There's often a reason to not use (1.): if $g$ is a generator (that is, any element of the group can be expressed as $g^x$ ), then from the value of $g^x$ it's easy to tell if $x$ is even or odd (by computing the [Legendre symbol][4]). That goes against hypothesis desirable for simple security proofs, against objectives in zero-knowledge proofs, and might allow attacks (e.g. ElGamal encryption fails [CPA][5]).

There's sometime a reason to not use (2.) and prefer (3.): in signature applications at least, one of the component of the signature has the same bit size as $q$, which needs to be large but can be much smaller than $p$. Thus using a small $q$ allows shorter signatures. That's the reason Schnorr groups have been introduced.
___
For (1.) and (2.), generating $p$ and $q$ boils down to generating prime $q$ with $p=2q+1$ also prime. It's efficient to make a relatively quick test that $q$ is not composite (e.g. a Fermat test to base 2: check that $2^{q-1}\bmod q=1$), then a thorough test that $p=2q+1$ is prime; then a thorough test that $q$ is prime. Further, for a small prime $s$, it holds $q\bmod s\not\in\{0,(s-1)/2\}$. Thus $q\bmod3=2$, $q\bmod5\in\{1,3,4\}$, thus $q\bmod30\in\{11,23,29\}$, which narrows down the search. Considering slightly larger $s$ can be put to use for a sieve or other speedup.

Also, it can be shown that for $p$ and $q$ as for (1.) or (2.), $g=2$ is a generator for (1.) if and only if $q\bmod4=1$ (equivalently $p\bmod8=3$ ); and a generator for (2.) otherwise. Thus to use $g=2$ (the smallest possible generator, and one that slightly simplifies some computations), we want $q\bmod60\in\{29,41,53\}$ for (1.) and $q\bmod60\in\{11,23,59\}$ for (2.).

For (3.), we can pick random prime $q$ of the desired size, then random even $r$ of adequate size to have $p=q\,r+1$ prime. A generator is $g=2^r$ (unless that's $1$, which is at least overwhelmingly unlikely; I wonder if it's even possible):

If we want a random generator, we can pick a random secret $x$ in $[1,q]$ and use $g$ as follows, after checking $g\ne1$ ($g=1$ is overwhelmingly unlikely):
- for (1.) and $p\bmod8=3$: $g\gets2^{2x+1}$
- for (1.) and $p\bmod8=7$: by trial and error we find $y$ with $y^{(p-1)/2}\bmod p\ne 1$ (equivalently, with $y^{(p-1)/2}\bmod p=p-1$ ); the expected number of attempts is about two if we explore consecutive odd $u$ starting $u=3$; then $g\gets u^{2x+1}$.
- for (2.), $g\gets4^x$
- for (3.), $g\gets2^{(p/q)\,x}$

> There is no known algorithm to find $g$ in polynomial time.

That's true if the factorization of the group order is unknown, or if we restrict to _deterministic_ algorithms. But in practice the order is known, sometime we can predict a generator, and simple probabilistic algorithms will quickly yield one otherwise.

> Is there any big number packages?

Python has native support for big numbers, and it's a simple exercise to build the parameters $(p,g)$ in pure Python. The only mild difficulties are the prime test and the random prime generation. But they can be done in pure Python, and there's support for these in several packages including SymPy. Since $(p,g)$ are usually public, side channels are not an issue.


  [1]: https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm_for_logarithms
  [2]: https://en.wikipedia.org/wiki/Index_calculus_algorithm#The_Index_Calculus_family
  [3]: https://datatracker.ietf.org/doc/html/rfc3526
  [4]: https://en.wikipedia.org/wiki/Legendre_symbol
  [5]: https://en.wikipedia.org/wiki/Chosen-plaintext_attack